pade_test(nperiod,b0,b1)

%% Presentation
%
%  Pade : 6th order finite difference derivative 
%  
%  .m file in the package:
%  pade_init.m      : initialization of the Pade coefficients
%  pade_firstder.m  : first order derivative
%  pade_secder.m    : second order derivative
%  pase_test.m      : this file as an exemple
%
%% Author
%  Julien Reveillon -- Julien.Reveillon(at)coria.fr
%  CORIA, University of Rouen, France
%
%% History
%  January 27, 2007 : Version 1.0 = 1D, regular grid
%                       
%% Remarks
%
%  1) A sixth order FD derivative can not be used with too strong gradients
%  
%  2) The grid (xp) is generated inside the pade_init function.
%  A slight modification would authorize external grids
%  but beware to the boundary conditions
%% How to Use it ?
%  runpade(3,0,0) : to test periodic derivatives of cos(x) (3 wave lengths)
%  runpade(4,1,1) : with left and right boundaries
%  runpade(2,2,1) : with left symetric and right boundaries
%   etc
%  all combinations are possible expect of course with the periodic BC that
%  must be 0 left and right.
function pade_test(nperiod,b0,b1)
%
% pase_test is a test function of the pade package.
%
% nperiod   : number of periods in the periodical test functions
% b0        : left boundary condition
% b1        : right boundary condition
%     
%--------------------------------
% User Data
%--------------------------------
% pade.nx   : Number of nodes (isotropic grid)
pade.nx = 400;
% pade.l0   : starting coordinate of the domain
pade.l0 = 0;
% pade.lx   : Length of the domain
pade.lx = 1;
%  pade.bdO and .bd1 : Boundary condition
% 0 : periodic
% 1 : bounded
% 2 : non centered symmetric
% 3 : centered symmetric
pade.bd0= b0   % left boundary condition
pade.bd1= b1   % right boundary condition (periodic : 0 left and right)
%--------------------------------
% Pade Private Data : just call it once
%--------------------------------
[pade.ax pade.bx pade.cx pade.xp pade.dx] = ...
    pade_init(pade.nx,pade.l0,pade.lx,pade.bd0,pade.bd1);
%
% Function initialization to test the subroutine
%
for i= 1:pade.nx
 x      = pade.xp(i);
 % function to derive
 ASF(i)   = pade.l0+sin(nperiod*x*2*pi/(pade.lx-pade.l0)); % anti-symmetrical (AS)
 SF(i)    = pade.l0+cos(nperiod*x*2*pi/(pade.lx-pade.l0)); % symmetrical (S)
 % Exact solution
 DASFE(i) =  nperiod*2*pi*cos(nperiod*x*2*pi/(pade.lx-pade.l0))/(pade.lx-pade.l0);        %D1, AS
 D2ASFE(i)= -nperiod^2*4*pi^2*sin(nperiod*x*2*pi/(pade.lx-pade.l0))/(pade.lx-pade.l0)^2; %D2, AS
 DSFE(i)  = -nperiod*2*pi*sin(nperiod*x*2*pi/(pade.lx-pade.l0))/(pade.lx-pade.l0);       %D1, S
 D2SFE(i) = -nperiod^2*4*pi^2*cos(nperiod*x*2*pi/(pade.lx-pade.l0))/(pade.lx-pade.l0)^2; %D2, S
end
%
%
% Computation of the first derivative
[DASF] = pade_firstder(ASF,pade,2); % Antisymmetrical function
[DSF]  = pade_firstder(SF,pade,1);  % Symmetrical function
%
figure(1);
plot(pade.xp,DASFE,'ok',pade.xp,DSFE,'ob');
hold on;
plot(pade.xp,DASF,'-g',pade.xp,DSF,'-r');
hold off;
%
%
% computation of the second derivative
[D2ASF] = pade_secder(ASF,pade,2);  % Antisymmetrical function
[D2SF]  = pade_secder(SF,pade,1);    % Symmetrical function
%
figure(2);
plot(pade.xp,D2ASFE,'ok',pade.xp,D2SFE,'ob');
hold on;
plot(pade.xp,D2ASF,'-g',pade.xp,D2SF,'-r');
hold off;
%
%figure(3);
%plot(pade.xp,ASF,'ok',pade.xp,SF,'ob');
%hold off
%
% maximal error : Anti-symmetrical
%
diff1max = -1e30;
diff2max = -1e30;
damp     = nperiod*2*pi/(pade.lx-pade.l0);       %max amplitude
d2amp    = nperiod^2*4*pi^2/(pade.lx-pade.l0)^2; %max amplitude
for i= 1:pade.nx
    if(abs(DASF(i)-DASFE(i)) > diff1max)
        diff1max = abs(DASF(i)-DASFE(i));
        i1max    = i;
    end
    if(abs(D2ASF(i)-D2ASFE(i)) > diff2max)
        diff2max = abs(D2ASF(i)-D2ASFE(i));
        i2max = i;
    end
end
fprintf(1,'Anti-symetric test.  percent error max f '' : (position: %i, error d1: %d)\n',i1max,diff1max/damp*100);
fprintf(1,'Anti-symetric test.  percent error max f'''': (position: %i, error d2: %d)\n',i2max,diff2max/d2amp*100);
%
% maximal error : Symmetrical
%
diff1max = -1e30;
diff2max = -1e30;
for i= 1:pade.nx
    if(abs(DSF(i)-DSFE(i)) > diff1max)
        diff1max = abs(DSF(i)-DSFE(i));
        i1max    = i;
    end
    if(abs(D2SF(i)-D2SFE(i)) > diff2max)
        diff2max = abs(D2SF(i)-D2SFE(i));
        i2max = i;
    end
end
fprintf(1,'Symetric test. percent error max f'' : (position: %i, error d1: %d)\n',i1max,diff1max/damp*100);
fprintf(1,'Symetric test. percent error max f'''': (position: %i, error d2: %d)\n',i2max,diff2max/d2amp*100);
%
end