% 
% A test and verification of the routine hstcrt.f, using a 
% known analytical solution to the Cartesian Helmholtz equation.
% 
% Written by:
% -- 
% John L. Weatherwax                2005-10-29
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; 

A = 0; B = +2; M = 40; 
%MBDCND = 2; 
C = -1; D = +3; N = 80; 
%NBDCND = 0; 

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

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
figure; 
mesh( X, Y, Fout.' ); 
save -v4 'th_rhs.mat' 'X' 'Y' 'Fout'; 

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

[ 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 ); 

%U_approx = hwscrtg( A, B, M, MBDCND, BDA, BDB, C, D, N, NBDCND, BDC, BDD, ELMBDA, Fout ); 

% Plot solution: 
figure; mesh( X, Y, U_approx.' );
xlabel( 'X' ); ylabel( 'Y' ); zlabel( 'U_approx' ); 
save -v4 'th_uapprox.mat' 'X' 'Y' 'U_approx'; 

% Compare solution with truth: 
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
figure; 
mesh( X, Y, U_true.' ); 
xlabel( 'X' ); ylabel( 'Y' ); zlabel( 'U_true' ); 
save -v4 'th_true.mat' 'X' 'Y' 'U_true'; 

max( max( abs( U_true-U_approx ) ) ) 

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