pade_firstder(F,pade,Isym)

%% Presentation
%
%  Pade : 6th order finite difference derivative 
%  
%  .m file in the package:
%  pade_init.m      : initialization of the Pade coefficients
%  pade_firstder    : first order derivative
%  pade_secder      : second order derivative
%
%% 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
%  otherwise it leads to oscillations of the derivatives
%
%  2) The grid (xp) is generated inside the "pade_init" function.
%  A slight modification would authorize external grids
%  but beware to the boundary conditions 
%

function [DF] = pade_firstder(F,pade,Isym)
%% Head
%
% F                 : function to be derivated
% pade              : structure initilized in pade_init.m
% Isym              : symmetrical (=1) or anti-symmetrical (~=1) 
%                     function (if necessary)
%
N          = pade.nx;
DN         = 1/pade.dx;
NxM1       = N-1;
NxM2       = N-2;
NxM3       = N-3;
NxM4       = N-4;
ax(1:N)    = pade.ax(1:N,Isym,1);
bx(1:N)    = pade.bx(1:N,Isym,1);
cx(1:N)    = pade.cx(1:N,Isym,1);
A1         = 3;
B1         = (4*A1 + 2)*DN/3*0.5;
C1         = (4- A1)*DN/3*0.25;
A2         = 16*(2*A1 + 1)/(40- A1);
B2         = (4*A2 + 2)*DN/3*0.5;
C2         = (4- A2)*DN/3*0.25;
B3         = DN*3;
coef       = ax(1)+ax(NxM2)+bx(1);
if(Isym ~= 1) 
    Isym = 2;
end

%% Derivative process : left boundary
DF(3:NxM2)   = (F(4:NxM1)-F(2:NxM3))*B1 + (F(5:N)-F(1:NxM4))*C1;
if ( pade.bd0 == 0 ) 
    DF(1)        = (F(2)-F(NxM1))*B1    + (F(3)-F(NxM2))*C1;
    DF(2)        = (F(3)-F(1))*B1       + (F(4)-F(NxM1))*C1;
    DF(3)        = (F(4)-F(2))*B1       + (F(5)-F(1))   *C1;
elseif ( pade.bd0 == 1 )
    DF(1)        =((F(2)-F(1))*4        + (F(3)-F(1))  )*DN;
    DF(2)        = (F(3)-F(1))*B3;
    DF(3)        = (F(4)-F(2))*B2       + (F(5)-F(1))   *C2;
elseif( pade.bd0 == 2 )
    if( Isym == 2 ) 
        DF(1)    = (F(2)+F(1))*B1       + (F(3)+F(2))  *C1;
        DF(2)    = (F(3)-F(1))*B1       + (F(4)+F(1))  *C1;
    else 
        DF(1)    = (F(2)-F(1))*B1       + (F(3)-F(2))  *C1;
        DF(2)    = (F(3)-F(1))*B1       + (F(4)-F(1))  *C1;
    end
elseif( pade.bd0 == 3 ) 
    if( Isym == 2 ) 
        DF(1)    = (2*F(2))*B1          + (2*F(3))     *C1;
        DF(2)    = (F(3)-F(1))*B1       + (F(4)+F(2))  *C1;
    else 
        DF(1)    = 0;
        DF(2)    = (F(3)-F(1))*B1       + (F(4)-F(2))  *C1;
    end
else
    fprintf(1, 'error FirstDerivative \n');
    pause;
end

%% Derivative process : right boundary
if ( pade.bd1 == 0 ) 
    DF(NxM1)     = (F(1)-F(NxM2))*B1    + (F(2)-F(NxM3))*C1;
    DF(NxM2)     = (F(NxM1)-F(NxM3))*B1 + (F(1)-F(NxM4))*C1;
elseif ( pade.bd1 == 1 );
    DF(N)        =((F(N)-F(NxM1))*4     - (F(NxM2)-F(N)))*DN;
    DF(NxM1)     = (F(N)-F(NxM2))*B3;
    DF(NxM2)     = (F(NxM1)-F(NxM3))*B2 + (F(N)-F(NxM4)) *C2;
elseif( pade.bd1 == 2 )
    if( Isym == 2 ) 
        DF(N)    = (-F(N)-F(NxM1))*B1   + (-F(NxM1)-F(NxM2))*C1;
        DF(NxM1) = ( F(N)-F(NxM2))*B1   + (-F(N)-F(NxM3))   *C1;
    else 
        DF(N)    =( F(N)-F(NxM1))*B1    + ( F(NxM1)-F(NxM2))*C1;
        DF(NxM1) =( F(N)-F(NxM2))*B1    + ( F(N)-F(NxM3))   *C1;
    end
elseif( pade.bd1 == 3 );
    if( Isym == 2 ) 
        DF(N)    = -2*F(NxM1)*B1        - 2*F(NxM2)*C1;
        DF(NxM1) = (F(N)-F(NxM2))*B1    + (-F(NxM1)-F(NxM3))*C1;
    else 
        DF(N)    = 0;
        DF(NxM1) = (F(N)-F(NxM2))*B1    + (F(NxM1)-F(NxM3))*C1;
    end
else
        fprintf(1, 'error FirstDerivative \n');
        pause;
end

%% Thomas Algorithm phase II
if ( pade.bd0 == 0 )
    % periodic 
    DF(1)    =  DF(1)*cx(1);
    for I = 2:NxM2
    DF(I)   = (DF(I)-DF(I-1))*cx(I);
    end
    for I = NxM3:-1:1
        DF(I)   =  DF(I)-DF(I+1)*cx(I); 
    end
    %
    DF(NxM1) = (DF(NxM1)-DF(1)-DF(NxM2))/coef;
  %
    for I = 1:NxM2
        DF(I)   = DF(I)+DF(NxM1)*ax(I);
    end
  %
    DF(N)    =  DF(1);
else
  % non periodic 
    DF(1)    =  DF(1)/bx(1);
    for I = 2:NxM1
        DF(I)   = (DF(I)-DF(I-1))/bx(I);
    end
  %
    DF(N)    = (DF(N)-DF(NxM1)*ax(N))/bx(N);
  %
    for I = NxM1:-1:1
        DF(I)   = DF(I)-DF(I+1)*cx(I);
    end
end
%
end