/*------------------------------------------------------------------------------
**  psspline.src
**  author: Paul Soderlind mailto:Paul.Soderlind@hhs.se
**  date:   18 May 2000
**  Code provided gratis without any performance warrantee.
**
**
**
**
**
**
**
**
**
**
------------------------------------------------------------------------------*/



/*------------------------------------------------------------------------------
**  CubicSplinePs
**
**  Purpose:   Cubic spline interpolation of some function f(x)
**
**
**
**
**  Usage:     y2M = CubicSplinePs(y,x,D2y_1,D2y_N,x2,LinearIt);
**
**  Inputs:    y          Nx1 vector of function values, y=f(x)
**             x          Nx1 vector
**             D2y_1      scalar, f''(x[1]), often 0
**             D2y_N      scalar, f''(x[N]), often 0
**             x2         MxK matrix, points where f() should be approximated
**             LinearIt   scalar, 1: make linear interpolation instead
**
**  Output:    y2         MxK matrix, approximately f(x2)
**
**
**
**  Note 1:    This proc calls on SplinePs for calculating the
**             second derivatives, then on SplintPs to do the
**             interpolation
**
**
**
**
**
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), 18 May 2000
------------------------------------------------------------------------------*/

proc CubicSplinePs(y,x,D2y_1,D2y_N,x2,LinearIt);

  local D2y,y2;

  if LinearIt == 1;               /*only linear interpolation, impose f''(x)=0*/
    D2y = zeros( rows(y),1 );
  else;
    D2y = SplinePs(y,x,D2y_1,D2y_N);        /*calculate second order derivatives*/
  endif;

  y2  = SplintPs(y,x,D2y,x2);               /*make interpolation*/

  retp(y2);

endp;
/*----------------------------------------------------------------------------*/




/*------------------------------------------------------------------------------
**  SplintPs
**
**  Purpose:   Cubic spline interpolation of some function f(x), using already
**             available second derivatives, f''(x)
**
**
**
**
**  Usage:     y2 = SplintPs(y,x,D2y,x2);
**
**  Inputs:    y          Nx1 vector of function values, y=f(x)
**             x          Nx1 vector
**             D2y        Nx1, f''(x])
**             x2         MxK matrix, points where f() should be approximated
**
**  Output:    y2         MxK matrix, approximately f(x2)
**
**
**
**
**  Note 2:    This code uses the notation in Press et al,
**             Numerical Recipes in Fortran, 2nd ed, 107-109,
**             but is not a translation of the routines there
**
**
**
**
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), 18 May 2000
------------------------------------------------------------------------------*/

proc SplintPs(y,x,D2y,x2);


  local M,K,y2,i,kk,x2_i,j,A,B,C,D;

  M = rows(x2);
  K = cols(x2);


  y2 = miss( zeros(M,K),0 );    /*to put the results in: approximate y*/

  i = 1;
  do until i > M;                      /*loop over rows in x2*/
    kk = 1;
    do until kk > K;                   /*loop over columns in x2*/
      x2_i = x2[i,k];
      j = (x.<=x2_i) .and (x2_i.<lagn(x,-1));       /*x(j) <= x2_i < x(j+1)*/
      j = indexcat(j,1);                       /*find that j*/
      if scalmiss(j) == 0;                     /*if x2_i is within the range*/
        A = (x[j+1]-x2_i)/(x[j+1]-x[j]);
        B = 1 - A;
        C = (1/6)*(A^3-A)*(x[j+1]-x[j])^2;
        D = (1/6)*(B^3-B)*(x[j+1]-x[j])^2;
        y2[i,k] = A*y[j] + B*y[j+1] + C*D2y[j] + D*D2y[j+1];
      else;
        y2[i,k] = error(0);
      endif;
      kk = kk + 1;
    endo;                                   /*end loop over colums in x2*/
    i = i + 1;
  endo;                                     /*end loop over rows in x2*/


  retp(y2);

endp;
/*----------------------------------------------------------------------------*/



/*------------------------------------------------------------------------------
**  SplinePs
**
**  Usage:     calculates the 2nd derivatives, f''(x) for the cubic spline
**
**  Usage:   D2y = SplinePs(y,x,D2y_1,D2y_N);
**
**  Inputs:    y          Nx1 vector of function values, y=f(x)
**             x          Nx1 vector
**             D2y_1      scalar, f''(x[1]), often 0
**             D2y_N      scalar, f''(x[N]), often 0
**
**  Output:    Dy2        Nx1 vector, f''(x)
**
**
**
**  Note:      This code uses the notation in Press et al,
**             Numerical Recipes in Fortran, 2nd ed, 107-109,
**             but is not a translation of the routines there
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), 18 May 2000
------------------------------------------------------------------------------*/

proc SplinePs(y,x,D2y_1,D2y_N);

  local N,m,r,a,b,c,D2y;


  N = rows(y);
  m = N - 2;                    /*dimension of system*/

  r = (lagn(y,-1)-y)./(lagn(x,-1)-x) - (y - lagn(y,1))./(x-lagn(x,1));
  r[1] = r[1] - (x[2]-x[1])*D2y_1/6;
  r[m] = r[m] - (x[N]-x[N-1])*D2y_N/6;

  a = (x - lagn(x,1)          )/6;  /*vectors in tridiagonal matrix*/
  b = (lagn(x,-1) - lagn(x,1) )/3;
  c = (lagn(x,-1) - x         )/6;

  r = trimr(r,1,1);              /*j=2,..,N-1*/
  a = trimr(a,1,1);
  b = trimr(b,1,1);
  c = trimr(c,1,1);

  D2y = TriDiagSolvePs( a,b,c,r );

  D2y =  D2y_1|D2y|D2y_N;


  retp(D2y);


endp;
/*----------------------------------------------------------------------------*/



/*-----------------------------------------------------------------------
**  TriDiagSolvePs
**
**  Purpose:   Solving ABC*u = r for u, where
**
**             ABC   =    b1 c1 0     ...     ...         0
**                        a2 b2 c2 0  ...     ...         0
**                        0  a3 b3 c2 ...     ...         0
**                        . .      ...                    .
**                        .   .           a(N-1)  b(N-1)  c(N-1)
**                        .     .  ...    0       aN      bN
**
**            using translation of FORTRAN code (to Gauss code) from
**            Numerical Recipes, 2nd ed
**
**
**  Usage:     u = TriDiagSolvePs( a,b,c,r );
**
**  Input:     a is an Nx1 vector
**             b    ""
**             c    ""
**             r    ""
**
**  Output:    u is an Nx1 vector with the solution
**
**
**  Reference: Numerical Recepies, 2nd edition, p 43
**
**
**
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), May 2000
-----------------------------------------------------------------------*/

proc TriDiagSolvePs( a,b,c,r );

  local N,bet,gam,u,j;

  N   = rows(b);
  u   = miss( zeros(N,1),0 );
  gam = u;


  if b[1] == 0;
    "TriDiagSolvePs: rewrite equations";
    retp( error(0) );
  endif;


  bet = b[1];
  u[1] =r[1]/bet;
  j = 2;
  do until j > N;                         /*forward substitution*/
    gam[j] = c[j-1]/bet;
    bet    = b[j] - a[j]*gam[j];
    u[j]   = (r[j]-a[j]*u[j-1])/bet;
    j = j + 1;
  endo;

  j = n - 1;                              /*backsubstitution*/
  do until j < 1;
    u[j] = u[j] - gam[j+1]*u[j+1];
    j = j - 1;
  endo;


retp(u);

endp;
/*----------------------------------------------------------------------------*/



/*------------------------------------------------------------------------------
**  TriDiagToDensePs
**
**  Purpose:   Creating the denbse ABC matrix
**
**             ABC   =    b1 c1 0     ...     ...         0
**                        a2 b2 c2 0  ...     ...         0
**                        0  a3 b3 c2 ...     ...         0
**                        . .      ...                    .
**                        .   .           a(N-1)  b(N-1)  c(N-1)
**                        .     .  ...    0       aN      bN
**
**            from the vectors a,b, and c
**
**
**  Usage:     abc = TriDiagToDensePs( a,b,c );
**
**  Input:     a is an Nx1 vector
**             b    ""
**             c    ""
**
**
**  Output:    abc is an NxN matrix
**
**
**
**
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), 18 May 2000
------------------------------------------------------------------------------*/

proc TriDiagToDensePs( a,b,c );

  local N, abc,i;

  N   = rows(b);
  abc = zeros(N,N);

  i = 1;
  do until i > N;
    if i > 1;
      abc[i,i-1] = a[i];
    endif;
    abc[i,i] = b[i];
    if i < N;
      abc[i,i+1] = c[i];
    endif;
    i = i + 1;
  endo;


  retp(abc);

endp;
/*----------------------------------------------------------------------------*/



/*------------------------------------------------------------------------------
**  xyShadingPs
**
**  Purpose:   to create shading in xy graphs by creating a _pline matrix
**
**
**
**
**
**  Usage:     pline = xyShadingPs( yLo,yHi,x,xmin,xmax,Nline,theta,
**                                  StepSiz,LinearIt,
**                                  lintype,lincolor,linthick );
**
**  Inputs:    yLo              Nx1 vector, function yLo(x), lower boundary
**             yHi              Nx1 vector, function yHi(x), upper boundary
**             x                Nx1 vector, where yLo and yHi are evaluated
**             xmin             scalar, lower x boundary
**             xmax             scalar, upper x boundary
**             Nline            scalar, number of shading lines
**             theta            scalar, angle of lines, in radians
**             StepSiz          scalar, step size in trying sloping lines
**             LinearIt         scalar, 1: only linear interpolation of YLo and
**                              yHi, otherwise cubic spline
**             LinType          scalar, 6: solid; 4: closely spaced dots
**             LinColor         scalar, line colour, 8: dark grey
**             LinThick         scalar, line thickness, 0: default
**
**  Output:    pline            Nline x 9 matrix, to be used as _pline
**
**
**
**
**
**
**  Example:   x = seqa(0,(2+pi)/99,100);
**             y = cos(x) ~ sin(x);
**             graphset;
**             _pline = xyShadingPs( y[.,1],y[.,2],x,1,3,15,pi/2.5,0.001,
**                                   0,6,8,0.1 );
**             _pcolor = { 8 };                  @dark grey@
**             _pmcolor = { 8, 8, 8, 8, 8, 8, 8, 8, 15 };
**             xy(x,y);
**
**
**
**
**  Paul Soderlind (Paul.Soderlind@hhs.se), 19 May 2000
------------------------------------------------------------------------------*/

proc xyShadingPs( yLo,yHi,x,xmin,xmax,Nline,theta,StepSiz,LinearIt,
                  lintype,lincolor,linthick );


  local pline,x0,x1,y0,y1,D2y,b,i,StopIt,xx,yy,yCub;


  pline = 1 ~lintype ~0 ~0 ~1 ~1 ~1 ~lincolor ~linthick;
  pline = pline.*ones(Nline,1);


  x0 = seqa(xmin,(xmax-xmin)/(Nline-1),Nline);  /*lower left x points*/
    x1 = x0;                                      /*upper right x points*/
  if rows(yLo) == 1;                            /*lower left y points*/
    y0 = yLo.*ones(Nline,1);
  else;
    y0 = CubicSplinePs(yLo,x,0,0,x0,LinearIt);  /*cubic or linear spline*/
  endif;
  if rows(yHi) == 1;                            /*upper right y points*/
    y1 = yHi.*ones(Nline,1);
  else;
    y1 = CubicSplinePs(yHi,x,0,0,x1,LinearIt);
  endif;


  if (theta>=0) and  (theta.<pi/2);     /*SLOPING LINES*/
    if LinearIt == 1;               /*linear interpolation, impose f''(x)=0*/
      D2y = zeros( rows(yHi),1 );
    else;
      D2y = SplinePs(yHi,x,0,0);    /*calculate second order derivatives*/
    endif;
    b = tan(theta);                    /*slope of line*/
    i = 1;
    do until i > Nline;
      StopIt = 0;
      xx = x0[i];
      do until StopIt == 1;
        xx = xx + StepSiz;
        yy = y0[i] + b*(xx-x0[i]);
        yCub = SplintPs(yHi,x,D2y,xx);               /*make interpolation*/
        StopIt = (xx.>xmax) .or (yy.>yCub);   /*is xy outside, is yy outside*/
      endo;
      x1[i] = xx;
      y1[i] = yy;
      i = i + 1;
    endo;

  endif;                           /*END SLOPING LINES*/

  pline[.,3:6] = x0 ~ y0 ~ x1 ~ y1;   /*put lines in correct places in pline*/

  retp(pline);

endp;
/*----------------------------------------------------------------------------*/
