%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                        %
%% CSE 557 - HW #1                                                        %
%% [http://www.cse.psu.edu/~plassman/cse557/assignments/hw1.html]         %
%% by Anirudh Modi (anirudh@anirudh.net) on 1/28/2000-Fri                 %
%% (modi@cse.psu.edu)                                                     %
%%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_grid = 20;  % number of different grid sizes 
errLimit = 2.0; % convergence limit specified as -log10(norm)
version = 1;    % which approach to be used (1 or 2?)
plotflag = 0;   % whether the data should be plotted?

filename = sprintf('conv%d.out',version);
outf = fopen(filename,'w');

for count=1:num_grid
  k = 5*count;          % k = 5,10,....,95,100
  % print out what we are doing.
  fprintf('Version %d for HW #1 with gridsize %dx%d and log10(resid) = %g\n',...
      version,k,k,-errLimit);
  
  % initialize two arrays, for the relaxation iteration...
  oldM = zeros(k+2,k+2); newM = oldM;
  
  % x and y positions for plotting
  X = [0:k+1]./(k+1);
  Y = [0:k+1]./(k+1);
  
  % set the middle of oldM to ones, keeping the boundary at zero.
  oldM(2:k+1,2:k+1) = ones(k,k);
  if version == 2
        newM = oldM;
  end
  
  % plot the initial image
  if plotflag == 1
    graph_title = sprintf('Initial conditions (%dx%d grid)',k,k);
    figure(1); subplot(3,1,1); surf(Y,X,oldM); title(graph_title); 
  end
  
  % do a few iterations iterations of relaxation (version 1):
  it = 0;
  errlog = 10.0;
  errlog0 = norm(oldM);
  while (-errlog < errLimit)
     it = it + 1;
     if version == 1
        for i=2:k+1
           for j=2:k+1
            newM(i,j) = 0.25*(oldM(i-1,j)+oldM(i+1,j)+oldM(i,j-1)+oldM(i,j+1));
          end
        end
        oldM(2:k+1,2:k+1) = newM(2:k+1,2:k+1);
     elseif version == 2
        for i=2:k+1
           for j=2:k+1
            newM(i,j) = 0.25*(newM(i-1,j)+newM(i+1,j)+newM(i,j-1)+newM(i,j+1));
          end
        end
     end
    % copy interior values back to oldM
    errlog = log10(norm(newM)/errlog0);
  end      % while loop for convergence ends

  fprintf('k = %d, iter = %d, errlog = %g\n',k,it,-errlog);
  fprintf(outf,'%d %d\n',k,it);
  
  % plot the image after "it" iterations
  if plotflag == 1
    graph_title = sprintf('After %d iterations (Version %d)',it,version);
    figure(1); subplot(3,1,2); surf(Y,X,newM); title(graph_title); pause(0); 
  end

end % the loop for grid size ends
fclose(outf);