/***********************  xyShadingPs   *********************
**  Author: Paul Soderlind http://www.hhs.se/personal/psoderlind/
**  Date:   19 May 2000
**  Caveats: no performance warrantee or guarantee, explict or implied
**
**  Purpose:
**   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);
**
**
**
**
------------------------------------------------------------------------------*/

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;
/*----------------------------------------------------------------------------*/
