/* taylor2.c   expand 2D taylor series, check accuracy and algorithm
 *
 *                 jn   kn  hx^j  hy^k  d(j)  d(k) 
 * F(x+hx,y+hy) = sum  sum  ----  ----  ----- ----- F(x,y)
 *                j=0  k=0   j!    k!   dx(j) dy(k) 
 *
 *         d(j)  d(k) 
 * denote  ----- ----- F(x,y) as  Fxxxyy   for j=3, k=2 
 *         dx(j) dy(k) 
 */
 
/* the first few terms are: (remember hy^0 = 1, 0! = 1)
 * F(x+hx,y+hy) = F                        ! j=0, k=0  thru nj=0, nk=0
 *
 *              + hx Fx                    ! j=1, k=0
 *              + hy Fy                    ! j=0, k=1
 *              + hx hy Fxy                ! j=1, k=1  thru nj=1, nk=1
 *
 *              + hx^2/2! Fxx              ! j=2, k=0
 *              + hx^2/2! hy/1!  Fxxy      ! j=2, k=1
 *              + hx^2/2! hy^2/2! Fxxyy    ! j=2, k=2
 *              + hy^2/2! Fyy              ! j=0, k=2
 *              + hx/1! hy^2/2! Fxyy       ! j=1, k=2  thru nj=2, nk=2
 *
 *              + hx^3/3! Fxxx             ! j=3, k=0
 *              + hx^3/3! hy/1! Fxxxy      ! j=3, k=1
 *              + hx^3/3! hy^2/2! Fxxxyy   ! j=3, k=2
 *              + hx^3/3! hy^3/3! Fxxxyyy  ! j=3, k=3
 *              + hy^3/3! Fyyy             ! j=0, k=3
 *              + hx/1! hy^3/3! Fxyyy      ! j=1, k=3
 *              + hx^2/2! hy^3/3! Fxxyyy   ! j=2, k=3  thru nj=3, nk=3
 *
 */

/* The above can be placed in a two dimensional matrix, A :
 *
 *          j=0        j=1         j=2        j=3        j=nj
 * k=0    1/(0! 0!)  1/(1! 0!)   1/(2! 0!)  1/(3! 0!)
 * k=1    1/(0! 1!)  1/(1! 1!)   1/(2! 1!)  1/(3! 1!)
 * k=2    1/(0! 2!)  1/(1! 2!)   1/(2! 2!)  1/(3! 2!)
 * k=3    1/(0! 3!)  1/(1! 3!)   1/(2! 3!)  1/(3! 3!)
 *
 * k=nk                                                 1/(nj! nk!)
 *
 * two dimensional matrix, B :
 *
 *         j=0        j=1         j=2        j=3       
 * k=0   hx^0 hy^0   hx^1 hy^0   hx^2 hy^0   hx^3 hy^0
 * k=1   hx^0 hy^1   hx^1 hy^1   hx^2 hy^1   hx^3 hy^1
 * k=2   hx^0 hy^2   hx^1 hy^2   hx^2 hy^2   hx^3 hy^2
 * k=3   hx^0 hy^3   hx^1 hy^3   hx^2 hy^3   hx^3 hy^3
 *
 * Numeric values evaluated at numeric (x,y) two dimensional matrix, B :
 *
 *         j=0        j=1         j=2        j=3       
 * k=0      F          Fx          Fxx        Fxxx
 * k=1      Fy         Fxy         Fxxy       Fxxxy
 * k=2      Fyy        Fxyy        Fxxyy      Fxxxyy
 * k=3      Fyyy       Fxyyy       Fxxyyy     Fxxxyyy 
 *
 * then the sum of the element by element product, A*B*C, is the Taylor series
 * expansion of F(x+hx,y+hy) in terms of F(x,y) and its derivatives at (x,y)
 * use upper left diagonal matrix for order O(h^n)
 *   for(j=0; j<n; j++) for(k=0; k+j<n; k++)  ABC[j][k]
 */

#include <stdio.h>
#include <math.h>
#include "deriv.h"
#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))

struct poly2{int nj;  /* data structure for two variable polynomial */
             int nk;
             double xy[20][20];}; /* should be [nj][nk] */

/* build A and B into AB */
static void taylor2_build(int nj, int nk, double hx, double hy,
                          struct poly2 *AB);

/* build just A, the factorials, into A */
static void taylor2_a(int nj, int nk, struct poly2 *A);

/* the A and B are computed, and f(x,y) fx(x,y) ... fxxxyyy(x,y) */ 
static double taylor2_eval(double x, double y, double hx, double hy,
                           struct poly2 P);

/* get the value of P1=ABC using only order n terms, nj>=n, nk>=n */
static double poly2_evaln(double x, double y, struct poly2 P1, int n);

/* evaluate P1=ABC at (x,y) */
static double poly2_eval (double x, double y, struct poly2 P1);

static void   poly2_print(struct poly2 P1);
static void   poly2_deriv(int jk, struct poly2 P1, struct poly2 *P2);
static void   poly2_clear(int nj, int nk, struct poly2 *P1);
static void   poly2_copy (struct poly2 P1, struct poly2 *P2);
static void   poly2_add  (struct poly2 P1, struct poly2 P2, struct poly2 *P3);
static void   poly2_sub  (struct poly2 P1, struct poly2 P2, struct poly2 *P3);
static void   poly2_mulc (double c, struct poly2 *P1);
static void   poly2_mul  (struct poly2 P1, struct poly2 P2, struct poly2 *P3);
static void   poly2_div  (double c, struct poly2 *P1);

/* value of partial derivative, P1 is for evaluating, not differentiating */
static double pderiv(int jorder, int korder, double x, double y,
                     struct poly2 P1);

int main(int argc, char * argv[]) /* test partial derivatives */
{
  struct poly2 P1; /* this is test polynomial */
  struct poly2 P2;
  struct poly2 P3;
  struct poly2 P4;
  struct poly2 P5;
  struct poly2 P6;
  struct poly2 P7;
  struct poly2 P8;
  struct poly2 P9;
  struct poly2 P10;
  /* F(x,y), Fx(x,y) ... is just a poly2 F */
  struct poly2 F; /* for taylor2_eval */
  /* full Taylor F + AB is FT */
  struct poly2 FT;
  struct poly2 A;
  struct poly2 AB;
  struct poly2 ABF;
  double x, y, hx, hy, val, val1;
  double xp, yp;
  double C[100];
  int i, j, k, n, nn, nj, nk;

  printf("taylor2.c running \n");
  n = 4;
  nn = (n+1)*(n+1); /* maximum constants */
  for(i=0; i<nn; i++)
  {
    C[i]= (double)(nn-i+1);
  }
  x = 1.0;
  y = 1.0;
  hx = 0.25;
  hy = 0.5;
  /* initialize P1 to test function */
  P1.nj = n;
  P1.nk = n;
  for(j=0; j<n; j++)
  {
    for(k=0; k<n; k++)
    {
      P1.xy[k][j]=C[((k+j)*(k+j)+k+3*j)/2]; /* has extra terms */
    }
  }

  poly2_print(P1);
  printf("trim \n");
  for(j=1; j<n; j++)
    for(k=n-1; k>=n-j; k--)
      P1.xy[k][j] = 0.0;
  /* P1 initialized */
  poly2_print(P1);
  val = poly2_eval(x, y, P1);
  printf("poly2_eval P1 at x=%g, y=%g, is %g \n",
           x, y, val);
  val = poly2_evaln(x, y, P1, n);
  printf("poly2_evaln P1 at x=%g, y=%g, is %g \n",
           x, y, val);
  xp = x + hx;
  yp = y + hy;
  val = poly2_eval(xp, yp, P1);
  printf("poly2_eval P1 at x+hx=%g, y+hy=%g, is %g \n",
           xp, yp, val);
  val = poly2_evaln(xp, yp, P1, n);
  printf("poly2_evaln P1 at x+hx=%g, y+hy=%g, is %g \n",
           xp, yp, val);

  val = pderiv(1, 0, x, y, P1);
  printf("derivative with respect to x = %g, using pderiv \n", val);
  poly2_deriv(0, P1, &P2);
  val = poly2_eval(x, y, P2);
  printf("derivative with respect to x = %g, is P2 \n", val);
  poly2_print(P2);

  val = pderiv(0, 1, x, y, P1);
  printf("derivative with respect to y = %g, using pderiv \n", val);
  poly2_deriv(1, P1, &P3);
  val = poly2_eval(x, y, P3);
  printf("derivative with respect to y = %g, is P3\n", val);
  poly2_print(P3);

  val = pderiv(2, 0, x, y, P1);
  printf("2nd derivative with respect to x = %g, using pderiv \n", val);
  poly2_deriv(0, P2, &P4);
  val = poly2_eval(x, y, P4);
  printf("2nd derivative with respect to x = %g, is P4\n", val);
  poly2_print(P4);

  val = pderiv(0, 2, x, y, P1);
  printf("2nd derivative with respect to y = %g, using pderiv \n", val);
  poly2_deriv(1, P3, &P5);
  val = poly2_eval(x, y, P5);
  printf("2nd derivative with respect to y = %g, is P5\n", val);
  poly2_print(P5);

  val = pderiv(1, 1, x, y, P1);
  printf("derivative with respect to x and y = %g, using pderiv \n", val);
  poly2_deriv(1, P2, &P6);
  val = poly2_eval(x, y, P6);
  printf("derivative with respect to x and y = %g, is P6\n", val);
  poly2_print(P6);

  val = pderiv(3, 0, x, y, P1);
  printf("3rd derivative with respect to x = %g, using pderiv \n", val);
  poly2_deriv(0, P4, &P7);
  val = poly2_eval(x, y, P7);
  printf("3rd derivative with respect to x = %g, is P7\n", val);
  poly2_print(P7);

  val = pderiv(0, 3, x, y, P1);
  printf("3rd derivative with respect to y = %g, using pderiv \n", val);
  poly2_deriv(1, P5, &P8);
  val = poly2_eval(x, y, P8);
  printf("3rd derivative with respect to y = %g, is P8\n", val);
  poly2_print(P8);

  val = pderiv(2, 1, x, y, P1);
  printf("derivative with respect to xx and y = %g, using pderiv \n", val);
  poly2_deriv(0, P6, &P9);
  val = poly2_eval(x, y, P9);
  printf("derivative with respect to xx and y = %g, is P9\n", val);
  poly2_print(P9);

  val = pderiv(1, 2, x, y, P1);
  printf("derivative with respect to x and yy = %g, using pderiv \n", val);
  poly2_deriv(1, P6, &P10);
  val = poly2_eval(x, y, P10);
  printf("derivative with respect to x and yy = %g, is P10\n", val);
  poly2_print(P10);

  printf("all fourth derivative matrices should be zero \n");
  printf("build test poly2 about (x,y) \n");

  taylor2_a(n, n, &A);
  printf("just factorials A \n");
  poly2_print(A);

  taylor2_build(n, n, hx, hy, &AB);
  printf("factorials and hx, hy  AB \n");
  poly2_print(AB);

  printf("build F manually, automatic in taylor2_eval \n");
  poly2_clear(n, n, &F);
  F.xy[0][0] = poly2_eval(x, y, P1); /* F(x,y) */
  F.xy[1][0] = poly2_eval(x, y, P2); /* Fx(x,y) */
  F.xy[0][1] = poly2_eval(x, y, P3); /* Fx(x,y) */
  F.xy[2][0] = poly2_eval(x, y, P4); /* Fxx(x,y) */
  F.xy[0][2] = poly2_eval(x, y, P5); /* Fyy(x,y) */
  F.xy[1][1] = poly2_eval(x, y, P6); /* Fxy(x,y) */
  F.xy[3][0] = poly2_eval(x, y, P7); /* Fxxx(x,y) */
  F.xy[0][3] = poly2_eval(x, y, P8); /* Fyyy(x,y) */
  F.xy[2][1] = poly2_eval(x, y, P9); /* Fxxy(x,y) */
  F.xy[1][2] = poly2_eval(x, y, P10); /* Fxyy(x,y) */

  printf("F(x,y), Fx(x,y) ... in  F \n");
  poly2_print(F);
  
  poly2_mul(AB, F, &ABF);
  printf("full matrix ABF \n");
  poly2_print(ABF);

  val = poly2_eval(x, y, ABF);
  val1 = poly2_eval(xp, yp, P1);
  printf("ABF(x,y,xh,yh) = %g same as poly2_eval(x+xh,y+yh) = %g \n\n",
         val, val1);

  val = taylor2_eval(x, y, hx, hy, P1); /* does work above */
  printf("taylor2_eval(x,y,xh,yh) = %g same as poly2_eval(x+xh,y+yh) = %g\n\n",
         val, val1);

  printf("taylor2.c finished \n");
  return 0;
} /* end main of taylor2.c  utility functions follow */

/* build A and B into T */
/* taylor2_build     expand two variable taylor series into matrix, no F */
static void taylor2_build(int nj, int nk, double hx, double hy,
                          struct poly2 *AB)
{
  int j, k; /* hx or hy may be positive, negative or zero */
  double pwrhx, pwrhy, jfac, kfac;
  
  AB->nj = nj;
  AB->nk = nk;
  pwrhx = 1.0;
  jfac  = 1.0;
  for(j=0; j<nj; j++)
  {
    pwrhy = 1.0;
    kfac  = 1.0;
    for(k=0; k<nk; k++)
    {
      AB->xy[j][k] = pwrhx*pwrhy/(jfac*kfac);
      pwrhy = pwrhy * hy;
      kfac = kfac * (k+1); /* for next time */
    }
    pwrhx = pwrhx * hx;
    jfac = jfac * (j+1);
  }
} /* end taylor2_build */

/* build just A, the factorials, into A */
/* taylor2_a  expand two variable taylor series into matrix, just factorials */
static void taylor2_a(int nj, int nk, struct poly2 *A)
{
  int j, k;
  double jfac, kfac;
  
  A->nj = nj;
  A->nk = nk;
  jfac  = 1.0;
  for(j=0; j<nj; j++)
  {
    kfac  = 1.0;
    for(k=0; k<nk; k++)
    {
      A->xy[j][k] = 1.0/(jfac*kfac);
      kfac = kfac * (k+1); /* for next time */
    }
    jfac = jfac * (j+1);
  }
} /* end taylor2_a */

/* the A and B are computed, and f(x,y) fx(x,y) ... fxxxyyy(x,y) */ 
/* taylor2_eval  evaluate two variable taylor series about x,y at x+hx,y+hy */
static double taylor2_eval(double x, double y, double hx, double hy,
                           struct poly2 P)
{
  int j, k, nj, nk; /* hx or hy may be positive, negative or zero */
  double pwrhx, pwrhy, jfac, kfac, val, dval, fval;
  struct poly2 D[10][10];
  
  nj = P.nj;
  nk = P.nk;
  if(nj>9 || nk>9)
  {
    printf("too big for taylor2_eval \n");
    return 0.0;
  }
  val = 0.0;
  pwrhx = 1.0;
  jfac = 1.0;
  poly2_copy(P, &D[0][0]);
  poly2_deriv(0, P, &D[1][0]);
  poly2_deriv(1, P, &D[0][1]);

  for(j=2; j<nj; j++)
  {
    poly2_deriv(0, D[j-1][0], &D[j][0]);
    poly2_deriv(1, D[0][j-1], &D[0][j]);
  }
  for(k=1; k<nk; k++)
  {
    for(j=1; j<nj; j++)
    {
      poly2_deriv(0, D[j-1][k], &D[j][k]);
    }
  }

  for(j=0; j<nj; j++)
  {
    pwrhy = 1.0;
    kfac  = 1.0;
    for(k=0; k<nk; k++)
    {
      dval = poly2_eval(x, y, D[j][k]); /* evaluate derivative at x,y */
      fval = pwrhx*pwrhy*dval/(jfac*kfac);
      val += fval;
      pwrhy = pwrhy * hy;
      kfac = kfac * (k+1); /* for next time */
    }
    pwrhx = pwrhx * hx;
    jfac = jfac * (j+1);
  }
  return val;
} /* end taylor2_eval */


/* get the value of P1=ABC using only order n terms, nj>=n, nk>=n */
static double poly2_evaln(double x, double y, struct poly2 P1, int n)
{
  double val;
  int j, k, nj, nk;
  double xpwr, ypwr;
  xpwr = 1.0;
  ypwr = 1.0;
  nj = P1.nj;
  nk = P1.nk;
  val = 0.0;
  if(nj<n || nk<n)
  {
    printf("poly2_evaln error n=%d, nj=%d, nk=%d \n", n, nj, nk);
    return val;
  }
  for(j=0; j<n; j++)
  {
    ypwr = 1.0;
    for(k=0; k+j<n; k++)
    {
      val = val + P1.xy[j][k]*xpwr*ypwr;
      ypwr = ypwr*y;
    }
    xpwr = xpwr*x;
  }
  return val;
} /* end poly2_evaln */

/* evaluate P1=ABC at (x,y) */
/* poly2_eval  evaluate a two variable polynomial, P1, at  x,y */
static double poly2_eval(double x, double y, struct poly2 P1)
{
  double val;
  int j, k, nj, nk;
  double xpwr, ypwr;
  xpwr = 1.0;
  ypwr = 1.0;
  nj = P1.nj;
  nk = P1.nk;
  val = 0.0;
  for(j=0; j<nj; j++)
  {
    ypwr = 1.0;
    for(k=0; k<nk; k++)
    {
      val = val + P1.xy[j][k]*xpwr*ypwr;
      ypwr = ypwr*y;
    }
    xpwr = xpwr*x;
  }
  return val;
} /* end poly2_eval */

/* poly2_print  print a two variable polynomial */
static void poly2_print(struct poly2 P1)
{
  int j, k, nj, nk;

  nj = P1.nj;
  nk = P1.nk;
  printf("poly2_print nj=%d, nk=%d \n", P1.nj, P1.nk);
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      printf("%9.6f  ", P1.xy[j][k]);  
    }
    printf("\n");
  }
  printf("\n");
} /* end poly2_print */


/* poly2_deriv derivative of two variable polynomial, jk=0 for first */
/*                                                    jk=1 for second */
static void poly2_deriv(int jk, struct poly2 P1, struct poly2 *P2)
{
  int j, k, nj, nk;
  
  if(jk==0) /* derivative with respect to x */
  {
    if(P1.nj>0) P2->nj = P1.nj - 1;
    else        P2->nj = 0;
    P2->nk = P1.nk;
    for(k=0; k<P1.nk; k++)
    {
      for(j=1; j<P1.nj; j++)
      {
        P2->xy[j-1][k] = P1.xy[j][k]*j;
      }
    }
  }
  else  /* derivative with respect to y */
  {
    P2->nj = P1.nj;
    if(P1.nk>0) P2->nk = P1.nk - 1;
    else        P2->nk = 0;
    for(j=0; j<P1.nj; j++)
    {
      for(k=1; k<P1.nk; k++)
      {
        P2->xy[j][k-1] = P1.xy[j][k]*k;
      }
    }
  }
} /* end poly2_deriv */

static void poly2_clear(int nj, int nk, struct poly2 *P1)
{
  int j, k;
  
  P1->nj = nj; /* also initializes, could use calloc */
  P1->nk = nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P1->xy[j][k] = 0.0;
    }
  }
} /* end poly2_clear */

static void   poly2_copy(struct poly2 P1, struct poly2 *P2)
{
  int j, k, nj, nk;
  
  nj = P1.nj;
  nk = P1.nk;
  P2->nj = nj;
  P2->nk = nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P2->xy[j][k] = P1.xy[j][k];
    }
  }
} /* end poly2_copy */

static void   poly2_add  (struct poly2 P1, struct poly2 P2, struct poly2 *P3)
{
  int j, k, nj, nk;
  
  if(P1.nj<P2.nj || P1.nk<P2.nk)
  {
    printf("poly2_add, P1 is smaller than P2 \n");
    return;
  }
  nj = P1.nj;
  nk = P1.nk;
  P3->nj = nj;
  P3->nk = nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P3->xy[j][k] = P1.xy[j][k];
    }
  }
  nj = P2.nj; /* may be smaller */
  nk = P2.nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P3->xy[j][k] += P2.xy[j][k];
    }
  }
} /* end poly2_add */

static void   poly2_sub  (struct poly2 P1, struct poly2 P2, struct poly2 *P3)
{
  int j, k, nj, nk;
  
  if(P1.nj<P2.nj || P1.nk<P2.nk)
  {
    printf("poly2_add, P1 is smaller than P2 \n");
    return;
  }
  nj = P1.nj;
  nk = P1.nk;
  P3->nj = nj;
  P3->nk = nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P3->xy[j][k] = P1.xy[j][k];
    }
  }
  nj = P2.nj; /* P2 may be smaller */
  nk = P2.nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P3->xy[j][k] -= P2.xy[j][k];
    }
  }
} /* end poly2_sub */

static void   poly2_mulc  (double c, struct poly2 *P1)
{
  int j, k, nj, nk;
  
  nj = P1->nj;
  nk = P1->nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P1->xy[j][k] *= c;
    }
  }
} /* end poly2_mul */

static void   poly2_mul  (struct poly2 P1, struct poly2 P2, struct poly2 *P3)
{
  int j, k, nj, nk;
  
  if(P1.nj!=P2.nj || P1.nk!=P2.nk)
  {
    printf("poly2_mul, P1 not same size as P2 \n");
    return;
  }
  nj = P1.nj;
  nk = P1.nk;
  P3->nj = nj;
  P3->nk = nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P3->xy[j][k] = P1.xy[j][k]*P2.xy[j][k];
    }
  }
} /* end poly2_mul */

static void   poly2_div  (double c, struct poly2 *P1)
{
  int j, k, nj, nk;
  
  nj = P1->nj;
  nk = P1->nk;
  for(j=0; j<nj; j++)
  {
    for(k=0; k<nk; k++)
    {
      P1->xy[j][k] /= c;
    }
  }
} /* end poly2_div */

static double pderiv(int jorder, int korder, double x, double y,
                     struct poly2 P1)
{
  double val, pval, hval, hx, hy;
  double xval[20];
  double denom;
  int npoints, point, ntry;
  int a[20];
  int i, j, bb;

  npoints = jorder+korder+4; /* ultra conservative */
  if((npoints/2)*2==npoints) npoints++;
  val = 0.0;
  point = (npoints-1)/2;
  hx = 0.01;
  hy = 0.01;
  ntry = 0;

  if(korder == 0)
  {
    deriv(jorder, npoints, point, a, &bb);
    hval = 0.0;
    for(i=0; i<npoints; i++)
      hval = hval + (double)a[i]*poly2_eval(x+(i-point)*hx, y, P1);
    hval = hval/((double)bb*pow(hx,jorder));
    hx = hx/2.0;
    val = 0.0;
    for(i=0; i<npoints; i++)
      val = val + (double)a[i]*poly2_eval(x+(i-point)*hx, y, P1);
    val = val/((double)bb*pow(hx,jorder));
    if(abs(val-hval)<0.00001) return val;
    npoints = npoints+2;
    point = (npoints-1)/2;
    deriv(jorder, npoints, point, a, &bb);
    pval = 0.0;
    for(i=0; i<npoints; i++)
      pval = pval + (double)a[i]*poly2_eval(x+(i-point)*hx, y, P1);
    pval = pval/((double)bb*pow(hx,jorder));
    if(abs(val-pval)<0.00001) return pval;
    printf("hval=%g, val=%g, pval=%g \n", hval, val, pval);
    return pval;
  }
  else if(jorder == 0)
  {
    deriv(korder, npoints, point, a, &bb);
    hval = 0.0;
    for(i=0; i<npoints; i++)
      hval = hval + (double)a[i]*poly2_eval(x, y+(i-point)*hy, P1);
    hval = hval/((double)bb*pow(hy,korder));
    hy = hy/2.0;
    val = 0.0;
    for(i=0; i<npoints; i++)
      val = val + (double)a[i]*poly2_eval(x, y+(i-point)*hy, P1);
    val = val/((double)bb*pow(hy,korder));
    if(abs(val-hval)<0.00001) return val;
    npoints = npoints+2;
    point = (npoints-1)/2;
    deriv(korder, npoints, point, a, &bb);
    pval = 0.0;
    for(i=0; i<npoints; i++)
      pval = pval + (double)a[i]*poly2_eval(x, y+(i-point)*hy, P1);
    pval = pval/((double)bb*pow(hx,korder));
    if(abs(val-pval)<0.00001) return pval;
    printf("hval=%g, val=%g, pval=%g \n", hval, val, pval);
    return pval;
  }
  else  /* must compute x derivative at npoints, then y derivative */
  {
    deriv(jorder, npoints, point, a, &bb);
    for(j=0; j<npoints; j++)
    {
      hval = 0.0;
      for(i=0; i<npoints; i++)
        hval = hval + (double)a[i]*
                      poly2_eval(x+(i-point)*hx, y+(j-point)*hy, P1);
      xval[j] = hval/((double)bb*pow(hx,jorder));
    }
    deriv(korder, npoints, point, a, &bb);
    val = 0.0;
    for(i=0; i<npoints; i++)
      val = val + (double)a[i]*xval[i];
    val = val/((double)bb*pow(hy,korder));
  }
  return val;
} /* end pderiv */


/* end taylor2.c */