// test_bihar3d.java  test setup for a specific Biharmonic PDE
//
// solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) + 
//        2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)
//
//  f(x,y,z) := .24
//
//  in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<xmax  ,
//
//  Dirichlet boundary conditions are computed using:
//  ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
//               0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 
//               0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1
//
//  Thus the derivatives in x and y are:
//
//  ux(x,y,z) := 0.04*x*x*x - 0.08*x*y*y - 0.12*z*z*x + 0.3 *x*x*y +
//               0.2*z*z*z - y - z
//
//  uy(x,y,z) := 0.08*y*y*y - 0.8*x*x*y - 0.1*y*z*z*z + 0.1*x*x*x +
//               0.9*y*y*z - x - z
//
//  uz(x,y,z) := 0.12*z*z*z - 0.1*y*y*z - 0.12*x*x*z +0.3*y*y*y +
//               0.6*x*z*z - y - x
//
//  continuing, for testing, thus additional derivatives are:
//  uxx(x,y,z) := 0.12*x*x -0.08*y*y - 0.12*z*z + 0.6*x*y
//  uxxx(x,y,z) := 0.24*x + 0.6*y
//  uxxxx(x,y,z) := 0.24
//  uyy(x,y,z) := 0.24*y*y -0.08*x*x - 0.1*z*z + 1.8*y*z
//  uyyy(x,y,z) := 0.48*y + 1.8*z
//  uyyyy(x,y,z) := 0.48
//  uzz(x,y,z) := 0.36*z*z - 0.10*y*y - 0.12*x*x + 1.2*x*z
//  uzzz(x,y,z) := 0.72*z + 1.2*x
//  uzzzz(x,y,z) := 0.72
//  uxxyy(x,y,z) := -0.16
//  uxxzz(x,y,z) := -0.24
//  uyyzz(x,y,z) := -0.20
//

import java.io.*;
import java.*;
import java.awt.*;
import java.awt.event.*;

public class test_bihar3d extends Frame
{
  // mouse and display variables
  int xm, ym, bm; // mouse position 
  int height = 700;
  int width = 700;
  double xmin = 0.0;
  double xmax = 2.0;
  double ymin = 0.0;
  double ymax = 2.0;
  double zmin = 0.0;
  double zmax = 2.0;
  double umin;
  double umax;
  double zs; // move plot
  double xw, yh, xoff, yoff;
  int nx = 9;
  int ny = 9;
  int nz = 9;
  double hx, hy, hz, hu;
  double ya = 0.5; //  orthographic typ 0.2 to 0.7, must fix text 
  double xs = (0.666*(width-100));
  double ys = (0.666*(height-100));
  int xp, yp; // for passing back ortho point

  public test_bihar3d()
  {
    System.out.println("test_bihar3d.java running.");
    System.out.println("solve  uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 2*uxxyy(x,y,z) + ");
    System.out.println("       2*uxxzz(x,y,z) + 2*uyyzz(x,y,z) = f(x,y,z)");
    System.out.println(" ");
    System.out.println("f(x,y,z) := .24");
    System.out.println(" ");
    System.out.println("in the cube   xmin<x<xmax and ymin<y<ymax and zmin<z<zmax  ,");
    System.out.println(" ");
    System.out.println("   subject to boundary conditions of the first kind,");
    System.out.println("   Dirichlet, must be specified at the boundary.");
    System.out.println("ub(x,y,z) := 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -");
    System.out.println("             0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y +");
    System.out.println("             0.3*y*y*y*z + 0.2*z*z*z*x - x*y - y*z - x*z + 1");
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax+
                       ", zmin="+zmin+", zmax="+zmax);

    setTitle("test_bihar3d");
    setSize(width,height);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    this.addMouseListener (new mousePressHandler());

    zs = zmax; // starting z for plot

    do_derivative_check();

    do_program();
  } // test_bihar3d constructor

  void do_derivative_check()
  {
    double x, y, z, err, maxerr, avgerr, rmserr;
    double sumerr, sumsqerr, deriv, eval_deriv;
    // for derivative check
    double xg[] = new double[nx];
    double yg[] = new double[ny];
    double zg[] = new double[nz];
    double ug[][][] = new double[nx][ny][nz];

    double cxx[][] = new double[nx][nx];
    double cxxxx[][] = new double[nx][nx];
    double cyy[][] = new double[ny][ny];
    double cyyyy[][] = new double[ny][ny];
    double czz[][] = new double[nz][nz];
    double czzzz[][] = new double[nz][nz];
    double tcxxyy[] = new double[ny]; // for partial uxxyy
    double tcxxzz[] = new double[nz]; // for partial uxxzz
    double tcyyzz[] = new double[nz]; // for partial uyyzz
    int i_max, j_max, k_max, cnt;

    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    hz = (zmax-zmin)/(double)(nz-1);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", hx="+hx+", hy="+hy+", hz="+hz);

    // check that analytic solution numerically satisfies PDE
    maxerr = 0.0;
    avgerr = 0.0;
    i_max = 0;
    j_max = 0;
    k_max = 0;
    for(int i=0; i<nx; i++)
    {
      x = xmin + (double)i*hx;
      xg[i] = x;
      for(int j=0; j<ny; j++)
      {
	y = ymin + (double)j*hy;
        yg[j] = y;
        for(int k=0; k<ny; k++)
        {
	  z = zmin + (double)k*hz;
          zg[k] = z;
          err = uxxxx(x,y,z) + uyyyy(x,y,z) + uzzzz(x,y,z) + 
	        2.0*uxxyy(x,y,z) + 2.0*uxxzz(x,y,z) + 2.0*uyyzz(x,y,z) - f(x,y,z);
	  ug[i][j][k] = ub(x,y,z);
          err = Math.abs(err); // PDE check against f(x,y,z)
          avgerr += err;
	  if(err>maxerr) 
	  {
	    maxerr = err;
	    i_max = i;
	    j_max = j;
            k_max = k;
	  }
	}
      }
    }
    x = xmin+i_max*hx;
    y = ymin+j_max*hy;
    z = zmin+k_max*hz;
    System.out.println("maxerr PDE-f(x,y,z) is "+maxerr+", at i_max="+
                       i_max+", j_max="+j_max+", k_max="+k_max);
    System.out.println("x="+x+", y="+y+", z="+z+", f(x,y,z)="+f(x,y,z));

    // check numerical accuracy of computation of derivatives
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    for(int i=1; i<nx-1; i++)
    {
      // compute derivative coefficients in x direction
      rderiv(2,nx,i,hx,cxx[i]);
      rderiv(4,nx,i,hx,cxxxx[i]);
    }
    for(int j=1; j<ny-1; j++)
    {
      // compute derivative coefficients in y direction
      rderiv(2,ny,j,hy,cyy[j]);
      rderiv(4,ny,j,hy,cyyyy[j]);
    }
    for(int k=1; k<nz-1; k++)
    {
      // compute derivative coefficients in z direction
      rderiv(2,nz,k,hz,czz[k]);
      rderiv(4,nz,k,hz,czzzz[k]);
    }

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxxxx
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  eval_deriv = 0.0; // numeric
          deriv = uxxxx(xg[i],yg[j],zg[k]); // analytic
          for(int m=0; m<nx; m++)
          {
	      eval_deriv += cxxxx[i][m]*ub(xg[m],yg[j],zg[k]);
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cxxxx  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uyyyy
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  eval_deriv = 0.0; // numeric
          deriv = uyyyy(xg[i],yg[j],zg[k]); // analytic
          for(int m=0; m<ny; m++)
          {
	      eval_deriv += cyyyy[j][m]*ub(xg[i],yg[m],zg[k]);
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cyyyy  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uzzzz
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  eval_deriv = 0.0; // numeric
          deriv = uzzzz(xg[i],yg[j],xg[k]); // analytic
          for(int m=0; m<ny; m++)
          {
	    eval_deriv += czzzz[j][m]*ub(xg[i],yg[j],zg[m]);
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cyyyy  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxxyy
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  deriv = uxxyy(xg[i],yg[j],zg[k]); // analytic
          for(int jj=0; jj<ny; jj++) // walk jj in y, computing uxx
	  {
	    tcxxyy[jj] = 0.0; // numeric  uxx part of uxxyy
            for(int m=0; m<nx; m++)
            {
	      tcxxyy[jj] += cxx[i][m]*ub(xg[m],yg[jj],zg[k]);
	    }
            // if(Math.abs(uxx(xg[i],yg[jj])-tcxxyy[jj])>0.001)
	    //   System.out.println("tcxxy["+jj+"]="+tcxxyy[jj]+
            //			  " != uxx="+uxx(xg[i],yg[jj]));
	  }
          eval_deriv = 0.0;
          for(int m=0; m<ny; m++)  // uyy part of uxxyy
          {
	    eval_deriv += cyy[j][m]*tcxxyy[m];
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cxxyy  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxxzz
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  deriv = uxxzz(xg[i],yg[j],zg[k]); // analytic
          for(int kk=0; kk<nz; kk++) // walk kk in z, computing uxx
	  {
	    tcxxzz[kk] = 0.0; // numeric  uxx part of uxxzz
            for(int m=0; m<nx; m++)
            {
	      tcxxzz[kk] += cxx[i][m]*ub(xg[m],yg[j],zg[kk]);
	    }
	  }
          eval_deriv = 0.0;
          for(int m=0; m<nz; m++)  // uzz part of uxxzz
          {
	    eval_deriv += czz[k][m]*tcxxzz[m];
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cxxzz  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uyyzz
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  deriv = uyyzz(xg[i],yg[j],zg[k]); // analytic
          for(int kk=0; kk<nz; kk++) // walk kk in z, computing uyy
	  {
	    tcyyzz[kk] = 0.0; // numeric  uyy part of uyyzz
            for(int m=0; m<ny; m++)
            {
	      tcyyzz[kk] += cyy[j][m]*ub(xg[i],yg[m],zg[kk]);
	    }
	  }
          eval_deriv = 0.0;
          for(int m=0; m<nz; m++)  // uzz part of uyyzz
          {
	    eval_deriv += czz[k][m]*tcyyzz[m];
	  }
          err = deriv - eval_deriv;
          maxerr = Math.max(maxerr,Math.abs(err));
	}
      }
    }
    System.out.println("check cyyzz  maxerr="+maxerr);

    maxerr = 0.0;
    sumerr = 0.0;
    sumsqerr = 0.0;
    for(int i=1; i<nx-1; i++) // full equation
    {
      x = xmin + (double)i*hx;
      for(int j=1; j<ny-1; j++)
      {
	y = ymin + (double)j*hy;
        for(int k=1; k<nz-1; k++)
        {
	  z = zmin + (double)k*hz;
	  eval_deriv = 0.0; // no PDE u term
        
          for(int m=0; m<nx; m++)
          {
	    eval_deriv += cxxxx[i][m]*ug[m][j][k];   // PDE  uxxxx term
	  } 
          for(int m=0; m<ny; m++)
          {
	    eval_deriv += cyyyy[j][m]*ug[i][m][k];    // PDE uyyyy term
	  }
          for(int m=0; m<nz; m++)
          {
	    eval_deriv += czzzz[k][m]*ug[i][j][m];    // PDE uyyyy term
	  }
          //  + 2*uxxyy
          for(int jj=0; jj<ny; jj++) // walk jj in y, computing uxx
	  {
	    tcxxyy[jj] = 0.0; // numeric  uxx part of uxxyy
            for(int m=0; m<nx; m++)
            {
	      tcxxyy[jj] += cxx[i][m]*ug[m][jj][k];
	    }
          }
          for(int m=0; m<ny; m++)  // uyy part of uxxyy
          {
	    eval_deriv += 2.0*cyy[j][m]*tcxxyy[m];  // PDE uxxyy term
	  }
          //  + 2*uxxzz
          for(int kk=0; kk<nz; kk++) // walk kk in z, computing uxx
	  {
	    tcxxzz[kk] = 0.0; // numeric  uxx part of uxxzz
            for(int m=0; m<nx; m++)
            {
	      tcxxzz[kk] += cxx[i][m]*ug[m][j][kk];
	    }
	  }
          for(int m=0; m<nx; m++)  // uzz part of uxxzz
          {
	    eval_deriv += 2.0*czz[k][m]*tcxxzz[m];
	  }
          //  + 2*uyyzz
          for(int kk=0; kk<nz; kk++) // walk kk in z, computing uyy
	  {
	    tcyyzz[kk] = 0.0; // numeric  uyy part of uyyzz
            for(int m=0; m<ny; m++)
            {
	      tcyyzz[kk] += cyy[j][m]*ug[i][m][kk];
	    }
	  }
          for(int m=0; m<nz; m++)  // uzz part of uyyzz
          {
	    eval_deriv += 2.0*czz[k][m]*tcyyzz[m];
	  }

          err = eval_deriv-f(xg[i],yg[j],zg[k]);
          sumerr += Math.abs(err);
          sumsqerr += err*err;
          maxerr = Math.max(Math.abs(err),maxerr);
          cnt++;
	} // end z
      } // end y
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("check all terms together in this PDE");
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+
                       ", hx="+hx+", hy="+hy+", hz="+hz);
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end do_derivative_check

  void do_program() // nx and ny can be changed
  {
    double x, y, z;

    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    hz = (zmax-zmin)/(double)(nz-1);
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz);
    System.out.println("hx="+hx+", hy="+hy+", hz="+hz);

    // for plotting
    xw = (xmax-xmin)*1.2;
    yh = (ymax-ymin)*1.2;
    xoff = 0.083333*xw;
    yoff = 0.083333*yh;

    umin = 100.0;
    umax = -100.0;
    for(int i=0; i<nx; i++)
    {
      x = xmin + (double)i*hx;
      for(int j=0; j<ny; j++)
      {
	y = ymin + (double)j*hy;
        for(int k=0; k<nz; k++)
        {
	  z = zmin + (double)k*hz;
          umax = Math.max(umax,u(x,y,z));
          umin = Math.min(umin,u(x,y,z));
	}
      }
    }
    System.out.println("umin="+umin+", umax="+umax);

    // plot
    repaint();

  } // end do_program

  class mousePressHandler extends MouseAdapter
  {
    public void mousePressed (MouseEvent e)
    {
      xm = e.getX();
      ym = e.getY();
      bm = e.getButton();
      // System.out.println("down xm="+xm+", ym="+ym+", bm="+bm); // debug

      if(in_rect(width-70, height-20, 30, 20))
      {
	// do something, run
	requestFocus();
        repaint();
      }
      else if(in_rect(width-105, height-20, 30, 20))
      {
	// Z up
	requestFocus();
        if(zs<zmax)
	{
	  zs = zs + hz;
	}
        do_program();
        repaint();
      }
      else if(in_rect(width-145, height-20, 35, 20))
      {
	// Z dn
	requestFocus();
        if(zs>zmin)
	{
	  zs = zs - hz;
	}
        do_program();
        repaint();
      }
      else if(in_rect(width-185, height-20, 30, 20))
      {
	// fine
	requestFocus();
	nx = nx+2;
	ny = ny+2;
        nz = nz+2;
        do_program();
        repaint();
      }
      else if(in_rect(width-225, height-20, 35, 20))
      {
	// big
	requestFocus();
        if(nx>5)
        {
	  nx = nx-2;
	  ny = ny-2;
          nz = nz-2;
	}
        do_program();
        repaint();
      }
      if(in_rect(width-35, height-20, 30, 20)) System.exit(0);
    } // end mousePressed 
  } // end mousePressHandler

  boolean in_rect(int x, int y, int w, int h)
  {
    return xm>x && xm<x+w && ym>y && ym<y+h;
  } // end in_rect

  boolean not_in(int i, int a[], int n)
  {
    boolean found = false;
    int j;

    for(j=0; j<n; j++)
    {
      if(i==a[j]) { found=true; break; }
    }
    return !found;
  } // end not_in

  double f(double x, double y, double z)
  {
      return 0.24;
  }

  double ub(double x, double y, double z)
  {
    return 0.01*x*x*x*x + 0.02*y*y*y*y + 0.03*z*z*z*z - 0.04*x*x*y*y -
	   0.05*y*y*z*z - 0.06*x*x*z*z + 0.1*x*x*x*y + 0.2*z*z*z*x -
           0.03*y*y*y*z - x*y - y*z - x*z + 1.0;
  }

  double u(double x, double y, double z)
  {
    return ub(x,y,z);
  }

  double uxxxx(double x, double y, double z)
  {
    return 0.24;
  }

  double uyyyy(double x, double y, double z)
  {
    return 0.48;
  }

  double uzzzz(double x, double y, double z)
  {
    return 0.72;
  }

  double uxxyy(double x, double y, double z)
  {
    return -0.16;
  }

  double uxxzz(double x, double y, double z)
  {
    return -0.24;
  }

  double uyyzz(double x, double y, double z)
  {
    return -0.20;
  }
 
  public void paint(Graphics g)
  {
    double x, y, u; // zs is changed by mouse
    double u1, u2, uz;
    int xp1, yp1, xp2, yp2; // line segment
    hx = (xmax-xmin)/nx; // for plotting with x+hx
    hy = (ymax-ymin)/ny;

    // buttons
    g.setColor(Color.red);
    g.fillRect(width-35,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("quit", width-30, height-10);
    g.setColor(Color.green);
    g.fillRect(width-70,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("run", width-65, height-10);
    g.setColor(Color.yellow);
    g.fillRect(width-105,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("Z up", width-100, height-10);
    g.setColor(Color.orange);
    g.fillRect(width-145,height-20,35,15);
    g.setColor(Color.black);
    g.drawString("Z dn", width-140, height-10);
    g.setColor(Color.yellow);
    g.fillRect(width-185,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("fine", width-180, height-10);
    g.setColor(Color.yellow);
    g.fillRect(width-225,height-20,35,15);
    g.setColor(Color.black);
    g.drawString("big ", width-220, height-10);

    // axis labels
    g.drawString("xmin", 25, height-15);
    g.drawString("xmax", (int)(0.66*width-70), height-15);
    g.drawString("ymin", (int)(0.66*width-30), height-25);
    g.drawString("ymax", (int)(width-80), (int)(0.66*height+25));
    g.drawString("umin", (int)(width-75), (int)(0.66*height+5));
    g.drawString("umax", (int)(width-75), (int)(85));

    // height color scale
    int hgt = 210;
    int ncol = 10;
    for(int i=0; i<ncol; i++)
    {
      set_color(i, g);
      g.fillRect(20,hgt-20*i,50,15);
    }
    g.drawString("u value", 80, 40);

    // draw box, bottom
    set_color(0, g);
    ortho(xmin, ymin, umin); xp1=xp; yp1=yp;
    ortho(xmax, ymin, umin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, umin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, umin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, umin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    set_color(ncol-1, g); // top
    ortho(xmin, ymin, umax); xp1=xp; yp1=yp;
    ortho(xmax, ymin, umax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, umax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, umax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, umax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    // draw box edges with u value colors
    u1 = umin;
    hu = (umax-umin)/(ncol);
    for(int i=0; i<ncol; i++)
    {
      u2 = u1+hu;
      set_color(i,g);
      ortho(xmin, ymin, u1); xp1=xp; yp1=yp;
      ortho(xmin, ymin, u2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmin, ymax, u1); xp1=xp; yp1=yp;
      ortho(xmin, ymax, u2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmax, ymax, u1); xp1=xp; yp1=yp;
      ortho(xmax, ymax, u2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmax, ymin, u1); xp1=xp; yp1=yp;
      ortho(xmax, ymin, u2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      u1 = u2;
    }

    // plot function or data with z value colors
    for(int i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(int j=0; j<ny; j++)
      {
	y = ymin+(double)j*hy;
        u1 = ub(x, y, zs);
        ortho(x, y, u1); xp1=xp; yp1=yp;
        u2 = ub(x+hx, y, zs);
        ortho(x+hx, y, u2); xp2=xp; yp2=yp;
        set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5), g);
        g.drawLine(xp1, yp1, xp2, yp2);
        u2 = ub(x, y+hy, zs);
        ortho(x, y+hy, u2); xp2=xp; yp2=yp;
        set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5), g);
        g.drawLine(xp1, yp1, xp2, yp2);
      }
    }
    for(int i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      y = ymax;
      u1 = ub(x, y, zs);
      ortho(x, y, u1); xp1=xp; yp1=yp;
      u2 = ub(x+hx, y, zs);
      ortho(x+hx, y, u2); xp2=xp; yp2=yp;
      set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5), g);
      g.drawLine(xp1, yp1, xp2, yp2);
    }
    for(int j=0; j<ny; j++)
    {
      y = ymin+(double)j*hy;
      x = xmax;
      u1 = ub(x, y, zs);
      ortho(x, y, u1); xp1=xp; yp1=yp;
      u2 = ub(x, y+hy, zs);
      ortho(x, y+hy, u2); xp2=xp; yp2=yp;
      set_color((int)((ncol*(((u1+u2)/2.0)-umin)/(umax-umin))+0.5), g);
      g.drawLine(xp1, yp1, xp2, yp2);
    }
  } // end paint

  // orthographic projection, sets xp,yp  "z" here may be u, etc
  void ortho(double x, double y, double u)
  {
    double x1, y1, u1, x2, y2, x3, y3;
    x1 = (x-xmin)/(xmax-xmin);   // zero base, normalized
    y1 = (y-ymin)/(ymax-ymin); 
    u1 = (u-umin)/(umax-umin);  
    x2 = x1+y1*ya; // ortho u is up, y is back
    y2 = u1+y1*ya;
    x3 = x2*xs;    // scale 
    y3 = y2*ys;
    x3 = x3+25;        // away from boarder
    y3 = height-y3-25; // away from boarder
    xp = (int)x3;
    yp = (int)y3;
    //System.out.println("ortho xp="+xp+", yp="+yp);
    //System.out.println("x1="+x1+", y1="+y1+", u1="+u1);
    //System.out.println(" ");
  } // end ortho

  void rderiv(int order, int npts, int point, double h, double c[])
  {
    nderiv D = new nderiv();
    int a[] = new int[npts]; // basic discretization from nderiv
    int b[] = new int[1];    // out parameter, denominator

    D.deriv(order,npts,point,a,b);
    for(int i=0; i<npts; i++)
    {
      c[i] = (double)a[i]/((double)b[0]);
      for(int j=0; j<order; j++) c[i] = c[i]/h;
    }
  } // end rderiv

  void set_color(int i, Graphics g)
  {
    if(i<0) {g.setColor(Color.black); return;}
    if(i>9) {g.setColor(Color.gray); return; }
    switch(i)
    {
      case 0: g.setColor(Color.black); break;
      case 1: g.setColor(Color.pink); break; // should be brown
      case 2: g.setColor(Color.red); break;
      case 3: g.setColor(Color.orange); break;
      case 4: g.setColor(Color.yellow); break;
      case 5: g.setColor(Color.green); break;
      case 6: g.setColor(Color.blue); break;
      case 7: g.setColor(Color.magenta); break;
      case 8: g.setColor(Color.cyan); break;
      case 9: g.setColor(Color.gray); break;
    }
  } // end set_color
 
  public static void main(String[] args)
  {
    new test_bihar3d();
  }
} // end class  test_bihar3d.java