// deriv2D.java  compute derivative D=U'(x,y) given U(x,y) and x,y arrays
//               uses nuderiv.java

class deriv2D
{
  deriv2D(int orderx, int ordery, int npx, int nx, double x[], int npy, int ny,
          double y[],  double U[][], double D[][])
  {
    double cx[] = new double[npx];
    double xg[] = new double[npx];
    double cy[] = new double[npy];
    double yg[] = new double[npy];
    double deriv;
    boolean debug1 = false;
    boolean debug2 = true;
    boolean debug = true;
    if(debug)
    {
      System.out.println("deriv2D orderx="+orderx+", npx="+npx+", nx="+nx);
      System.out.println("deriv2D ordery="+ordery+", npy="+npy+", ny="+ny);
      System.out.println("x[0]="+x[0]+", x[1]="+x[1]+", x[2]="+x[2]);
      System.out.println("y[0]="+y[0]+", y[1]="+y[1]+", y[2]="+y[2]);
      System.out.println("U[0][0]="+U[0][0]+", U[0][1]="+U[0][1]);
      System.out.println("U[1][0]="+U[1][0]+", U[1][1]="+U[1][1]);
    }
    if(ordery==0) // only x derivatives
    {
      for(int ii=0; ii<ny; ii++) // compute for all y
      {
        for(int j=0; j<npx; j++) xg[j] = x[j]; // xg used by nuderiv
        for(int i=0; i<(npx+1)/2; i++) // initial region
        {
          new nuderiv(orderx, npx, i, xg, cx);
          deriv = 0.0;
          for(int k=0; k<npx; k++) deriv += cx[k]*U[k][ii]; // raw k OK
          D[i][ii] = deriv; 
          if(debug1) System.out.println("i="+i+",ii="+ii+", D="+D[i][ii]);
        } // i

        for(int i=(npx+1)/2; i<nx-(npx/2); i++) // mid region
        {
          for(int j=0; j<npx; j++)
          {
            xg[j]=x[i+j-(npx+1)/2];
            if(debug1 && i==(npx+1)/2) System.out.println("j="+j+
                                                          ", xg="+xg[j]);
          }

          for(int j=0; j<npx; j++)
          {
	    new nuderiv(orderx, npx, (npx+1)/2, xg, cx);
            if(debug1 && j==0 && i<npx)
	    {
	      if(debug1) System.out.println("i="+i+", j="+j+
                                            ", pt="+((npx+1)/2));
	      for(int k=0; k<npx; k++)
	        if(debug1) System.out.println("k="+k+",cx[k]="+cx[k]+", U="+
                                          U[i+k-(npx-1)/2][ii]+", xg="+xg[k]);
	    }
            deriv = 0.0;
	    for(int k=0; k<npx; k++) deriv += cx[k]*U[i+k-(npx+1)/2][ii];
            D[i][ii] = deriv; 
          } // j
        } // i

        for(int j=0; j<npx; j++) xg[j] = x[nx-npx+j]; // xg used by nuderiv
        if(debug1)
          for(int j=0; j<npx; j++) System.out.println("ix="+(nx-npx+j)+", x="+
                                                 x[nx-npx+j]);
        for(int i=nx-(npx/2); i<nx; i++) // final region
        {
	  if(debug1) System.out.println("i="+i+", point="+(i-nx+npx));

          new nuderiv(orderx, npx, i-nx+npx, xg, cx);
          deriv = 0.0;
          for(int k=0; k<npx; k++) deriv += cx[k]*U[nx-npx+k][ii];
          D[i][ii] = deriv; 
          if(debug1) System.out.println("i="+i+", D="+D[i][ii]);
        } // i
      } // end ii for all y
    } // end if ordery==0

    if(orderx==0) // only y derivatives
    {
      for(int i=0; i<nx; i++)
      {
        for(int jj=0; jj<npy; jj++) yg[jj] = y[jj]; // yg used by nuderiv
        for(int ii=0; ii<(npy+1)/2; ii++) // initial region
        {
          new nuderiv(ordery, npy, ii, yg, cy);
          deriv = 0.0;
          for(int kk=0; kk<npy; kk++) deriv += cy[kk]*U[i][kk]; // raw k OK
          D[i][ii] = deriv; 
          if(debug1) System.out.println("i="+i+",ii="+ii+", D="+D[i][ii]);
        } // ii

        for(int ii=(npy+1)/2; ii<ny-(npy/2); ii++) // mid region
        {
          for(int jj=0; jj<npy; jj++)
          {
            yg[jj]=y[ii+jj-(npy+1)/2];
            if(debug1 && ii==(npy+1)/2) System.out.println("jj="+jj+", yg="+yg[jj]);
          }

          for(int jj=0; jj<npy; jj++)
          {
	    new nuderiv(ordery, npy, (npy+1)/2, yg, cy);
            if(debug1 && jj==0 && ii<npy)
	    {
	      if(debug1) System.out.println("ii="+ii+", jj="+jj+", pt="+
                                            ((npy+1)/2));
	      for(int kk=0; kk<npy; kk++)
	        if(debug1) System.out.println("kk="+kk+",cy[kk]="+cy[kk]+
                                 ", U="+U[i][ii+kk-(npy-1)/2]+", yg="+yg[kk]);
	    }
            deriv = 0.0;
	    for(int kk=0; kk<npy; kk++) deriv += cy[kk]*U[i][ii+kk-(npy+1)/2];
            D[i][ii] = deriv; 
          } // j
        } // i

        for(int jj=0; jj<npy; jj++) yg[jj] = y[ny-npy+jj]; //yg used by nuderiv
        if(debug1)
          for(int jj=0; jj<npy; jj++) System.out.println("iy="+(ny-npy+jj)+
                                       ", y="+y[ny-npy+jj]);
        for(int ii=ny-(npy/2); ii<ny; ii++) // final region
        {
	  if(debug1) System.out.println("ii="+ii+", point="+(ii-ny+npy));

          new nuderiv(ordery, npy, ii-ny+npy, yg, cy);
          deriv = 0.0;
          for(int kk=0; kk<npy; kk++) deriv += cy[kk]*U[i][ny-npy+kk];
          D[i][ii] = deriv; 
          if(debug1) System.out.println("ii="+ii+", D="+D[i][ii]);
        } // i
      } // end if nx==0
    } // end if orderx == 0

    // general case x and y derivatives
 
    for(int j=0; j<npx; j++) xg[j] = x[j]; // xg used by nuderiv
    for(int jj=0; jj<npy; jj++) yg[jj] = y[jj]; // yg used by nuderiv
    for(int i=0; i<(npx+1)/2; i++) // initial region
    {
      new nuderiv(orderx, npx, i, xg, cx);
      for(int ii=0; ii<(npy+1)/2; ii++)
      {
        new nuderiv(ordery, npy, ii, yg, cy);
        deriv = 0.0;
        for(int k=0; k<npx; k++)
	{
          for(int kk=0; kk<npy; kk++)
	  {
            deriv += cx[k]*cy[kk]*U[k][kk];
	  } // kk
	} // k
        D[i][ii] = deriv; 
        if(debug1) System.out.println("i="+i+",ii="+ii+", D="+D[i][ii]);
      } // ii
    } // i

    for(int i=(npx+1)/2; i<nx-(npx/2); i++) // mid region
    {
      for(int j=0; j<npx; j++)
      {
        xg[j]=x[i+j-(npx+1)/2];
      }
      new nuderiv(orderx, npx, (npx+1)/2, xg, cx);
      for(int ii=(npy+1)/2; ii<ny-(npy/2); ii++) // mid region
      {
        for(int jj=0; jj<npy; jj++)
        {
          yg[jj]=y[ii+jj-(npy+1)/2];
        }
	new nuderiv(ordery, npy, (npy+1)/2, yg, cy);
        deriv = 0.0;
        for(int k=0; k<npx; k++)
	{
          for(int kk=0; kk<npy; kk++)
	  {
            deriv += cx[k]*cy[kk]*U[i+k-(npx+1)/2][ii+kk-(npy+1)/2];
	  } // kk
	} // k
        D[i][ii] = deriv; 
      } // ii
    } // i

    for(int j=0; j<npx; j++) xg[j] = x[nx-npx+j]; // xg used by nuderiv
    for(int jj=0; jj<npy; jj++) yg[jj] = y[ny-npy+jj]; //yg used by nuderiv
    for(int i=nx-(npx/2); i<nx; i++) // final region
    {
      new nuderiv(orderx, npx, i-nx+npx, xg, cx);
      for(int ii=ny-(npy/2); ii<ny; ii++)
      {
        new nuderiv(ordery, npy, ii-ny+npy, yg, cy);
        deriv = 0.0;
        for(int k=0; k<npx; k++)
	{
          for(int kk=0; kk<npy; kk++)
	  {
            deriv += cx[k]*cy[kk]*U[nx-npx+k][ny-npy+kk];
	  } // kk
	} // k
        D[i][ii] = deriv; 
      } // ii
    } // i


  } // end deriv2D
} // end class deriv2D