% 
% A test and verification of the convergence properties (should be second order)
% of the routine hstcrt.f, using a known analytical solution to the 
% Cartesian Helmholtz equation.
% 
% Written by:
% -- 
% John L. Weatherwax                2006-07-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; 

A = 0; B = +2;
C = -1; D = +3;
ELMBDA = -4; 

nCells = 40*( 2.^[ 0:6 ] ); lInfty_error=zeros(1,length(nCells));
for gi=1:length(nCells),
  
  % Assign a uniform grid in both dimensions: 
  % 
  M=nCells(gi); N=nCells(gi); 

  X = linspace(A,B,M+1); 
  Y = linspace(C,D,N+1); 

  % Compute the required right hand side: 
  % 
  Fout = zeros(M+1,N+1); 
  for ii=1:M+1,
    x = A + (ii-1.0)*(B-A)/M;
    for jj=1:N+1,
      y = C + (jj-1.0)*(D-C)/N;
      Fout(ii,jj) = F(x,y); 
    end
  end

  % Compute the required boundary conditions: 
  % 
  BDB = 4*cos((Y+1)*(pi/2)); 
  % BDA, BDC, and BDD are dummy variables (in this example): 
  BDA = 0; BDC = 0; BDD = 0; 

  % Call our Helmholtz solver:
  % 
  [ U_approx, opts ] = hwscrt( Fout, BDA, BDB, BDC, BDD, ...
			       'A', A, 'B', B, 'M', M, 'xbc', 'dn', ...
			       'C', C, 'D', D, 'N', N, 'ybc', 'periodic', ...
			       'elmbda', ELMBDA ); 
  
  % Compare solution with truth: 
  clear U_true; 
  for ii=1:M+1,
    x = A + (ii-1.0)*(B-A)/M;
    for jj=1:N+1,
      y = C + (jj-1.0)*(D-C)/N;
      U_true(ii,jj) = (x^2)*cos((y+1)*(pi/2));
    end
  end
  
  % Save the output to be plotted by Matlab:
  save -v4 'test_conv.mat' 'X' 'Y' 'U_true' 'U_approx'; 

  % Compute the L^1 norm of the error:
  % 
  lInfty_error(gi) = max( max( abs( U_true-U_approx ) ) );

  if( 0 ) 
  figure; 
  mesh( X, Y, abs(U_true-U_approx).' ); 
  xlabel( 'X' ); ylabel( 'Y' ); zlabel( 'U_true-U_approx' ); 
  end

end

figure(1); plot( fliplr(1./nCells), fliplr(lInfty_error) ); 
xlabel( 'numerical cell size' ); ylabel( 'solution L^{infty} error' ); 
grid on; 
figure(2); loglog( fliplr(1./nCells), fliplr(lInfty_error) ); 
xlabel( 'numerical cell size' ); ylabel( 'solution L^{infty} error' ); 
grid on;