/*
 * MLINTP.C - interpolation routines (N-dimensional and otherwise)
 *
 * Source Version: 2.0
 * Software Release #92-0043
 *
 */

#include "cpyright.h"

#include "pml.h"

static REAL
 *SC_DECLARE(_PM_redist_nodes_logical, (REAL *f,
				     int km, int lm, int kmax, int lmax,
				     char *emap));

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_CONNECTIVITY - return a pointer to the mesh topology for mappings on
 *                 - Arbitrarily-Connected sets or
 *                 - to the maximum indices for Logical-Rectangular sets
 */

byte *PM_connectivity(f)
   PM_mapping *f;
   {byte *cnnct;

    if (strcmp(f->category, PM_LR_S) == 0)
       cnnct = (byte *) f->domain->max_index;

    else if (strcmp(f->category, PM_AC_S) == 0)
       cnnct = (byte *) f->domain->topology;

    else
       cnnct = NULL;

    return(cnnct);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_ZONE_NODE_LR_2D - given a Logical-Rectangular zone centered mesh array
 *                    - return a node centered version
 *                    - by redistributing the values to the node
 *                    - using a uniform fractional value (1/4)
 */

REAL *PM_zone_node_lr_2d(f, cnnct, alist)
   REAL *f;
   byte *cnnct;
   pcons *alist;
   {int *maxes, n, km, lm, kmax, lmax, eflag;
    char *emap;
    REAL *ret;

    maxes = (int *) cnnct;
    kmax  = maxes[0];
    lmax  = maxes[1];

    emap  = NULL;
    SC_assoc_info(alist,
		  "EXISTENCE", &emap,
		  NULL);

    n = SC_arrlen(f)/sizeof(REAL);
    if (n == (kmax - 1)*(lmax - 1))
       {km = kmax - 1;
	lm = lmax - 1;}

    else
       return(NULL);

    eflag = FALSE;
    if (emap != NULL)
       PM_CHECK_EMAP(alist, n, eflag, emap);

    ret = _PM_redist_nodes_logical(f, km, lm, kmax, lmax, emap);

    if (eflag)
       SFREE(emap);

    return(ret);}
                       
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_NODE_ZONE_LR_2D - given a Logical-Rectangular node centered mesh
 *                    - mesh array return a zone centered version
 *                    - by redistributing the values to the zone
 *                    - using a uniform fractional value (1/4)
 */

REAL *PM_node_zone_lr_2d(f, cnnct, alist)
   REAL *f;
   byte *cnnct;
   pcons *alist;
   {int *maxes, kmax, lmax, *pc, corner;
    int i, j, km, lm, n, eflag, npts, delta;
    char *emap;
    REAL *fp, *f1, *f2, *f3, *f4;

    maxes = (int *) cnnct;
    kmax  = maxes[0];
    lmax  = maxes[1];

    emap  = NULL;
    SC_assoc_info(alist,
		  "EXISTENCE", &emap,
		  "CORNER", &pc,
		  NULL);

    corner = (pc == NULL) ? 2 : *pc;

    npts  = kmax*lmax;
    eflag = (emap == NULL);
    if (eflag)
       {emap = FMAKE_N(char, npts, "PM_NODE_ZONE_LR_2D:emap");
	memset(emap, 1, npts);}
    else
       PM_CHECK_EMAP(alist, npts, eflag, emap);

    switch (corner)
       {case 1 :
	     delta = 1;
	     break;

	default :
        case 2 :
	     delta = kmax + 1;
	     break;

        case 3 :
	     delta = kmax;
	     break;

        case 4 :
	     delta = 0;
	     break;};

    emap += delta;

    fp = FMAKE_N(REAL, npts, "PM_NODE_ZONE_LR_2D:fp");
    if ((_SC_zero_space != 1) && (_SC_zero_space != 2))
       PM_set_value(fp, npts, 0.0);
    PM_LOGICAL_ZONE(f, f1, f2, f3, f4, kmax);

    km = kmax - 1;
    lm = lmax - 1;
    n  = km*lm;
    for (j = 0; j < n; j++)
        {i = j + j/km;
	 if (emap[i] != 0)
            fp[i] = 0.25*(f1[i] + f2[i] + f3[i] + f4[i]);};

    emap -= delta;

    if (eflag)
       SFREE(emap);

    return(fp);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_ZONE_NODE_AC_2D - given a Arbitrarily-Connected zone centered mesh
 *                    - mesh array return a node centered version
 *                    - by redistributing the values to the node
 *                    - using a uniform fractional value
 */

REAL *PM_zone_node_ac_2d(f, cnnct, alist)
   REAL *f;
   byte *cnnct;
   pcons *alist;
   {PM_mesh_topology *mt;
    REAL *fp, fv;
    long **cells, *zones, *sides;
    int *nc, nz, *np, nzp, nsp, nn;
    int in, iz, is, is1, is2, os, oz;

    mt = (PM_mesh_topology *) cnnct;

    cells = mt->boundaries;
    zones = cells[2];
    sides = cells[1];

    nc = mt->n_cells;
    nz = nc[2];
    nn = nc[0];

    np  = mt->n_bound_params;
    nzp = np[2];
    nsp = np[1];

    fp = FMAKE_N(REAL, nn, "PM_ZONE_NODE_AC_2D:fp");
    if ((_SC_zero_space != 1) && (_SC_zero_space != 2))
       PM_set_value(fp, nn, 0.0);
    np = FMAKE_N(int, nn, "PM_ZONE_NODE_AC_2D:np");
    if ((_SC_zero_space != 1) && (_SC_zero_space != 2))
       memset(np, 0, nn*sizeof(int));

/* accumulate nodal values from the zones */
    for (iz = 0; iz < nz; iz++)
        {oz  = iz*nzp;
	 is1 = zones[oz];
	 is2 = zones[oz+1];

         fv  = f[iz];

	 for (is = is1; is <= is2; is++)
	     {os = is*nsp;
	      in = sides[os];
              fp[in] += fv;
              np[in]++;};};
    
/* normalize the nodal values */
    for (in = 0; in < nn; in++)
        {fv = (REAL) np[in] + SMALL;
	 fp[in] /= fv;};

    SFREE(np);

    return(fp);}
                       
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_NODE_ZONE_AC_2D - given a Arbitrarily-Connected node centered mesh
 *                    - mesh array return a zone centered version
 *                    - by redistributing the values to the zone
 *                    - using a uniform fractional value
 */

REAL *PM_node_zone_ac_2d(f, cnnct, alist)
   REAL *f;
   byte *cnnct;
   pcons *alist;
   {PM_mesh_topology *mt;
    REAL *fp, fv;
    long **cells, *zones, *sides;
    int *nc, nz, *np, nzp, nsp;
    int in, iz, is, is1, is2, os, oz;

    mt = (PM_mesh_topology *) cnnct;

    cells = mt->boundaries;
    zones = cells[2];
    sides = cells[1];

    nc = mt->n_cells;
    nz = nc[2];

    np  = mt->n_bound_params;
    nzp = np[2];
    nsp = np[1];

    fp = FMAKE_N(REAL, nz, "PM_NODE_ZONE_AC_2D:fp");

/* accumulate nodal values from the zones */
    for (iz = 0; iz < nz; iz++)
        {oz  = iz*nzp;
	 is1 = zones[oz];
	 is2 = zones[oz+1];

         fv  = 0.0;
	 for (is = is1; is <= is2; is++)
	     {os  = is*nsp;
	      in  = sides[os];
              fv += f[in];};

         fp[iz] = fv/((REAL) (is2 - is1 + 1));};

    return(fp);}
                       
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* _PM_REDIST_NODES_LOGICAL - given a logical rectangular zone centered
 *                          - mesh array return a node centered version
 *                          - by redistributing the values to the node
 *                          - using a uniform fractional value (1/4)
 */

static REAL *_PM_redist_nodes_logical(f, km, lm, kmax, lmax, emap)
   REAL *f;
   int km, lm, kmax, lmax;
   char *emap;
   {int i, j, k, l, nn, nz, eflag;
    REAL val;
    REAL *fp, *fp1, *fp2, *fp3, *fp4;
    REAL *ip, *ip1, *ip2, *ip3, *ip4;

    nn    = kmax*lmax;
    eflag = (emap == NULL);
    if (eflag)
       {emap = FMAKE_N(char, nn, "_PM_REDIST_NODES_LOGICAL:emap");
	memset(emap, 1, nn);};

    ip = FMAKE_N(REAL, nn, "_PM_REDIST_NODES_LOGICAL:ip");
    fp = FMAKE_N(REAL, nn, "_PM_REDIST_NODES_LOGICAL:fp");
    if ((_SC_zero_space != 1) && (_SC_zero_space != 2))
       {PM_set_value(ip, nn, 0.0);
	PM_set_value(fp, nn, 0.0);};

    PM_LOGICAL_ZONE(fp, fp1, fp2, fp3, fp4, kmax);
    PM_LOGICAL_ZONE(ip, ip1, ip2, ip3, ip4, kmax);

    nz = km*lm;
    for (j = 0; j < nz; j++)
        {if (emap[j] != 0)
            {k = j % km;
	     l = j / km;
	     i = l*kmax + k;

	     val = f[j];

	     ip1[i]++;
	     ip2[i]++;
	     ip3[i]++;
	     ip4[i]++;

	     fp1[i] += val;
	     fp2[i] += val;
	     fp3[i] += val;
	     fp4[i] += val;};};

    for (j = 0; j < nn; j++)
        if (ip[j] > 0.0)
           fp[j] /= ip[j];
        else
           fp[j] = 0.0;

    SFREE(ip);
    if (eflag)
       SFREE(emap);

    return(fp);}
                       
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_ZONE_CENTERED_MESH_2D - compute a zone centered set of coordinates
 *                          - from the given logical rectangular mesh
 */

void PM_zone_centered_mesh_2d(px, py, rx, ry, kmax, lmax)
   REAL **px, **py, *rx, *ry;
   int kmax, lmax;
   {int i, n;
    REAL *xc, *x1, *x2, *x3, *x4;
    REAL *yc, *y1, *y2, *y3, *y4;

    n = kmax*lmax;

    xc = FMAKE_N(REAL, n, "PM_ZONE_CENTERED_MESH_2D:xc");
    yc = FMAKE_N(REAL, n, "PM_ZONE_CENTERED_MESH_2D:yc");

    PM_LOGICAL_ZONE(rx, x1, x2, x3, x4, kmax);
    PM_LOGICAL_ZONE(ry, y1, y2, y3, y4, kmax);

    for (i = 0; i < n; i++)
        {xc[i] = 0.25*(x1[i] + x2[i] + x3[i] + x4[i]);
         yc[i] = 0.25*(y1[i] + y2[i] + y3[i] + y4[i]);};

    *px = xc;
    *py = yc;

    return;}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_INTERPOL - given a mesh GRID, an array of points to interpolate to PTS,
 *             - and an array of functions to interpolate to the points FNCS
 *             - return an array of function values at the interpolation
 *             - points for each interpolation point we find the zone which
 *             - contains it by looping over zones and accepting the zone
 *             - for which the cross product of the vector from vertex to
 *             - point and vertex to next vertex is negative for each
 *             - vertex of the zone (note the vertices are ordered so as
 *             - to traverse the zone boundary in a counter-clockwise
 *             - fashion). FIND_VERTICES does this.
 *             - 
 *             - next, given the list of vertices defining the zone in which
 *             - the interpolation point resides, we find a set of weights for
 *             - the function values associated with the vertices so that the
 *             - value A0 at the interpolation point R0 is given by
 *             -
 *             -    A0 = sum(Ai x Wi)
 *             -
 *             - where the Ai are the values at the vertices and Wi are the 
 *             - weights. FIND_COEFFICIENTS does this.
 *             -
 *             - finally, given the set of weights which are only functions of
 *             - the geometry of the zones and are completely independent of
 *             - the function being interpolated, we interpolate each supplied
 *             - function to the interpolation point. INTERPOLATE_VALUE does
 *             - this.
 *             -
 *             - the rest of this routine handles the accessing and packaging
 *             - of the results.  this routine is completely general for 2d
 *             - meshes and as such is not particularly optimal for some
 *             - specific meshes.  it does best when given a large number of
 *             - functions to interpolate to each interpolation point.
 */

REAL **PM_interpol(grid, pts, n_pts, fncs, n_fncs)
   PM_lagrangian_mesh *grid;
   REAL **pts;
   int n_pts;
   REAL **fncs;
   int n_fncs;
   {int i, j;
    REAL *rix, *riy, **vals;
    coefficient *vertices;

/* allocate the return values */
    vals = FMAKE_N(REAL *, n_fncs, "PM_INTERPOL:vals");
    for (j = 0; j < n_fncs; j++)
        vals[j] = FMAKE_N(REAL, n_pts, "PM_INTERPOL:vals[]");     

/* get pointers to the points to interpolate to */
    rix = pts[0];
    riy = pts[1];

    vertices = PM_alloc_vertices(grid);
    for (i = 0; i < n_pts; i++)

/* find the vertices surrounding the IP */
        {PM_find_vertices(rix[i], riy[i], grid, vertices);

/* build the coefficients for the IP */
         PM_find_coefficients(rix[i], riy[i], grid, vertices);

/* interpolate all of the functions at the IP */
         for (j = 0; j < n_fncs; j++)
             vals[j][i] = PM_interpolate_value(vertices, fncs[j]);};

    return(vals);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_INSIDE - return TRUE if the point is "inside" the boundary segment
 *           -   (X, Y) is the interpolation point
 *           -   PX and PY are the arrays of mesh coordinates
 *           -   MAP is an array of indices into PX and PY which defines a
 *           -       zone (it's the list of vertices) the indices are ordered
 *           -       so as to traverse the boundary in a counter-clockwise
 *           -       direction
 *           -   N is the length of map, i.e. the number of vertices or sides
 *           -       to the zone
 *           -
 *           - the point is inside the zone if the cross product of the
 *           - vector from vertex to next vertex and vertex to point is
 *           - positive for each vertex of the zone
 *           - 
 *           - if any cross product is negative the point is outside the zone
 */

int PM_inside(x, y, px, py, map, n)
   double x, y;
   REAL *px, *py;
   int *map, n;
   {int ia, ib, ja, jb;
    REAL dxba, dyba, dx0a, dy0a, cross;

    for (ja = 0; ja < n; ja++)
        {jb = (ja + 1) % n;
         ia = map[ja];
         ib = map[jb];

         dxba = px[ib] - px[ia];
         dyba = py[ib] - py[ia];
         dx0a = x - px[ia];
         dy0a = y - py[ia];

         cross = (dxba*dy0a - dyba*dx0a);
         if (TOLERANCE > cross)
            return(FALSE);};

    return(TRUE);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_INSIDE_FIX - return TRUE if the point is "inside" the boundary
 *               -   (X, Y) is the interpolation point
 *               -   PX and PY are the arrays of boundary points
 *               -   N is the length of boundary
 *               -
 *               - the point is inside the zone if the cross product of the
 *               - vector from vertex to next vertex and vertex to point is
 *               - positive for each vertex of the boundary
 *               - 
 *               - if any cross product is negative the point is outside the
 *               - boundary
 */

int PM_inside_fix(x, y, px, py, n, direct)
   int x, y, *px, *py, n, direct;
   {int ja, jb;
    int dxba, dyba, dx0a, dy0a, cross;

    for (ja = 0; ja < n; ja++)
        {jb = (ja + 1) % n;

         dxba = px[jb] - px[ja];
         dyba = py[jb] - py[ja];
         dx0a = x - px[ja];
         dy0a = y - py[ja];

         cross = direct*(dxba*dy0a - dyba*dx0a);
         if (TOLERANCE > cross)
            return(FALSE);};

    return(TRUE);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_ALLOC_VERTICES - return an array containing the vertices */

coefficient *PM_alloc_vertices(grid)
   PM_lagrangian_mesh *grid;
   {int n_v;
    coefficient *vertices;

    grid = NULL;

    vertices = FMAKE(coefficient, "PM_ALLOC_VERTICES:vertices");

    vertices->n_points = n_v = 4;
    vertices->indexes  = FMAKE_N(int, n_v,
                         "PM_ALLOC_VERTICES:indexes");
    vertices->weights  = FMAKE_N(REAL, n_v,
                         "PM_ALLOC_VERTICES:weights");

    return(vertices);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_FIND_VERTICES - fill the array containing the vertices whose
 *                  - corresponding sides define the zone in which the
 *                  - given point resides
 *                  -
 *                  - it turns out that in this routine we assume that
 *                  - no connectivity information is supplied and proceed
 *                  - under the assumption that we have a logical
 *                  - rectangular mesh
 *                  -
 *                  - this being the case we compute the list of vertices
 *                  - (correctly ordered) for each zone as we go and test
 *                  - whether or not the point (X, Y) is inside
 *                  -
 *                  - when finished we will have a valid coefficient struct
 *                  - for the point (X, Y) which at this point will contain
 *                  - the vertices surrounding the interpolation point.
 */

int PM_find_vertices(x, y, grid, vertices)
   double x, y;
   PM_lagrangian_mesh *grid;
   coefficient *vertices;
   {int i, km, lm, n, n_v, *map;
    REAL *px, *py;

    px = grid->x;
    py = grid->y;
    km = grid->kmax;
    lm = grid->lmax;
    n  = km*lm;

    n_v = vertices->n_points;
    map = vertices->indexes;

    for (i = km+1; i < n; i++)
        {map[0] = i - km;
         map[1] = i;
         map[2] = i - 1;
         map[3] = i - km - 1;

         if (PM_inside(x, y, px, py, map, n_v))
            break;};

    return(TRUE);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_FIND_COEFFICIENTS - find the set of normalized weights which 
 *                      - specify the interpolating coefficients of the
 *                      - vertices
 *                      -
 *                      - this routine is entered with an interpolation
 *                      - point (X, Y), a mesh GRID, and a coefficient struct
 *                      - which contains the vertices surrounding the
 *                      - interpolation point
 *                      -
 *                      - when we leave the coefficient struct will also have
 *                      - a list of weights corresponding to the vertices
 *                      -
 *                      - for each vertex we compute the intersection of the
 *                      - ray from the vertex thru the interpolation point
 *                      - (IP) with one of the other zone boundary segments
 *                      - as defined by the other vertices.  we construct
 *                      - linear interpolation weights for this intersection
 *                      - point and then from those weights construct
 *                      - weights for the linear interpolation from the
 *                      - original vertex to the intersection point for
 *                      - the IP.  the partial weights are resolved into
 *                      - contributions from the three vertices involved
 *                      - and accumulated into values for all of the vertices.
 *                      - the contributions to the weights are further
 *                      - weighted by the inverse distance from the vertex
 *                      - to the IP.  this is done to insure that the
 *                      - contributions of the near vertices dominate over
 *                      - farther ones. consider the following situation
 *                      -
 *                      -                          .
 *                      -                         / \
 *                      -                        /   \
 *                      -                       /     \
 *                      -                      /       \
 *                      -      these should ->.    x    .
 *                      -      dominate        \       /
 *                      -                       \     /
 *                      -                        \   /
 *                      -                         \ /
 *                      -        not these ------> .
 *                      -
 *                      -
 *                      - when all of the vertices have been processed the
 *                      - weights are renormalized so that they sum to 1.
 */

int PM_find_coefficients(x, y, grid, vertices)
   double x, y;
   PM_lagrangian_mesh *grid;
   coefficient *vertices;
   {int ia, j, ja, ka, kb, la, lb, n;
    int i1, i2, i3, *map;
    REAL *px, *py;
    REAL xi, yi, dx0i, dy0i, d0i;
    REAL sx, sy, dsl;
    REAL dxj1j, dyj1j, dxij, dyij, cross, parallelp;
    REAL dsjx, dsjy, dsj, dj1j, dsign;
    REAL u, v, w;
    REAL *weights;

    px = grid->x;
    py = grid->y;
    i1 = 0;
    i2 = 0;
    i3 = 0;

    n       = vertices->n_points;
    map     = vertices->indexes;
    weights = vertices->weights;

    for (j = 0; j < n; j++)
        weights[j] = 0.0;

    for (ja = 0; ja < n; ja++)
        {ia = map[ja];

         xi = px[ia];
         yi = py[ia];

/* compute the distances from the vertices to the IP */
         dx0i = x - xi;
         dy0i = y - yi;
         d0i  = HYPOT(dx0i, dy0i);

/* find the side crossed by the line from the vertex to the IP */
         u  = HUGE;
         v  = HUGE;
         for (j = 1; j < n-1; j++)
             {ka = (ja + j) % n;
              kb = (ka + 1) % n;
              la = map[ka];
              lb = map[kb];

              dxj1j = px[lb] - px[la];
              dyj1j = py[lb] - py[la];

              parallelp = dx0i*dyj1j - dy0i*dxj1j;
              if (d0i < TOLERANCE)
                 {i1 = ja;
                  i2 = ka;
                  i3 = kb;
                  u = 1.0;
                  v = 0.0;}
              else if (ABS(parallelp) > TOLERANCE)
                 {dxij  = xi - px[la];
                  dyij  = yi - py[la];
                  cross = dxj1j*dyij - dyj1j*dxij;
                  sx    = xi + dx0i*cross/parallelp;
                  sy    = yi + dy0i*cross/parallelp;

                  dsl   = HYPOT(sx-xi, sy-yi);
                  dsjx  = sx - px[lb];
                  dsjy  = sy - py[lb];
                  dsj   = HYPOT(dsjx, dsjy);
                  dj1j  = HYPOT(dxj1j, dyj1j);
                  dsign = -(dsjx*dxj1j + dsjy*dyj1j);

/* dsl > d0i - says that the crossing is further away from xi than x0
 * dj1j >= dsj - d(Xj, Xj+1) >= d(cross, Xj)
 * dsign > 0 - says that the crossing is on the same side of Xj as Xj+1
 */
                  if ((dsl > d0i - TOLERANCE) &&
                      (dj1j >= dsj) && (dsign >= 0.0))
                     {i1 = ja;
                      i2 = ka;
                      i3 = kb;
                      u  = dsj/dj1j;
                      v  = d0i/dsl;
                      break;};};};

         if ((u > 1.0 + TOLERANCE) || (v > 1.0 + TOLERANCE))
            {io_printf(stderr, "INTERPOLATION ERROR AT (%11.3e, %11.3e)",
		       x, y);
             exit(3);};
         weights[i1] += 1.0 - v;
         weights[i2] += v*u;
         weights[i3] += v*(1.0 - u);};

/* renormalize the weights */
    w = 0.0;
    for (j = 0; j < n; j++)
        w += weights[j];

    w = 1.0/w;
    for (j = 0; j < n; j++)
        weights[j] *= w;

    return(TRUE);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_INTERPOLATE_VALUE - given the interpolating coefficients and the
 *                      - all of function values compute and return the
 *                      - value of the function at the interpolation point
 */

REAL PM_interpolate_value(vertices, f)
   coefficient *vertices;
   REAL *f;
   {int i, j, n, *map;
    REAL *weights, value;

    n       = vertices->n_points;
    map     = vertices->indexes;
    weights = vertices->weights;
    value   = 0.0;
    for (i = 0; i < n; i++)
        {j = map[i];
         value += (*(weights++))*f[j];};

    return(value);}

/*--------------------------------------------------------------------------*/

/*                       SPLINE INTERPOLATION ROUTINES                      */

/*--------------------------------------------------------------------------*/

/* _PM_SPLINE - set up cubic spline interpolation coefficients */

void _PM_spline(x, y, n, yp1, ypn, d2y)
   REAL *x, *y;
   double yp1, ypn;
   REAL *d2y;
   int n;
   {int i, k, nm;
    double p, qn, sig, un, *u;

    nm = n - 1;
    u  = FMAKE_N(double, nm, "_PM_SPLINE:u");
    if (yp1 == HUGE)
       d2y[0] = u[0] = 0.0;
    else
       {d2y[0] = -0.5;
	u[0]   = (3.0/(x[1] - x[0]))*((y[1] - y[0])/(x[1] - x[0]) - yp1);};

    for (i = 1; i < nm; i++)
        {sig    = (x[i] - x[i-1])/(x[i+1] - x[i-1]);
	 p      = sig*d2y[i-1] + 2.0;
	 d2y[i] = (sig - 1.0)/p;
	 u[i]   = (y[i+1] - y[i])/(x[i+1] - x[i]) -
	          (y[i] - y[i-1])/(x[i] - x[i-1]);
	 u[i]   = (6.0*u[i]/(x[i+1] - x[i-1]) - sig*u[i-1])/p;};

    if (ypn == HUGE)
       qn = un = 0.0;
    else
       {qn = 0.5;
	un = (3.0/(x[nm] - x[nm-1]))*
             (ypn - (y[nm] - y[nm-1])/(x[nm] - x[nm-1]));};

    d2y[nm] = (un - qn*u[nm-1])/(qn*d2y[nm-1] + 1.0);
    for (k = nm-1; k >= 0; k--)
        d2y[k] = d2y[k]*d2y[k+1] + u[k];

    SFREE(u);

    return;}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/

/* PM_CUBIC_SPLINE_INT - find the value of the function represented in
 *                     - FX and FY at the point X and return it in PY
 *                     - adapted from Numerical Recipes in C
 */

int PM_cubic_spline_int(fx, fy, d2y, n, x, py)
   REAL *fx, *fy, *d2y;
   double x;
   REAL *py;
   int n;
   {int k0, kn, k;
    double h, b, a;

/* find the appropriate bin */
    k0 = 0;
    kn = n - 1;
    while (kn-k0 > 1)
       {k = (kn + k0) >> 1;
	if (fx[k] > x)
	   kn = k;
	else
	   k0 = k;};

    h = fx[kn] - fx[k0];
    if (h == 0.0)
       return(FALSE);

    a = (fx[kn] - x)/h;
    b = (x - fx[k0])/h;

    h   = h*h/6.0;
    *py = a*(fy[k0] + (a*a - 1.0)*d2y[k0]*h) +
          b*(fy[kn] + (b*b - 1.0)*d2y[kn]*h);

    return(TRUE);}

/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/