%{
#Function  dif1D

%}

function [d,dd,w,s]=dif1D(type,s0,L,N,opt);

%{ 

This function creates a mesh s, differentiation matrices of first and second
order d and dd, as well as the weight w. 
 
 On output: s is a column-vector containing the mesh points, starting from the point at the left.
            d and dd are the derivation matrices (quare matrices, in full or sparse storing format)
            w is a line-vector containing the weights. 
            
Last parameter 'opt' is optional. For Finite-Difference it is interpreted as
the number of points of the stencil (default :3) and for
Stretched Chebyschev it is interpreted as the parameter b
controlling the stretching.

%}

switch type

%{
## Finite differences (optimized case)

For a domain [s0,s0+L]

Last parameter 'pts' is optional, if absent it is set to 3 (relevant value in almost all cases)
%}
  
   case 'fd' 
       
  if(nargin==4) 
    pts = 3;
  else
   pts = opt;
  end
  scale=L/2;
  [s,d] = fddif(N,1,pts); 
  [s,dd] = fddif(N,2,pts); 
  s=s*scale; 
  d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
  w=([diff(s'),0]+[0,diff(s')])/2;    
  
  
  %{ 
  ## Finite differences (basic implementation, equivalent to the previous)
  %}   
    
 case 'fds' 
    

  dx = L/(N-1); % grid spacing
  s = [  s0 : L/(N-1) : s0+L ]';  % creates an equispaced mesh with N point on the interval
  for i=2:N-1 % centered three-point formulas 
    d(i,i) = 0;
    d(i,i+1) =  1/(2*dx); 
    d(i,i-1) = -1/(2*dx);
    dd(i,i) = -2/dx^2;
    dd(i,i+1) =  1/(dx^2); 
    dd(i,i-1) =  1/(dx^2);        
  end
 % uncentred formulas for first and last grid points 
   d(1,1) = -3/(2*dx); d(1,2) = 4/(2*dx) ; d(1,3) = -1/(2*dx);
   d(N,N-2) =  1/(2*dx); d(N,N-1) = -4/(2*dx); d(N,N)   =  3/(2*dx);
   dd(1,1) = 2/dx^2; dd(1,2) = -5/dx^2; dd(1,3) = 4/dx^2; dd(1,4) = -1/dx^2; 
   dd(N,N) = 2/dx^2; dd(N,N-1) = -5/dx^2; dd(N,N-2) = 4/dx^2; dd(N,N-3) = -1/dx^2; 
  % "weight" to compute integrals using the trapeze rule
  for i=2:N-1
    w(i) = dx;
  end
  w(1) = dx/2; w(N) = dx/2;
  
%{ 
## Finite differences on a periodic domain 
 
 Domain is defined as [s0,s0+L]
%}
  
 case 'fp' % periodic finite differences 
  scale=L/2;
  [s,d] = fdper(N,1,opt); 
  [s,dd] = fdper(N,2,opt);
  s=s*scale; 
  d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
  w=ones(1,N)*(s(2)-s(1)); 

%{ 
## Chebyshev collocation 

on domain [s0,s0+L]
%}  

 case 'cheb' % chebychev 
  scale=-L/2;
  [s,DM] = chebdif(N,2);
  d=DM(:,:,1);  
  dd=DM(:,:,2);
  s=s*scale; 
  d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
  w=L*clencurt(N)/2; 
  
  
%{ 
## Fourier collocation 
   
   (requires a periodic domain [s0,s0+L])
  
%} 
  case 'fou' % Fourier 
  scale=L/(2*pi); 
  [s, d] = fourdif(N, 1);
  [s, dd] = fourdif(N, 2);
  s=s*scale; 
  w=ones(1,N)*(s(2)-s(1));
  d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;

%{ 
## Hermite collocation on infinite domain.
  
  Domain is centred on s0, L controls the density 
  
  (increase L to cluster the points towards s0) 

%} 

  case 'her' % hermite 
  scale=1;
  [s,DM] = herdif(N,2,L);
  d=DM(:,:,1);  
  dd=DM(:,:,2);
  s=s*scale; 
  d=d/scale; dd=dd/scale^2; s=s+s0;
  w=zeros(1,N); % weight is not implemented yet !!!
  
%{
## Stretched Chebyshev ; infinite ; exponential
  
  * s0 is the center of the domain.
  * L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
  * optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
%}
  
  case 'chebInfExp' % chebychev with coordinate stretching to treat large domains
      
 
   if(nargin==4) 
      b = 0.999;
   else
      b = opt;
   end
   
   [dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
   
  s = L*atanh(xi)+s0; % coordinate stretching (exponential)
  G(:, 1) =                          - (xi.^2 - 1) /L;      % G1 = ds / dxi
  G(:, 2) =                   2 *xi .* (xi.^2 - 1) /L^2;    % G2 = d^2 s /d xi^2
  G(:, 3) =    - 2 * (3 * xi.^2 - 1) .* (xi.^2 - 1) / L^3;
  G(:, 4) = 8 * xi .* (3 * xi.^2 - 2) .* (xi.^2 - 1) / L^4;
  
  d = G(:,1).*dxi;
  dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
  w = wxi./G(:,1)';
  
%{
## Stretched Chebyshev ; infinite ; algebraic
  
  * s0 is the center of the domain.
  * L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
  * optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
  %}
  
  case 'chebInfAlg' % chebychev with coordinate stretching to treat large domains
      
  if(nargin==4) 
      b = 0.999;
   else
      b = opt;
  end
  
  [dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
  
  s = L * xi ./ sqrt(1 - xi.^2); 
  K = sqrt(s.^2 + L^2);
  G(:, 1) =                                    L^2 ./ K.^3;
  G(:, 2) =                          - 3 * s * L^2 ./ K.^5;
  G(:, 3) =             3 * (4 * s.^2 - L^2) * L^2 ./ K.^7;
  G(:, 4) = - 15 * s .* (4 * s.^2 - 3 * L^2) * L^2 ./ K.^9;
 
  d = G(:,1).*dxi;
  dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
  w = wxi./G(:,1)';
  
 %{
## Stretched Chebyshev ; semi-infinite ; exponential
  
  * s0 is the left bound of the domain.
  * L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
  * optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
  %}
  
 case 'chebSemiInfExp' % chebychev with coordinate stretching to treat large domains
 
 if(nargin==4) 
      b = 0.999;
   else
      b = opt;
 end
     
  [dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,b]
  
  s = L*atanh(xi)+s0; % coordinate stretching (exponential)
  G(:, 1) =                          - (xi.^2 - 1) /L;      % G1 = ds / dxi
  G(:, 2) =                   2 * xi .* (xi.^2 - 1) /L^2;    % G2 = d^2 s /d xi^2
  d = G(:,1).*dxi;
  dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
  w = wxi./G(:,1)';

  %{
## Stretched Chebyshev ; semi-infinite ; algebraic
  
  * s0 is the left bound of the domain.
  * L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
  * optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
  %}
  
  case 'chebSemiInfAlg' % chebychev with coordinate stretching to treat large domains
      
  if(nargin==4)  
      b = 0.999;  
  else
      b = opt; 
  end
   
  [dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,1]

  s = L * xi ./ sqrt(1 - xi.^2); 
  K = sqrt(s.^2 + L^2);
  G(:, 1) =                                    L^2 ./ K.^3;
  G(:, 2) =                          - 3 * s .* L^2 ./ K.^5;
 
  d = G(:,1).*dxi;
  dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
  w = wxi./G(:,1)';
  

end