% a simple function written to look at the eigenvalues
% of the weighted jacobian iteration matrix as they
% changed as a function of omega (the weight) in order
% to determine the smoothing factor.
%
% The smoothing factor is that value of the eigenvalue
% of the iteration matrix when the curve is such that
% the eigenvalue at N/2 is the same in absolute value as
% the eignvalue at N, where N is the number of eigenmodes.
%
% for 1-D, one can see that smoothing factor=1/3 and occurs
% when omega is 2/3
%
% by Nasser M. Abbasi
% math 228A, UC Davis, fal 2010
%
% Written using Matlab 2010a
%
% HOW TO RUN:
% put this function in the matlab path, and type its name below
% no input arguments needed.  This file requires one additional
% exterma; file: This requires the file my_xticklabels.m from
% matlab file exchange used to setting the xlabel ticks. This file
% is included in the same folder as this file. Make sure both m files
% are in the Matlab path.
%
% TYPE
%        nma_JacobiIterationMatrixEigenModes
%
% from the Matlab command line to start the program
%
% Move the slider to change omega. Omega is restriced between 0..1
%
% last updated 11/11/2010
%
function nma_JacobiIterationMatrixEigenModes()

    function my_Callback(hObject, eventdata, h)
        
        nPoints    = 100;  %fixed value. no need parametize
        omega      = get(h.omega,'Value');
        theta      = arrayfun(@(n) n*pi/nPoints,1:nPoints);
        eigenvalue = omega*cos(theta)+(1-omega);
        
        plot(1:nPoints,eigenvalue,'r-','LineWidth',2);  hold on;
        line([0 nPoints+5],[0 0]);
        
        % display values at half N and N to see where optimal weight is
        y1  = eigenvalue(nPoints/2);
        y2  = eigenvalue(nPoints);
        line([nPoints/2 nPoints/2],[0 y1],'LineWidth',2);
        line([nPoints nPoints],[0 y2],'LineWidth',2);
        text(nPoints/2,y1+0.1,sprintf('%3.3f',eigenvalue(nPoints/2)),...
            'HorizontalAlignment','left');
        text(nPoints,y2+0.1,sprintf('%3.3f',eigenvalue(nPoints)),...
            'HorizontalAlignment','left');
        
        %put circle where curve crosses x-axis
        if(min(eigenvalue)<0)
            [~,I] = min(abs(eigenvalue));
            plot(I,0,'ro','LineWidth',2);
            text(I,-0.2,sprintf('k=%d',I),...
                'HorizontalAlignment','left');
        end
        
        % set up plot
        xlim([0,nPoints+5]);    ylim([-1.2,1.2]);
        ylabel('eigenvalue');
        xlabel({'';'';'';'';'                select omega'},'Fontsize',8);
        
        set(gca,'XTick',0:20:nPoints);
        
        xtl = {{'k=1';'\theta=\pi/100'} {'20';'0.2 \pi'} {'40';'0.4 \pi'} ...
            {'60';'0.6 \pi'} {'80';'0.8 \pi'} {'100';'\pi'}};
        my_xticklabels(gca,[0 20 40 60 80 100],xtl,'FontSize',8);
        
        if (get(h.grid,'Value') == get(h.grid,'Max'))
            grid on;
        else
            grid off;
        end
        
        title(sprintf...
            ('Weighted Jacobian iteration matrix eigenvalue\n1-D problem. omega=%3.3f', ...
            omega));
        hold off;
    end

h.f = figure('Visible','on','Name','weighted Jacobian iteration matrix',...
    'Position',[250,250,600,550],'Resize','off');
h.hax = axes('parent',h.f,'Position',[.2,.25,.65,.65]);

h.omega = uicontrol(h.f,'Style','slider','Max',1,'Min',0.0,'Value',2/3,...
    'SliderStep',[0.002 0.02],'Units','normalized',...
    'Position',[.4 .035 .4 .05],'String','Omega');

h.grid=uicontrol(h.f,'Style','checkbox','Max',1,'Value',1,...
    'Units','normalized','Position',[.2 .035 .1 .05],...
    'String','grid on');

set(h.grid,'Callback',{@my_Callback h});
set(h.omega,'Callback',{@my_Callback h});
my_Callback(0,0,h);
end