// plot44.java  development test for perspective Utah graphics .dat
//                rotations about  dimensions X>, Y>, Z>, W>
//                color surfaces sorted X+, Y+, Z+, W+
//                needs datread44.java  Matrix.java
//
// perspective Viewing, Resize Choices and Transformation Matrices.
// makes a 2D screen look like a 4D sceen.
// In the simplest case, computing the i,j screen position of
// a vertex at world coordinates x,y,z,w is accomplished by
// multiplying a 5 by 5 perspective matrix by a 5 by 5 model matrix
// times the x,y,z,w,p vertex vector, then placing the i,j in
// the viewport.
//
// Name the vertices on the diagram as:
// On the 'near' rectangle, the lower left is  xmin,ymin,zmin,qmin
//                          the upper left is  'xmin' 'ymax'
//                          the lower right is 'xmax' 'ymin'
//                          the upper right is 'xmax' 'ymax'
// The distance from the eye to near is 'near' a positive number
// The distance from the eye to far  is 'far'  a positive number
//
// The 5 by 5 perspective matrix is
// 4D, four dimensional rendering, uses 5 by 5 matrices:
// Note that there are now eight (8) degrees of freedom:
// Move in X, Y, Z, W, P and rotations about each axis (split into 6 matrices)
// Notation: x is left and right.
//           y is up and down
//           z is forward and back
//           w is in and out (a fourth spacial dimension)
//           p = 1.0 point or vertex
// The 5 by 5 perspective matrix PM is
//
// |2*near        0.0          0.0       xmax+xmin        0.0      |
// |---------                            ---------                 |
// |xmax-xmin                            xmax-xmin                 |
// |                                                               |
// |  0.0        2*near        0.0       ymax+ymin        0.0      |
// |            --------                 ---------                 |
// |            ymax-ymin                ymax-ymin                 |
// |                                                               |
// |  0.0         0.0         2*near     zmax+zmin        0.0      |
// |                         -------     -----------               |
// |                         zmax-zmin   zmax-zmin                 |
// |                                                               |
// |  0.0         0.0          0.0       -(far+near)  -2*far*near  |
// |                                     -----------  -----------  |
// |                                      far-near     far-near    |
// |                                                               |
// |  0.0         0.0          0.0          -1.0          0.0      |
//
//  PM = makePM(xmin, xmax, ymin, ymax, zmin, zmax, near, far)
//
// The model view matrix is the product of the needed matrices below.
//
// The translation matrix TM, to translate 0,0,0,0 to x,y,z,w is
// | 1.0  0.0  0.0  0.0   fx  |
// | 0.0  1.0  0.0  0.0   fy  |    unused translations are 0.0
// | 0.0  0.0  1.0  0.0   fz  |
// | 0.0  0.0  0.0  1.0   fw  |
// | 0.0  0.0  0.0  0.0   1.0 |
//
//  TM = makeTM(fx,fy,fz,fw)
//  new (x, y, z, w, p) = TM * (x, y, z, w, p) 
//
//
// The scaling matrix SM, to scale x by sx, y by sy, z by sz, w by sw is
// |  sx  0.0  0.0  0.0  0.0 |
// | 0.0   sy  0.0  0.0  0.0 |    unused scales are 1.0
// | 0.0  0.0   sz  0.0  0.0 |
// | 0.0  0.0  0.0   sw  0.0 |
// | 0.0  0.0  0.0  0.0  1.0 |
//
//  SM = make_SM(sx, sy, sz, sw)
//  new (x, y, z, w, p) =  SM * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the X,Y axis is
// | 1.0  0.0  0.0    0.0     0.0 |
// | 0.0  1.0  0.0    0.0     0.0 |
// | 0.0  0.0  cos a  -sin a  0.0 |
// | 0.0  0.0  sin a  cos a   0.0 |
// | 0.0  0.0  0.0    0.0     1.0 |
//
//  XY = makeXY(a)
//  new (x, y, z, w, p) = XY * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the X,Z axis is
// | 1.0  0.0    0.0  0.0     0.0 |
// | 0.0  cos a  0.0  -sin a  0.0 |
// | 0.0  0.0    1.0  0.0     0.0 |
// | 0.0  sin a  0.0  cos a   0.0 |
// | 0.0  0.0    0.0  0.0     1.0 |
//
//  XZ = makeXZ(a)
//  new (x, y, z, w, p) = XZ * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the X,W axis is
// | 1.0    0.0    0.0     0.0  0.0 |
// | 0.0    cos a  -sin a  0.0  0.0 |
// | 0.0    sin a  cos a   0.0  0.0 |
// | 0.0    0.0    0.0     1.0  0.0 |
// | 0.0    0.0    0.0     0.0  1.0 |
//
//  XW = makeXW(a)
//  new (x, y, z, w, p) = XW * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the Y,Z axis is
// | cos a   0.0   0.0   -sin a  0.0 |
// | 0.0     1.0   0.0   0.0     0.0 |
// | 0.0     0.0   1.0   0.0     0.0 |
// | sin a   0.0   0.0   cos a   0.0 |
// | 0.0     0.0   0.0   0.0     1.0 |
//
//  YZ = makeYZ(a)
//  new (x, y, z, w, p) = YZ * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the Y,W axis is
// | cos a   0.0   -sin a  0.0  0.0 |
// | 0.0     1.0   0.0     0.0  0.0 |
// | sin a   0.0   cos a   0.0  0.0 |
// | 0.0     0.0   0.0     1.0  0.0 |
// | 0.0     0.0   0.0     0.0  1.0 |
//
//  YW = makeYW(a)
//  new (x, y, z, w, p) = YW * (x, y, z, w, p)
//
//
// The rotation matrix by angle a about the Z,W axis is
// | cos a  -sin a  0.0    0.0  0.0 |
// | sin a  cos a   0.0    0.0  0.0 |
// | 0.0    0.0     1.0    0.0  0.0 |
// | 0.0    0.0     0.0    1.0  0.0 |
// | 0.0    0.0     0.0    0.0  1.0 |
//
//   ZW = makeZW(a)
//   new (x, y, z, w, p) = ZW * (x, y, z, w, p)
//
//
// To get a rotation about only the X axis by ax,
// use the matrix product of  X,Y  X,Z  X,W = XR
//  new (x, y, z, w, p) = XY(ax) * XZ(ax) * XW(ax) * (x, y, z, w, p)
//
// To get a rotation about only the Y axis by ay,
// use the matrix product of  X,Y  Y,Z  Y,W = YR
//  new (x, y, z, w, p) = XY(ay) * YZ(ay) * YW(ay) * (x, y, z, w, p)
//
// To get a rotation about only the Z axis by az,
// use the matrix product of  X,Z  Y,Z  Z,W = ZR
//  new (x, y, z, w, p) = XZ(az) * YZ(az) * ZW(az) * (x, y, z, w, p)
//
// To get a rotation about only the W axis,
// use the matrix product of  X,W  Y,W  Z,W =WR
//  new (x, y, z, w, p) = XW(aw) * YW(aw) * ZW(aw) * (x, y, z, w, p)
//
//
// A user world coordinate vertex p = x, y, z, w, p  (p=1.0)
// is transformed into  pp  by
//
// perspective matrix times selection of model view matrix times p is pp
//
// To get screen coordinates, given the screen width w, and
// screen height h,
//
// screen x = w * ((pp.x/pp.w)-xmin)/(xmax-xmin) ?
// screen y = h * ((pp.y/pp.w)-ymin)/(ymax-ymin) ?

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

public class plot44 extends Frame
{
  double PM[][] = new double[5][5]; // perspective matrix
  double TM[][] = new double[5][5]; // translation matrix
  double SM[][] = new double[5][5]; // scale matrix
  double XY[][] = new double[5][5]; // XY rotate matrix
  double XZ[][] = new double[5][5]; // XZ rotate matrix
  double XW[][] = new double[5][5]; // XW rotate matrix
  double YZ[][] = new double[5][5]; // YZ rotate matrix
  double YW[][] = new double[5][5]; // YW rotate matrix
  double ZW[][] = new double[5][5]; // ZW rotate matrix
  double XR[][] = new double[5][5]; // X rotate matrix  ax, -ax
  double XRT[][] = new double[5][5]; // X rotate matrix transpose
  double XRI[][] = new double[5][5]; // X rotate matrix inverse
  double YR[][] = new double[5][5]; // Y rotate matrix  ay, -ay
  double ZR[][] = new double[5][5]; // Z rotate matrix  az, -az
  double WR[][] = new double[5][5]; // W rotate matrix  aw, -aw
  double XYZWR[][] = new double[5][5]; // 4 angle rotate matrix
  double fx = 1.0; // offset for test, reset after test
  double fy = 2.0;
  double fz = 3.0;
  double fw = 4.0;
  double off = 0.0;
  double offx = 0.0;
  double offy = 0.0;
  double offz = 0.0;
  double offw = 0.0;
  double sx = 1.25; // scale for test, reset after test
  double sy = 1.33;
  double sz = 1.5;
  double sw = 1.7;
  double sca = 1.0;
  double scax = 1.0;
  double scay = 1.0;
  double scaz = 1.0;
  double scaw = 1.0;
  double Pi = Math.PI;
  double ax = (15.0/180.0)*Pi; // angle for testing, reset after test 
  double ay = (30.0/180.0)*Pi;
  double az = (45.0/180.0)*Pi;
  double aw = (60.0/180.0)*Pi;
  double ang = 0.0;
  double angx = 0.0; // for XYZWR
  double angy = 0.0;
  double angz = 0.0;
  double angw = 0.0;
  double p1234[] = {1.0, 2.0, 3.0, 4.0, 1.0}; // for testing
  double p1111[] = {1.0, 1.0, 1.0, 1.0, 1.0}; // for testing
  double p0000[] = {0.0, 0.0, 0.0, 0.0, 1.0}; // for testing
  double pxi[]   = {1.0, 0.0, 0.0, 0.0, 1.0};   
  double pyi[]   = {0.0, 1.0, 0.0, 0.0, 1.0};   
  double pzi[]   = {0.0, 0.0, 1.0, 0.0, 1.0};   
  double pwi[]   = {0.0, 0.0, 0.0, 1.0, 1.0};   
  double pnxi[]   = {-1.0, 0.0, 0.0, 0.0, 1.0};   
  double pnyi[]   = {0.0, -1.0, 0.0, 0.0, 1.0};   
  double pnzi[]   = {0.0, 0.0, -1.0, 0.0, 1.0};   
  double pnwi[]   = {0.0, 0.0, 0.0, -1.0, 1.0};   
  double p[]   = new double[5]; // temp for testing   
  double p1[]  = new double[5];   
  double p2[]  = new double[5];   
  double p3[]  = new double[5];
  double pr[]  = new double[5];   
  double M[][] = new double[5][5]; // temp
  double M1[][] = new double[5][5];
  double xmin, xmax, ymin, ymax, zmin, zmax, wmin, wmax;
  double near, far;
  double a, b, c, d; // for testing
    
  // mouse and display variables
  boolean debug = false;
  boolean debug1 = true; // set true to check matrix operations
  boolean debug2 = false; // set true to check matrix operations
  int xm, ym, bm; // mouse position 
  int height = 700;
  int width = 700; // must be same as height for ortho
  int xp, yp; // ortho returns
  int pts, pt, k;
  datread44 DR; // = new datread44(filename); // defines data structures
  String filein;
  int grid = 1; // box off/on
  int colora = 1; // black, X, Y, blue Z, red xW
  int points = 1; // lines, points only, both
  boolean fillc = false;
  double ss, za, ws, wf; // used in ortho
  String lcol[] = {" ", "x", "y", "z", "w"};
  double xyzw[][] = new double[2000][5]; // DR.data_points for rotation
  int num_points = 0; // DR.num_points for rotation
  int num_polys = 0;  // DR.num_polys for rotation
  DecimalFormat f1 = new DecimalFormat("0.00000");
  double nx, ny, nz, nw; // global for normal4
  float red, grn, blu;
  int pripoly[] = new int[2000]; // for sort
  int priindex[] = new int[2000];
  int kpoly;
  double prix[] = new double[2000]; // used for sort
  double priy[] = new double[2000];
  double priz[] = new double[2000];
  double priw[] = new double[2000];
  double avgx, avgy, avgz, avgw;
  boolean xplus = false;
  boolean yplus = false;
  boolean zplus = false;
  boolean wplus = false;
  boolean xsort = false;
  boolean ysort = false;
  boolean zsort = false;
  boolean wsort = false;
    
  public plot44(String filename)
  {
    System.out.println("plot44.java running");
    setTitle("plot44");
    setSize(width,height);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
	System.out.println("plot44.java window closed.");
        System.exit(0);
      }
    });
    setVisible(true);
    this.addMouseListener (new mousePressHandler());

    System.out.println("plot44.java rotate test running");
    System.out.println("test matrix functions");
    // test_matrix();
    ax = (15.0/180.0)*Pi; 
    ay = (15.0/180.0)*Pi;
    az = (15.0/180.0)*Pi;
    aw = (15.0/180.0)*Pi;
    angx = 0.0;
    angy = 0.0;
    angz = 0.0;
    angw = 0.0;
    scax = 1.1;
    scay = 1.1;
    scaz = 1.1;
    scaw = 1.1;
    offx = 0.1;
    offy = 0.1;
    offz = 0.1;
    offw = 0.1;
    // System.out.println("datread44.java reading "+filename);
    filein = filename;
    DR = new datread44(filename); // defines data structures
    num_points = DR.num_points;
    num_polys = DR.num_polys;
    System.out.println("num_points="+num_points+
                       ", num_polys="+num_polys);
    System.out.println("print vertices");
    xmin = DR.data_points[0].x;
    xmax = DR.data_points[0].x;
    ymin = DR.data_points[0].y;
    ymax = DR.data_points[0].y;
    zmin = DR.data_points[0].z;
    zmax = DR.data_points[0].z;
    wmin = DR.data_points[0].w;
    wmax = DR.data_points[0].w;
    for(int i=0; i<num_points; i++) // vertices
    {
      System.out.println(+(i+1)+"  "+DR.data_points[i].x+
                         "  "+DR.data_points[i].y+
			 "  "+DR.data_points[i].z+
			 "  "+DR.data_points[i].w);
      xyzw[i][0] = DR.data_points[i].x;
      xyzw[i][1] = DR.data_points[i].y;
      xyzw[i][2] = DR.data_points[i].z;
      xyzw[i][3] = DR.data_points[i].w;
      xyzw[i][4] = 1.0; // for rotation
      xmin = Math.min(xmin, DR.data_points[i].x);
      xmax = Math.max(xmax, DR.data_points[i].x);
      ymin = Math.min(ymin, DR.data_points[i].y);
      ymax = Math.max(ymax, DR.data_points[i].y);
      zmin = Math.min(zmin, DR.data_points[i].z);
      zmax = Math.max(zmax, DR.data_points[i].z);
      wmin = Math.min(wmin, DR.data_points[i].w);
      wmax = Math.max(wmax, DR.data_points[i].w);
    } // end print vertices

    System.out.println("print polygons");
    k = 0; // polygons
    for(int i=0; i<num_polys; i++)
    {
      pts = DR.data_polys[k++];
      System.out.print(+pts); // no cr
      for(int j=0; j<pts; j++)
      {
        pt = DR.data_polys[k++];
        System.out.print(" "+ pt);
      }
      System.out.println(" ");
    } // end print polygons

  } // plot44 constructor


  class mousePressHandler extends MouseAdapter
  {
    public void mousePressed (MouseEvent e)
    {
      int oldxm, oldym;
      xm = e.getX();
      ym = e.getY();
      bm = e.getButton();
      if(debug) System.out.println("down xm="+xm+", ym="+ym+", bm="+bm);
      if(in_rect(width-100, height-40, 40, 15)) // bigger is down
      {
	// color points and lines
	colora = 1;
      }
      if(in_rect(width-100, height-20, 40, 15)) // bigger is down
      {
	// no color points and lines
	colora = 0;
      }
      else if(in_rect(width-145, height-40, 35, 15))
      {
	// show  min, max  boundaries
	grid = 1;
      }
      else if(in_rect(width-145, height-20, 35, 15))
      {
	// no show  min, max  boundaries
	grid = 0;
      }
      else if(in_rect(width-185, height-40, 30, 15))
      {
	// show  points
        points = 1;
      }
      else if(in_rect(width-185, height-20, 30, 15))
      {
	// show  points
        points = 0;	
      }
      else if(in_rect(width-50, height-20, 40, 15))
      {
	System.out.println("plot44.java quit.");
        System.exit(0);
      }
      else if(in_rect(width-50, height-40, 40, 15)) // init original
      {
        for(int i=0; i<num_points; i++) // vertices
        {
          xyzw[i][0] = DR.data_points[i].x;
          xyzw[i][1] = DR.data_points[i].y;
          xyzw[i][2] = DR.data_points[i].z;
          xyzw[i][3] = DR.data_points[i].w;
          xyzw[i][4] = 1.0; // for rotation
	}
	maxmin(xyzw);
	angx = 0.0;
	angy = 0.0;
	angz = 0.0;
	angw = 0.0;
	fillc = false;
      }
      // rotations < and >
      else if(in_rect(width-170, height-60, 30, 15)) // <
      {
	System.out.println("<X");
	angx -= ax;
	makeXR(-angx,XR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(XR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-135, height-60, 30, 15))
      {
	System.out.println("<Y");
	angy -= ay;
	makeYR(-angy,YR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(YR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-100, height-60, 30, 15))
      {
	System.out.println("<Z");
	angz -= az;
	makeZR(-angz,ZR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(ZR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-65, height-60, 30, 15))
      {
	System.out.println("<W");
	angw -= aw;
	makeWR(-angw,WR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(WR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-30, height-60, 30, 15))
      {
	System.out.println("<R");
	angx -= ax;
	angy -= ay;
	angz -= az;
	angw -= aw;
	makeXYZWR(-angx, -angy, -angz, -angw, XYZWR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(XYZWR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-170, height-80, 30, 15)) // >
      {
	if(debug) System.out.println("X>");
	angx += ax;
	makeXR(angx,XR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(XR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-135, height-80, 30, 15))
      {
	if(debug) System.out.println("Y>");
	angy += ay;
	makeYR(angy,YR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(YR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-100, height-80, 30, 15))
      {
	if(debug) System.out.println("Z>");
	angz += az;
	makeZR(angz,ZR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(ZR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-65, height-80, 30, 15))
      {
	if(debug) System.out.println("W>");
	angw += aw;
	makeWR(angw,WR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(WR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-30, height-80, 30, 15))
      {
	if(debug) System.out.println("R>");
	angx += ax;
	angy += ay;
	angz += az;
	angw += aw;
	makeXYZWR(angx, angy, angz, angw, XYZWR);
	for(int i=0; i<num_points; i++)
	{
	  getpr(i, pr);
	  matvec(XYZWR, pr, p);
	  veccopy(p, xyzw[i]);
	}
	maxmin(xyzw);
      }
      else if(in_rect(width-30, height-120, 30, 15))
      {
	System.out.println("fillc");
        fillc = true;
      }
      else if(in_rect(width-30, height-100, 30, 15))
      {
	System.out.println("no fillc");
        fillc = false;
      }
      else if(in_rect(width-170, height-120, 30, 15)) // >
      {
	if(debug) System.out.println("X+");
        xsort = true;
        ysort = false;
        zsort = false;
        wsort = false;
      }
      else if(in_rect(width-135, height-120, 30, 15))
      {
	if(debug) System.out.println("Y+");
        ysort = true;
        xsort = false;
        zsort = false;
        wsort = false;
      }
      else if(in_rect(width-100, height-120, 30, 15))
      {
	if(debug) System.out.println("Z+");
        zsort = true;
        xsort = false;
        ysort = false;
        wsort = false;
      }
      else if(in_rect(width-65, height-120, 30, 15))
      {
	if(debug) System.out.println("W+");
        wsort = true;
        xsort = false;
        ysort = false;
        zsort = false;
      }
      else if(in_rect(width-170, height-100, 30, 15)) // >
      {
	if(debug) System.out.println("X+");
        xplus = true;
      }
      else if(in_rect(width-135, height-100, 30, 15))
      {
	if(debug) System.out.println("Y+");
        yplus = true;
      }
      else if(in_rect(width-100, height-100, 30, 15))
      {
	if(debug) System.out.println("Z+");
        zplus = true;
      }
      else if(in_rect(width-65, height-100, 30, 15))
      {
	if(debug) System.out.println("W+");
        wplus = true;
      }
      else
      { // do nothing
        return;
      }
      xm = -1; // prevent reclick
      ym = -1;
      requestFocus();
      repaint();
    } // 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;
  }
 
  public void paint(Graphics g)
  {
    double x, y, z, w;
    double x0=0.0, y0=0.0, z0=0.0, w0=0.0;
    double z1=0.0, z2=0.0, x1=0.0, x2=0.0, y1=0.0, y2=0.0, hx, hy, hz, hw;
    double w1=0.0, w2=0.0;
    int xp1=0, yp1=0, xp2=0, yp2=0; // line segment
    String fstring;
    int hgt = 240;
    int ncol = 4;
    int xvec[] = new int[12]; // for g.fillPolygon  after set color
    int yvec[] = new int[12]; // screen coordinates
    int npp = 0;                    // number points in polygon
    double P1[] = new double[4]; // for normal4, assume plane
    double P2[] = new double[4];
    double P3[] = new double[4];

    // buttons
    g.setColor(Color.red);
    g.fillRect(width-55,height-20,50,15);
    g.setColor(Color.black);
    g.drawString("quit", width-50, height-10);

    g.setColor(Color.green);
    g.fillRect(width-55,height-40,50,15);
    g.setColor(Color.black);
    g.drawString("initial", width-50, height-30);

    g.setColor(Color.green);
    g.fillRect(width-100,height-40,40,15);
    g.setColor(Color.black);
    g.drawString("color", width-95, height-30);

    g.setColor(Color.red);
    g.fillRect(width-100,height-20,40,15);
    g.setColor(Color.black);
    g.drawString("no", width-95, height-10);
    g.setColor(Color.green);
    g.fillRect(width-145,height-40,35,15);
    g.setColor(Color.black);
    g.drawString("grid", width-140, height-30);
    g.setColor(Color.red);
    g.fillRect(width-145,height-20,35,15);
    g.setColor(Color.black);
    g.drawString("no", width-140, height-10);
    g.setColor(Color.green);
    g.fillRect(width-185,height-40,30,15);
    g.setColor(Color.black);
    g.drawString("pts", width-180, height-30);
    g.setColor(Color.red);
    g.fillRect(width-185,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("no", width-180, height-10);

    g.setColor(Color.green);
    g.fillRect(width-170,height-80,30,15);
    g.setColor(Color.black);
    g.drawString("XR>", width-167, height-67);
    g.setColor(Color.red);
    g.fillRect(width-170,height-60,30,15);
    g.setColor(Color.black);
    g.drawString("<XR", width-167, height-47);

    g.setColor(Color.green);
    g.fillRect(width-135,height-80,30,15);
    g.setColor(Color.black);
    g.drawString("YR>", width-132, height-67);
    g.setColor(Color.red);
    g.fillRect(width-135,height-60,30,15);
    g.setColor(Color.black);
    g.drawString("<YR", width-132, height-47);

    g.setColor(Color.green);
    g.fillRect(width-100,height-80,30,15);
    g.setColor(Color.black);
    g.drawString("ZR>", width-97, height-67);
    g.setColor(Color.red);
    g.fillRect(width-100,height-60,30,15);
    g.setColor(Color.black);
    g.drawString("<ZR", width-97, height-47);

    g.setColor(Color.green);
    g.fillRect(width-65,height-80,30,15);
    g.setColor(Color.black);
    g.drawString("WR>", width-63, height-67);
    g.setColor(Color.red);
    g.fillRect(width-65,height-60,30,15);
    g.setColor(Color.black);
    g.drawString("<WR", width-63, height-47);

    g.setColor(Color.green);
    g.fillRect(width-30,height-80,30,15);
    g.setColor(Color.black);
    g.drawString("RR>", width-27, height-67);
    g.setColor(Color.red);
    g.fillRect(width-30,height-60,30,15);
    g.setColor(Color.black);
    g.drawString("<RR", width-27, height-47);

    g.setColor(Color.green);
    g.fillRect(width-170,height-120,30,15);
    g.setColor(Color.black);
    g.drawString(" X+", width-167, height-107);
    g.setColor(Color.red);
    g.fillRect(width-170,height-100,30,15);
    g.setColor(Color.black);
    g.drawString(" X-", width-167, height-87);

    g.setColor(Color.green);
    g.fillRect(width-135,height-120,30,15);
    g.setColor(Color.black);
    g.drawString(" Y+", width-132, height-107);
    g.setColor(Color.red);
    g.fillRect(width-135,height-100,30,15);
    g.setColor(Color.black);
    g.drawString(" Y-", width-132, height-87);

    g.setColor(Color.green);
    g.fillRect(width-100,height-120,30,15);
    g.setColor(Color.black);
    g.drawString(" Z+", width-97, height-107);
    g.setColor(Color.red);
    g.fillRect(width-100,height-100,30,15);
    g.setColor(Color.black);
    g.drawString(" Z-", width-97, height-87);

    g.setColor(Color.green);
    g.fillRect(width-65,height-120,30,15);
    g.setColor(Color.black);
    g.drawString(" W+>", width-63, height-107);
    g.setColor(Color.red);
    g.fillRect(width-65,height-100,30,15);
    g.setColor(Color.black);
    g.drawString(" W-", width-63, height-87);


    g.setColor(Color.green);
    g.fillRect(width-30,height-120,30,15);
    g.setColor(Color.black);
    g.drawString("COL", width-27, height-107);
    g.setColor(Color.red);
    g.fillRect(width-30,height-100,30,15);
    g.setColor(Color.black);
    g.drawString("COL", width-27, height-87);



    //  color scale
    for(int i=1; i<=ncol; i++)
    {
      set_color(i, g);
      g.fillRect(20,60+20*i,50,15);
      fstring = Integer.toString(i);
      g.setColor(Color.black);
      g.drawString(lcol[i], 80, 60+12+20*i);
    } // end i

    // min and max
    g.setColor(Color.black);
    g.drawString("reading file "+filein, 20, 60);
    g.drawString("npts=", 95, 75);
    fstring = Integer.toString(DR.num_points);
    g.drawString(fstring, 145, 75);
    g.drawString("npoly=", 95, 90);
    fstring = Integer.toString(DR.num_polys);
    g.drawString(fstring, 145, 90);

    g.drawString("x min=", 190, 60);
    // fstring = Double.toString(ez(xmin));
    fstring = f1.format(xmin);
    g.drawString(fstring, 240, 60);
    g.drawString("x max=", 400, 60);
    //fstring = Double.toString(ez(xmax));
    fstring = f1.format(xmax);
    g.drawString(fstring, 450, 60);
    g.drawString("angx=", 600, 60); // angx
    // fstring = Integer.toString((int)((angx/Pi)*180.0+0.5));
    fstring = intang(angx);
    g.drawString(fstring, 650, 60);
    g.drawString("y min=", 190, 75);
    // fstring = Double.toString(ez(ymin));
    fstring = f1.format(ymin);
    g.drawString(fstring, 240, 75);
    g.drawString("y max=", 400, 75);
    // fstring = Double.toString(ez(ymax));
    fstring = f1.format(ymax);
    g.drawString(fstring, 450, 75);
    g.drawString("angy=", 600, 75); // angy
    // fstring = Integer.toString((int)((angy/Pi)*180.0+0.5));
    fstring = intang(angy);
    g.drawString(fstring, 650, 75);
    g.drawString("z min=", 190, 90);
    //fstring = Double.toString(ez(zmin));
    fstring = f1.format(zmin);
    g.drawString(fstring, 240, 90);
    g.drawString("z max=", 400, 90);
    // fstring = Double.toString(ez(zmax));
    fstring = f1.format(zmax);
    g.drawString(fstring, 450, 90);
    g.drawString("angz=", 600, 90); // angz
    // fstring = Integer.toString((int)((angz/Pi)*180.0+0.5));
    fstring = intang(angz);
    g.drawString(fstring, 650, 90);
    g.drawString("w min=", 190, 105);
    // fstring = Double.toString(ez(wmin));
    fstring = f1.format(wmin);
    g.drawString(fstring, 240, 105);
    g.drawString("w max=", 400,105);
    // fstring = Double.toString(ez(wmax));
    fstring = f1.format(wmax);
    g.drawString(fstring, 450, 105);
    g.drawString("angw=", 600, 105); // angw
    // fstring = Integer.toString((int)((angw/Pi)*180.0+0.5));
    fstring = intang(angw);
    g.drawString(fstring, 650, 105);
    ortho(xmin, ymin, zmin, wmin); // just setup
    if(grid>0)
    {
      // axis labels x1, x2, x3, x5
      g.drawString("x min", 25, height-10);
      g.drawString("x max", (int)(ss*width-65), height-10);
      g.drawString("y min", (int)(ss*width),  (int)(height-25));
      g.drawString("y max", (int)(width-100), (int)(550));
      g.drawString("z min", (int)(width-90),  (int)(535));
      g.drawString("z max", (int)(width-90),  (int)(120));
      g.drawString("w min", (int)(190),       (int)(485));
      g.drawString("w max", (int)(220),       (int)(450));

      // draw box x and y and z and w                  side number
      ortho(xmin, ymin, zmin, wmin); xp1=xp; yp1=yp; // 1
      ortho(xmin, ymax, zmin, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.black);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      ortho(xmin, ymin, zmin, wmin); xp1=xp; yp1=yp; // 2
      ortho(xmax, ymin, zmin, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.black);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymin, zmin, wmin); xp1=xp; yp1=yp; // 3
      ortho(xmax, ymax, zmin, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.black);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      ortho(xmin, ymax, zmin, wmin); xp1=xp; yp1=yp; // 4
      ortho(xmax, ymax, zmin, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.black);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymax, zmin, wmin); xp1=xp; yp1=yp; // 5
      ortho(xmin, ymax, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymax, zmax, wmin); xp1=xp; yp1=yp; // 6
      ortho(xmax, ymax, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymax, zmin, wmin); xp1=xp; yp1=yp; // 7
      ortho(xmax, ymax, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymin, zmin, wmin); xp1=xp; yp1=yp; // 8
      ortho(xmax, ymin, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmax, ymin, zmax, wmin); xp1=xp; yp1=yp; // 9
      ortho(xmax, ymax, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      ortho(xmin, ymin, zmin, wmin); xp1=xp; yp1=yp; // 10
      ortho(xmin, ymin, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymin, zmax, wmin); xp1=xp; yp1=yp; // 11
      ortho(xmin, ymax, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmin, ymin, zmax, wmin); xp1=xp; yp1=yp; // 12
      ortho(xmax, ymin, zmax, wmin); xp2=xp; yp2=yp;
      g.setColor(Color.blue);
      g.drawLine(xp1, yp1, xp2, yp2);
//    }
//    if(grid>0 && (Math.abs(wmax-wmin)>0.01) || wmax>0.01) // option
//    {                                                     // w box
      // draw box x4                                  side number
      g.setColor(Color.red);
      ortho(xmin, ymin, zmin, wmax); xp1=xp; yp1=yp; // 13
      ortho(xmin, ymax, zmin, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymin, zmin, wmax); xp1=xp; yp1=yp; // 14
      ortho(xmax, ymin, zmin, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymin, zmin, wmax); xp1=xp; yp1=yp; // 15
      ortho(xmax, ymax, zmin, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymax, zmin, wmax); xp1=xp; yp1=yp; // 16
      ortho(xmax, ymax, zmin, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymax, zmin, wmax); xp1=xp; yp1=yp; // 17
      ortho(xmin, ymax, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymax, zmax, wmax); xp1=xp; yp1=yp; // 18
      ortho(xmax, ymax, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymax, zmin, wmax); xp1=xp; yp1=yp; // 19
      ortho(xmax, ymax, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymin, zmin, wmax); xp1=xp; yp1=yp; // 20
      ortho(xmax, ymin, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmax, ymin, zmax, wmax); xp1=xp; yp1=yp; // 21
      ortho(xmax, ymax, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymin, zmin, wmax); xp1=xp; yp1=yp; // 22
      ortho(xmin, ymin, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymin, zmax, wmax); xp1=xp; yp1=yp; // 23
      ortho(xmin, ymax, zmax, wmax); xp2=xp; yp2=yp;
      g.fillOval(xp2-4, yp2-4, 8, 8);
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
      ortho(xmin, ymin, zmax, wmax); xp1=xp; yp1=yp; // 24
      ortho(xmax, ymin, zmax, wmax); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      g.fillOval(xp1-4, yp1-4, 8, 8);
      g.fillOval(xp2-4, yp2-4, 8, 8);
    }

    // plot function or data with optional colors, points

    k = 0; // polygons auto incr
    avgx = 0.0; // for sort
    avgy = 0.0;
    avgz = 0.0;
    avgw = 0.0;
    for(int i=0; i<DR.num_polys; i++)
    {
      kpoly = k;
      pts = DR.data_polys[k++];
      if(debug) System.out.println("k="+(k-1)+", ipoly="+i+", pts="+pts);
      npp = 0;
      for(int j=0; j<pts-1; j++)
      {
	if(j<1)
	{
          pt = DR.data_polys[k++];
          x1 = xyzw[pt-1][0]; // DR.data_points[pt-1].x;
          y1 = xyzw[pt-1][1]; // DR.data_points[pt-1].y;
          z1 = xyzw[pt-1][2]; // DR.data_points[pt-1].z;
          w1 = xyzw[pt-1][3]; // DR.data_points[pt-1].w;
	  if(debug) System.out.println("pt="+(pt-1));
          ortho(x1, y1, z1, w1); xp1=xp; yp1=yp;
	  xvec[npp] = xp;
	  yvec[npp] = yp;
	  npp++;
          makeP(x1, y1, z1, w1, P1);
	  doavg(x1, y1, z1, w1);
          g.setColor(Color.black);
	  g.drawString(Integer.toString(pt),xp1+6,yp1+2);
          x0 = x1;
          y0 = y1;
          z0 = z1;
	  w0 = w1;
	} // end start point
        pt = DR.data_polys[k++];
        x2 = xyzw[pt-1][0]; // DR.data_points[pt-1].x;
        y2 = xyzw[pt-1][1]; // DR.data_points[pt-1].y;
        z2 = xyzw[pt-1][2]; // DR.data_points[pt-1].z;
        w2 = xyzw[pt-1][3]; // DR.data_points[pt-1].w;
	if(debug) System.out.println("pt="+(pt-1));
        ortho(x2, y2, z2, w2); xp2=xp; yp2=yp;
	xvec[npp] = xp;
	yvec[npp] = yp;
	npp++;
        if(npp==2)  makeP(x2, y2, z2, w2, P2); 
        if(npp==3)  makeP(x2, y2, z2, w2, P3); 
	doavg(x2, y2, z2, w2);
	g.drawString(Integer.toString(pt),xp2+6,yp2+2);
        g.setColor(Color.black);
        g.drawLine(xp1, yp1, xp2, yp2);
        if(points>0) g.fillOval(xp1-4, yp1-4, 8, 8);
        if(points>0) g.fillOval(xp2-4, yp2-4, 8, 8);
        xp1 = xp2;
        yp1 = yp2;
      } // end this poly, now close
      if(pts>2)
      {
        ortho(x0, y0, z0, w0); xp2=xp; yp2=yp;
        g.setColor(Color.black);
        g.drawLine(xp1, yp1, xp2, yp2);
      }
      avgx = avgx/(double)pts;
      prix[i] = avgx;
      pripoly[i] = kpoly;
      avgy = avgy/(double)pts;
      priy[i] = avgy;
      avgz = avgz/(double)pts;
      priz[i] = avgz;
      avgw = avgw/(double)pts;
      priw[i] = avgw;
    } // end polygons
      
    if(fillc) // fill colored polygon
    {
      if(xsort)
      {
        System.out.println("X color sort");
        sort(prix, num_polys, priindex);
      }
      if(ysort)
      {
        System.out.println("Y color sort");
        sort(priy, num_polys, priindex);
      }
      if(zsort)
      {
        System.out.println("Z color sort");
        sort(priz, num_polys, priindex);
      }
      if(wsort)
      {
        System.out.println("W color sort");
        sort(priw, num_polys, priindex);
      }
      k = 0;
      for(int i=0; i<num_polys; i++)
      {
	// System.out.println("poly="+priindex[i]+", k="+pripoly[priindex[i]]+
        //                    ", x="+prix[i]);
	k = pripoly[priindex[i]];
        pts = DR.data_polys[k++];
        if(debug) System.out.println("k="+(k-1)+", ipoly="+i+", pts="+pts);
        npp = 0;
        for(int j=0; j<pts-1; j++)
        {
	  if(j<1)
	  {
            pt = DR.data_polys[k++];
            x1 = xyzw[pt-1][0]; // DR.data_points[pt-1].x;
            y1 = xyzw[pt-1][1]; // DR.data_points[pt-1].y;
            z1 = xyzw[pt-1][2]; // DR.data_points[pt-1].z;
            w1 = xyzw[pt-1][3]; // DR.data_points[pt-1].w;
	    if(debug) System.out.println("pt="+(pt-1));
            ortho(x1, y1, z1, w1); xp1=xp; yp1=yp;
	    xvec[npp] = xp;
	    yvec[npp] = yp;
	    npp++;
            makeP(x1, y1, z1, w1, P1);
	  } // end start point
          pt = DR.data_polys[k++];
          x2 = xyzw[pt-1][0]; // DR.data_points[pt-1].x;
          y2 = xyzw[pt-1][1]; // DR.data_points[pt-1].y;
          z2 = xyzw[pt-1][2]; // DR.data_points[pt-1].z;
          w2 = xyzw[pt-1][3]; // DR.data_points[pt-1].w;
          ortho(x2, y2, z2, w2); xp2=xp; yp2=yp;
	  xvec[npp] = xp;
	  yvec[npp] = yp;
	  npp++;
          if(npp==2)  makeP(x2, y2, z2, w2, P2); 
          if(npp==3)  makeP(x2, y2, z2, w2, P3); 
        }
        normal4(P1, P2, P3);
	if(xsort)
	{
          red = (float)((ny+1.0)/2.0); // change color based on normal vector
          grn = (float)((nz+1.0)/2.0); // normal from -1 to +1, color 0 to 1
          blu = (float)((nw+1.0)/2.0); // must be at least 2
	}
	if(ysort)
	{
          red = (float)((nx+1.0)/2.0); // change color based on normal vector
          grn = (float)((nz+1.0)/2.0); // normal from -1 to +1, color 0 to 1
          blu = (float)((nw+1.0)/2.0); // must be at least 2
	}
	if(zsort)
	{
          red = (float)((nx+1.0)/2.0); // change color based on normal vector
          grn = (float)((ny+1.0)/2.0); // normal from -1 to +1, color 0 to 1
          blu = (float)((nw+1.0)/2.0); // must be at least 2
	}
	if(wsort)
	{
          red = (float)((nx+1.0)/2.0); // change color based on normal vector
          grn = (float)((ny+1.0)/2.0); // normal from -1 to +1, color 0 to 1
          blu = (float)((nz+1.0)/2.0); // must be at least 2
	}
        g.setColor(new Color(red,grn,blu));
        g.fillPolygon(xvec, yvec, npp);
      }
    } //end fill colored polygon
  } // end paint

  void makeP(double x, double y, double z, double w, double p[])
  {
    p[0] = x;
    p[1] = y;
    p[2] = z;
    p[3] = w;
  } // end makeP

  void doavg(double x, double y, double z, double w)
  {
    avgx += x;
    avgy += y;
    avgz += z;
    avgw += w;
  } // end doavg  later /num_poly
    
  // orthographic projection, sets xp,yp
  void ortho(double x, double y, double z, double w)
  {
      //    x   y   z   w
    double x1, x2, x3, x4;  // dimension space
    double xc, yc; // coordinates centered for x4 scaling
    double x4fx, x4fy, x4fz; // scaling x4 smaller
    
    if(debug) System.out.println("ortho x="+x+" y="+y+" z="+z+" w="+w);
    ss = 0.6666; //  1,2,3 dimension scale
    za = 0.3333; //  z makes angle adds to x and y 45 degree here 
    ws = ss/2.0; //  4th dimension scale

    if(Math.abs(xmax-xmin)<0.001) { x1 = 0.0; }
    else { x1 = (x-xmin)/(xmax-xmin); }  // 0.0 to 1.0
    if(Math.abs(ymax-ymin)<0.001) { x2 = 0.0; }
    else { x2 = (y-ymin)/(ymax-ymin); }
    if(Math.abs(zmax-zmin)<0.001) { x3 = 0.0; }
    else { x3 = (z-zmin)/(zmax-zmin); }
    if(Math.abs(wmax-wmin)<0.001) { x4 = 0.0; }
    else { x4 = (w-wmin)/(wmax-wmin); } 
    if(debug) System.out.println("x="+x1+" y="+x2+" z="+x3+" w="+x4);
    xc = 0.5;  // center,  everything 0 to 1
    yc = 0.5;

    x4fx = 1.0-x4*ss;     // +0.1;
    x4fy = 1.0-x4*0.60;   // 0.7071; // -0.1;
    x4fz = 1.0-x4*za;
    if(debug) System.out.println("x4fx="+x4fx+" x4fy="+x4fy);
	
    x1 = x1*ss; // scale
    x2 = x2*ss; //         
    x3 = x3*ss; //
        
    x1 = x1+x3*za*x4fz; // ortho y is up, z is back, w is offset and small
    x2 = x2+x3*za*x4fz;

    x1 = x1-xc; // center
    x2 = x2-yc;
    // System.out.println("center x1="+x1+" x2="+x2);
    
    x1 = x1*x4fx; // scale x4
    x2 = x2*x4fy;
    // System.out.println("x4 scale x1="+x1+" x2="+x2);

    x1 = x1 + xc; // uncenter
    x2 = x2 + yc;
    if(debug) System.out.println("0-1 pos x1="+x1+" x2="+x2);
    
    x1 = x1*(double)width*0.9; // plot scale
    x2 = x2*(double)height*0.9;

    x1 = x1+25.0; // plot offset
    x2 = (double)height-x2-25.0; // away from boarder and text
    xp = (int)x1; // output global
    yp = (int)x2;
  } // end ortho

  void set_color(int i, Graphics g)
  {
    if(i<0) {g.setColor(Color.black); return;}
    if(i>4) {g.setColor(Color.yellow); return; }
    switch(i)
    {
      case 0: g.setColor(Color.black); break;
      case 1: g.setColor(Color.black); break;  // x1 > 0.01  x
      case 2: g.setColor(Color.black); break;  // x2 > 0.01  y
      case 3: g.setColor(Color.blue); break;   // x3 > 0.01  z
      case 4: g.setColor(Color.red); break;    // x4 > 0.01  w
    }
  } // end set_color

  // matrix utility funcions
    
  void prtpt(String name, double A[])
  {
    int len1 = A.length;
    System.out.print(name+"[]");
    for(int i=0; i<len1; i++)
    {
      System.out.print(" "+A[i]);
    }
    System.out.println(" ");
  } // end prtpt

  void prtmat(String name, double A[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      System.out.print(name+"["+i+"]");
      for(int j=0; j<len2; j++)
      {
	System.out.print(" "+A[i][j]);
      }
      System.out.println(" ");
    }
    System.out.println(" ");
  } // end prtmat

  void matmul(double A[][], double B[][], double C[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    if(len1 != len2)
    {
      System.out.println("bad for here, non square matrix");
      return;
    }
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len1; j++)
      {
	C[i][j] = 0.0;
	for(int k=0; k<len1; k++)
	{
	  C[i][j] += A[i][k]*B[k][j];
	} // end k
      } // end j
    } // end i
  } // end matmul

  void matadd(double A[][], double B[][], double C[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len2; j++)
      {
	C[i][j] += A[i][j]+B[i][j];
      } // end j
    } // end i
  } // end matadd

  void mattrans(double A[][], double C[][]) // simple transpose
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len2; j++)
      {
	C[i][j] = A[j][i];
      } // end j
    } // end i
  } // end mattrans
    
  void matvec(double A[][], double V[], double W[])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    if(len1 != len2)
    {
      System.out.println("bad for here, non square matrix");
      return;
    }
    for(int i=0; i<len1; i++)
    {
      W[i] = 0.0;
      for(int j=0; j<len1; j++)
      {
        W[i] += A[i][j]*V[j];
      } // end j
    } // end i
  } // end matvec

  void matzero(double A[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len2; j++)
      {
	A[i][j] = 0.0;
      } // end j
    } // end i
  } // end matzero

  void matident(double A[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len2; j++)
      {
	A[i][j] = 0.0;
	if(i == j) A[i][j] = 1.0;
      } // end j
    } // end i
  } // end matident

  void matcopy(double A[][], double C[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    for(int i=0; i<len1; i++)
    {
      for(int j=0; j<len2; j++)
      {
	C[i][j] = A[i][j];
      } // end j
    } // end i
  } // end matcopy

  void veccopy(double A[], double C[])
  {
    int len1 = A.length;
    for(int i=0; i<len1; i++)
    {
      C[i] = A[i];
    } // end i
  } // end veccopy

  double matdet(double A[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    double B[][] = new double[len1][len2];
    matcopy(A, B);
    return Matrix.determinant(B);  
  } // end matdet

  void matinv(double A[][], double C[][])
  {
    int len1 = A.length;
    int len2 = A[0].length;
    matcopy(A, C);
    Matrix.invert(C);  
  } // end matinv

  // data manipulation functions, create matrix for matvec on point
    
  void makePM(double xmin, double xmax, double ymin, double ymax,
	      double zmin, double zmax, double wmin, double wmax,
	      double near, double far, double PM[][])
  {
    PM[0][0] = 2.0*near/(xmax-xmin);
    PM[0][1] = 0.0;
    PM[0][2] = 0.0;
    PM[0][3] = (xmax+xmin)/(xmax-xmin);
    PM[0][4] = 0.0;
    PM[1][0] = 0.0;
    PM[1][1] = 2.0*near/(ymax-ymin);
    PM[1][2] = 0.0;
    PM[1][3] = (ymax+ymin)/(ymax-ymin);
    PM[1][4] = 0.0;
    PM[2][0] = 0.0;
    PM[2][1] = 0.0;
    PM[2][2] = 2.0*near/(zmax-zmin);
    PM[2][3] = (zmax+zmin)/(zmax-zmin);
    PM[2][4] = 0.0;
    PM[3][0] = 0.0;
    PM[3][1] = 0.0;
    PM[3][2] = 0.0;
    PM[3][3] = -(far+near)/(far-near);
    PM[3][4] = -2.0*far*near*(far-near); // ?? (far+near)*
    PM[4][0] = 0.0;
    PM[4][1] = 0.0;
    PM[4][2] = 0.0;
    PM[4][3] = -1.0;
    PM[4][4] = 0.0;
  } // end makePM

  void makeTM(double fx, double fy, double fz, double fw, double TM[][])
  {
    TM[0][0] = 1.0;
    TM[0][1] = 0.0;
    TM[0][2] = 0.0;
    TM[0][3] = 0.0;
    TM[0][4] = fx;
    TM[1][0] = 0.0;
    TM[1][1] = 1.0;
    TM[1][2] = 0.0;
    TM[1][3] = 0.0;
    TM[1][4] = fy;
    TM[2][0] = 0.0;
    TM[2][1] = 0.0;
    TM[2][2] = 1.0;
    TM[2][3] = 0.0;
    TM[2][4] = fz;
    TM[3][0] = 0.0;
    TM[3][1] = 0.0;
    TM[3][2] = 0.0;
    TM[3][3] = 1.0;
    TM[3][4] = fw;
    TM[4][0] = 0.0;
    TM[4][1] = 0.0;
    TM[4][2] = 0.0;
    TM[4][3] = 0.0;
    TM[4][4] = 1.0;
  } // end makeTM

  void makeSM(double sx, double sy, double sz, double sw, double SM[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	SM[i][j] = 0.0;
      }
    }
    SM[0][0] = sx;
    SM[1][1] = sy;
    SM[2][2] = sz;
    SM[3][3] = sw;
    SM[4][4] = 1.0;
  } // end makeSM

  void makeXY(double ax, double XY[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	XY[i][j] = 0.0;
      }
    }
    XY[0][0] = 1.0;
    XY[1][1] = 1.0;
    XY[2][2] = Math.cos(ax);
    XY[2][3] = -Math.sin(ax);
    XY[3][2] = Math.sin(ax);
    XY[3][3] = Math.cos(ax);
    XY[4][4] = 1.0;
  } // end makeXY

  void makeXZ(double ax, double XZ[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	XZ[i][j] = 0.0;
      }
    }
    XZ[0][0] = 1.0;
    XZ[1][1] = Math.cos(ax);
    XZ[1][3] = -Math.sin(ax);
    XZ[2][2] = 1.0;
    XZ[3][1] = Math.sin(ax);
    XZ[3][3] = Math.cos(ax);
    XZ[4][4] = 1.0;
  } // end makeXZ

  void makeXW(double ax, double XW[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	XW[i][j] = 0.0;
      }
    }
    XW[0][0] = 1.0;
    XW[1][1] = Math.cos(ax);
    XW[1][2] = -Math.sin(ax);
    XW[2][1] = Math.sin(ax);
    XW[2][2] = Math.cos(ax);
    XW[3][3] = 1.0;
    XW[4][4] = 1.0;
  } // end makeXW

  void makeYZ(double ay, double YZ[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	YZ[i][j] = 0.0;
      }
    }
    YZ[0][0] = Math.cos(ay);
    YZ[0][3] = -Math.sin(ay);
    YZ[1][1] = 1.0;
    YZ[2][2] = 1.0;
    YZ[3][0] = Math.sin(ay);
    YZ[3][3] = Math.cos(ay);
    YZ[4][4] = 1.0;
  } // end makeYZ

  void makeYW(double ay, double YW[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	YW[i][j] = 0.0;
      }
    }
    YW[0][0] = Math.cos(ay);
    YW[0][2] = -Math.sin(ay);
    YW[1][1] = 1.0;
    YW[2][0] = Math.sin(ay);
    YW[2][2] = Math.cos(ay);
    YW[3][3] = 1.0;
    YW[4][4] = 1.0;
  } // end makeYW

  void makeZW(double az, double ZW[][])
  {
    for(int i=0; i<5; i++)
    {
      for(int j=0; j<5; j++)
      {
	ZW[i][j] = 0.0;
      }
    }
    ZW[0][0] = Math.cos(az);
    ZW[0][1] = -Math.sin(az);
    ZW[1][0] = Math.sin(az);
    ZW[1][1] = Math.cos(az);
    ZW[2][2] = 1.0;
    ZW[3][3] = 1.0;
    ZW[4][4] = 1.0;
  } // end makeZW

  void makeXR(double ax, double XR[][])
  {
    double XY[][] = new double[5][5];
    double XZ[][] = new double[5][5];
    double XW[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    if(debug2) System.out.println("makeXR  ax="+ax);
    if(ax>=0.0)
    {
      makeXY(ax,XY);
      if(debug2) prtmat("XY",XY);
      makeXZ(ax,XZ);
      if(debug2) prtmat("XZ",XZ);
      makeXW(ax,XW);
      if(debug2) prtmat("XW",XW);
      matmul(XY, XZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, XW, XR);
    }
    else
    {
      makeXY(-ax,XY);
      if(debug2) prtmat("XY",XY);
      makeXZ(-ax,XZ);
      if(debug2) prtmat("XZ",XZ);
      makeXW(-ax,XW);
      if(debug2) prtmat("XW",XW);
      matmul(XY, XZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, XW, M1);
      mattrans(M1, XR);
    }
    if(debug2) prtmat("XR",XR);
    if(debug2) System.out.println(" ");
  } // end makeXR

  void makeXRT(double ax, double XRT[][]) // using transpose
  {
    double XY[][] = new double[5][5];
    double XZ[][] = new double[5][5];
    double XZT[][] = new double[5][5];
    double XW[][] = new double[5][5];
    double XWT[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    if(debug2) System.out.println("makeXRT ax="+ax+"  positive");
    makeXY(ax,XY);
    if(debug2) prtmat("XY",XY);
    makeXZ(ax,XZ);
    if(debug2) prtmat("XZ",XZ);
    matmul(XY, XZ, M);
    if(debug2) prtmat("M",M);
    makeXW(ax,XW);
    if(debug2) prtmat("XW",XW);
    matmul(M, XW, M1);
    mattrans(M1, XRT);
    if(debug2) prtmat("XRT",XRT);
    if(debug2) System.out.println(" ");
  } // end makeXRT

  void makeYR(double ay, double YR[][])
  {
    double XY[][] = new double[5][5];
    double YZ[][] = new double[5][5];
    double YW[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    if(debug2) System.out.println("makeYR  ay="+ay);
    if(ay>=0.0)
    {
      makeXY(ay,XY);
      if(debug2) prtmat("XY",XY);
      makeYZ(ay,YZ);
      if(debug2) prtmat("YZ",YZ);
      makeYW(ay,YW);
      if(debug2) prtmat("YW",YW);
      matmul(XY, YZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, YW, YR);
    }
    else
    {
      makeXY(-ay,XY);
      if(debug2) prtmat("XY",XY);
      makeYZ(-ay,YZ);
      if(debug2) prtmat("YZ",YZ);
      makeYW(-ay,YW);
      if(debug2) prtmat("YW",YW);
      matmul(XY, YZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, YW, M1);
      mattrans(M1, YR);
    }
    if(debug2) prtmat("YR",YR);
    if(debug2) System.out.println(" ");
  } // end makeYR

  void makeZR(double az, double ZR[][])
  {
    double XZ[][] = new double[5][5];
    double YZ[][] = new double[5][5];
    double ZW[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    if(debug2) System.out.println("makeZR az="+az);
    if(az>=0.0)
    {
      makeXZ(az,XZ);
      if(debug2) prtmat("XZ",XZ);
      makeYZ(az,YZ);
      if(debug2) prtmat("YZ",YZ);
      makeZW(az,ZW);
      if(debug2) prtmat("ZW",ZW);
      matmul(XZ, YZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, ZW, ZR);
    }
    else
    {
      makeXZ(-az,XZ);
      if(debug2) prtmat("XZ",XZ);
      makeYZ(-az,YZ);
      if(debug2) prtmat("YZ",YZ);
      makeZW(-az,ZW);
      if(debug2) prtmat("ZW",ZW);
      matmul(XZ, YZ, M);
      if(debug2) prtmat("M",M);
      matmul(M, ZW, M1);
      mattrans(M1, ZR);
    }
    if(debug2) prtmat("ZR",ZR);
    if(debug2) System.out.println(" ");
  } // end makeZR

  void makeWR(double aw, double WR[][])
  {
    double XW[][] = new double[5][5];
    double YW[][] = new double[5][5];
    double ZW[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    if(debug2) System.out.println("makeWR aw="+aw);
    if(aw>=0.0)
    {
      makeXW(aw,XW);
      if(debug2) prtmat("XW",XW);
      makeYW(aw,YW);
      if(debug2) prtmat("YW",YW);
      makeZW(aw,ZW);
      if(debug2) prtmat("ZW",ZW);
      matmul(XW, YW, M);
      if(debug2) prtmat("M",M);
      matmul(M, ZW, WR);
    }
    else
    {
      makeXW(-aw,XW);
      if(debug2) prtmat("XW",XW);
      makeYW(-aw,YW);
      if(debug2) prtmat("YW",YW);
      makeZW(-aw,ZW);
      if(debug2) prtmat("ZW",ZW);
      matmul(XW, YW, M);
      if(debug2) prtmat("M",M);
      matmul(M, ZW, M1);
      mattrans(M1, WR);
    }
    if(debug2) prtmat("WR",WR);
    if(debug2) System.out.println(" ");
  } // end makeWR

  void makeXYZWR(double ax, double ay, double az, double aw, double XYZWR[][])
  {
    double XR[][] = new double[5][5];
    double YR[][] = new double[5][5];
    double ZR[][] = new double[5][5];
    double WR[][] = new double[5][5];
    double M[][]  = new double[5][5];
    double M1[][]  = new double[5][5];
    double M2[][]  = new double[5][5];
    if(debug2) System.out.println("makeXYZWR ax="+ax+", ay="+ay+
				  ", az="+az+", aw="+aw);
    if(aw>=0.0)
    {
      makeXR(ax,XR);
      if(debug2) prtmat("XR",XR);
      makeYR(ay,YR);
      if(debug2) prtmat("YR",YR);
      makeZR(az,ZR);
      if(debug2) prtmat("ZR",ZR);
      makeWR(aw,WR);
      if(debug2) prtmat("WR",WR);
      matmul(XR, YR, M);
      if(debug2) prtmat("M",M);
      matmul(ZR, M, M1);
      if(debug2) prtmat("M1",M1);
      matmul(WR, M1, XYZWR);
    }
    else
    {
      makeXR(-ax,XR);
      if(debug2) prtmat("XR",XR);
      makeYR(-ay,YR);
      if(debug2) prtmat("YR",YR);
      makeZR(-az,ZR);
      if(debug2) prtmat("ZR",ZR);
      makeWR(-aw,WR);
      if(debug2) prtmat("WR",WR);
      matmul(XR, YR, M);
      if(debug2) prtmat("M",M);
      matmul(ZR, M, M1);
      if(debug2) prtmat("M1",M1);
      matmul(WR, M1, M2);
      mattrans(M2, XYZWR);
    }
    if(debug2) prtmat("XYZWR",XYZWR);
    if(debug2) System.out.println(" ");
  } // end makeXYZWR

  void getpr(int i, double pr[])
  {
    pr[0] = DR.data_points[i].x;
    pr[1] = DR.data_points[i].y;
    pr[2] = DR.data_points[i].z;
    pr[3] = DR.data_points[i].w;
    pr[4] = 1.0; // for rotation
  } // end getpr

  String intang(double ang)
  {
    double deg = (ang/Pi)*180.0; // radians to degrees
    if(deg>0.0) deg = deg + 0.5; // round
    if(deg<0.0) deg = deg - 0.5;
    return Integer.toString((int)deg);
  } // end intang

  void maxmin(double M[][])
  {
    xmin = M[0][0];
    xmax = M[0][0];
    ymin = M[0][1];
    ymax = M[0][1];
    zmin = M[0][2];
    zmax = M[0][2];
    wmin = M[0][3];
    wmax = M[0][3];
    for(int i=1; i<num_points; i++) // vertices
    {
      xmin = Math.min(xmin, M[i][0]);
      xmax = Math.max(xmax, M[i][0]);
      ymin = Math.min(ymin, M[i][1]);
      ymax = Math.max(ymax, M[i][1]);
      zmin = Math.min(zmin, M[i][2]);
      zmax = Math.max(zmax, M[i][2]);
      wmin = Math.min(wmin, M[i][3]);
      wmax = Math.max(wmax, M[i][3]);
    }
  } // end maxmin

  double ez(double a)
  {
    if(Math.abs(a)<1.0e-6) return 0.0;
    if(Math.abs(a-1.0)<1.0e-6) return 1.0;
    if(Math.abs(a+1.0)<1.0e-6) return -1.0;
    if(Math.abs(a-2.0)<1.0e-6) return 2.0;
    if(Math.abs(a+2.0)<1.0e-6) return -2.0;
    return a;
  } // end ez

  // from  normal44c.java   normal4(p1, p2, p3);
  //    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nw="+nw);
  void normal4(double a[], double b[], double c[])  // 3 point 4D triangle
  {
    // compute vector of side a,b  vector of side a,c  vector of side b,c
    // cross product of these 3 vectors is vector of normal to triangle
    double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, ss;
    
    a1 = a[0]; // ax
    a2 = a[1]; // ay
    a3 = a[2]; // az
    a4 = a[3]; // aw
    b1 = b[0]; // bx
    b2 = b[1]; // by
    b3 = b[2]; // bz
    b4 = b[3]; // bw
    c1 = c[0]; // cx
    c2 = c[1]; // cy
    c3 = c[2]; // cz
    c4 = c[3]; // cw
    //       |a2 a3 a4|         |a1 a3 a4|        |a1 a2 a4|         |a1 a2 a3|
    // nx=det|b2 b3 b4|  ny=-det|b1 b3 b4|  nz=det|b1 b2 b4|  nw=-det|b1 b2 b3|
    //       |c2 c3 c4|         |c1 c3 c4|        |c1 c2 c4|         |c1 c2 c3|

    nx = a2*b3*c4 + a3*a4*c2 + a4*c3*b2 -
	 c2*b3*a4 - c3*b4*a2 - c4*a3*b2;
    ny = a1*b3*c4 + a3*a4*c1 + a4*c3*b1 -
	 c1*b3*a4 - c3*b4*a1 - c4*a3*b1;
    ny = -ny;
    nz = a1*b2*c4 + a2*b4*c2 + a4*c2*b1 -
	 c1*b2*a4 - c2*b4*a1 - c4*a2*b1;
    nw = a1*b2*c3 + a2*b3*c1 + c3*b1*a2 -
	 c1*b2*a3 - c2*b3*a1 - c3*a2*b1;
    nw = -nw;
	 
    ss = Math.sqrt(nx*nx+ny*ny+nz*nz+nw*nw);
    if(Math.abs(ss)>0.01)
    {
      nx = nx / ss;
      ny = ny / ss;
      nz = nz / ss;
      nw = nw / ss;
    }
    // r=(nx+1)/2 g=(ny+2)/2 b=(nz+1)/2 add (nw+1)/2?  g.setcolor(r,g,b)
    // g.fillPolygon(xvec, yvec, np)
  } // end normal4

  void test_matrix()
  {

    // generate and test many matrix
    matzero(M);
    prtmat("matzero",M);
    matident(M);
    prtmat("matident",M);
    matcopy(M,M1);
    prtmat("M1=matcopyM",M1);
    System.out.println(" ");
    makeXY(ax,XY);
    prtmat("XY",XY);
    d = matdet(XY);
    System.out.println("matdet(XY)="+d);
    mattrans(XY,M);
    prtmat("XY^t",M);
    matinv(XY,M1);
    prtmat("XY^-1",M1);
    System.out.println("check XY transpose = XY inverse");
    System.out.println(" ");

	   
    near = (wmax-wmin); // may be 2*  or 5*
    far = 2.0*near; // may be 5 or even 10 or 20
    makePM(xmin, xmax, ymin, ymax, zmin, zmax, wmin, wmax, near, far, PM);
    prtmat("PM",PM);
    matvec(PM,pxi,p);
    prtpt("PM*pxi",p);
    matvec(PM,pyi,p);
    prtpt("PM*pyi",p);
    matvec(PM,pzi,p);
    prtpt("PM*pzi",p);
    matvec(PM,pwi,p);
    prtpt("PM*pwi",p);
    matvec(PM,pnxi,p);
    prtpt("PM*pnxi",p);
    matvec(PM,pnyi,p);
    prtpt("PM*pnyi",p);
    matvec(PM,pnzi,p);
    prtpt("PM*pnzi",p);
    matvec(PM,pnwi,p);
    prtpt("PM*pnwi",p);
    matvec(PM,p0000,p);
    prtpt("PM*p0000",p);
    matvec(PM,p1111,p);
    prtpt("PM*p1111",p);
    matvec(PM,p1234,p);
    prtpt("PM*p1234",p);
    System.out.println(" ");

    makeTM(fx, fy, fz, fw, TM); // translation matrix
    prtmat("TM",TM);
    prtpt("p1234", p1234);
    matvec(TM, p1234, p);
    prtpt("TM*p1234", p);
    System.out.println("reverse translate");
    makeTM(-fx, -fy, -fz, -fw, TM);
    prtmat("TM",TM);
    prtpt("p", p);
    matvec(TM, p, p1);
    prtpt("TM*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");
    
    makeSM(sx, sy, sz, sw, SM); // scale matrix
    prtmat("SM",SM);
    prtpt("p1234", p1234);
    matvec(SM, p1234, p);
    prtpt("SM*p1234", p);
    System.out.println("reverse scale");
    makeSM(1.0/sx, 1.0/sy, 1.0/sz, 1.0/sw, SM);
    prtmat("SM",SM);
    prtpt("p", p);
    matvec(SM, p, p1);
    prtpt("SM*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    makeXY(ax,XY); // rotate about X,Y axis
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    prtmat("XY",XY);
    prtpt("p1234", p1234);
    matvec(XY, p1234, p);
    prtpt("XY*p1234", p);
    System.out.println("reverse rotate XY");
    makeXY(-ax,XY);
    prtmat("XY",XY);
    prtpt("p", p);
    matvec(XY, p, p1);
    prtpt("XY*p", p1);
    System.out.println("check == p1234");
    makeXY(ax,XY);
    matvec(XY, p1234, p);
    prtpt("15 XY*p1234", p);
    System.out.println("rotate by XY again");
    matvec(XY, p, p1);
    prtpt("30 XY*p1", p1);
    makeXY(2.0*ax,XY);
    matvec(XY, p1234, p2);
    prtpt("30 XY*p1234", p2);
    System.out.println("check p1== p2");
    System.out.println(" ");
    System.out.println(" ");

    makeXZ(ax,XZ);  // rotate about X,Z axis
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    prtmat("XZ",XZ);
    prtpt("p1234", p1234);
    matvec(XZ, p1234, p);
    prtpt("XZ*p1234", p);
    System.out.println("reverse rotate XZ");
    makeXZ(-ax,XZ);
    prtmat("XZ",XZ);
    prtpt("p", p);
    matvec(XZ, p, p1);
    prtpt("XZ*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    makeXW(ax,XW);  // rotate about X,W axis
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    prtmat("XW",XW);
    prtpt("p1234", p1234);
    matvec(XW, p1234, p);
    prtpt("XW*p1234", p);
    System.out.println("reverse rotate XW");
    makeXW(-ax,XW);
    prtmat("XW",XW);
    prtpt("p", p);
    matvec(XW, p, p1);
    prtpt("XW*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    makeYZ(ay,YZ); // rotate about Y,Z axis
    System.out.println("ay="+ay+" rad,  ay="+((ay/Pi)*180.0)+" deg");  
    prtmat("YZ",YZ);
    prtpt("p1234", p1234);
    matvec(YZ, p1234, p);
    prtpt("YZ*p1234", p);
    System.out.println("reverse rotate YZ");
    makeYZ(-ay,YZ);
    prtmat("YZ",YZ);
    prtpt("p", p);
    matvec(YZ, p, p1);
    prtpt("YZ*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    makeYW(ay,YW);  // rotate about Y,W axis
    System.out.println("ay="+ay+" rad,  ay="+((ay/Pi)*180.0)+" deg");  
    prtmat("YW",YW);
    prtpt("p1234", p1234);
    matvec(YW, p1234, p);
    prtpt("YW*p1234", p);
    System.out.println("reverse rotate YW");
    makeYW(-ay,YW);
    prtmat("YW",YW);
    prtpt("p", p);
    matvec(YW, p, p1);
    prtpt("YW*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    makeZW(az,ZW); // rotate about Z,W axis
    System.out.println("az="+az+" rad,  az="+((az/Pi)*180.0)+" deg");  
    prtmat("ZW", ZW);
    prtpt("p1234", p1234);
    matvec(ZW, p1234, p);
    prtpt("ZW*p1234", p);
    System.out.println("reverse rotate ZW");
    makeZW(-az,ZW);
    prtmat("ZW", ZW);
    prtpt("p", p);
    matvec(ZW, p, p1);
    prtpt("ZW*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("rotate around X axis 3 matvec");
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    makeXY(ax,XY);
    prtmat("XY",XY);
    prtpt("p1234", p1234);
    matvec(XY, p1234, p);
    prtpt("XY*p1234", p);
    makeXZ(ax,XZ);
    prtmat("XZ",XZ);
    prtpt("p", p);
    matvec(XZ, p, p1);
    prtpt("XZ*p", p1);
    makeXW(ax,XW);
    prtmat("XW",XW);
    prtpt("p1", p1);
    matvec(XW, p1, p2);
    prtpt("p2 = XW*p1", p2);
    System.out.println("reverse rotate 3 matvec, -ax wont work");
    makeXW(-ax,XW);                // inverse = trnaspose easier
    prtmat("XW",XW);
    prtpt("p2", p2);
    matvec(XW, p2, p1);
    prtpt("p1 = XW*p2", p1);
    makeXZ(-ax,XZ);
    prtmat("XZ",XZ);
    prtpt("p1", p1);
    matvec(XZ, p1, p);
    prtpt("p = XZ*p1", p);
    makeXY(-ax,XY);
    prtmat("XY",XY);
    prtpt("p", p);
    matvec(XY, p, p3);
    prtpt("p3 = XY*p", p3);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");


    System.out.println("rotate around X axis using XR, transpose");    
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    makeXR(ax, XR);
    prtmat("XR", XR);
    matvec(XR, p1234, p);
    prtpt("XR*p1234", p);
    System.out.println("reverse rotate using XRT");
    makeXRT(ax, XRT); // same angle, transpose matrix
    prtmat("XRT", XRT);
    matvec(XRT, p, p1);
    prtpt("XRT*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("check makeXRT == mattrans(XR)");
    makeXRT(ax, XRT);
    prtmat("XRT", XRT);
    makeXR(ax, XR);
    mattrans(XR, XRI);
    prtmat("transXR", XRI);
    System.out.println("check XRT = transXR");
    System.out.println(" ");
    System.out.println(" ");
    
    System.out.println("rotate around X axis using XR, inverse");    
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    makeXRT(ax, XR);
    prtmat("XR", XR);
    matvec(XR, p1234, p);
    prtpt("XR*p1234", p);
    System.out.println("reverse rotate using inverse XR");
    matinv(XR, XRI); // same angle, inverse matrix
    prtmat("XRI", XRI);
    matvec(XRI, p, p1);
    prtpt("XRI*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("rotate around X axis using XR, negative ax");    
    System.out.println("ax="+ax+" rad,  ax="+((ax/Pi)*180.0)+" deg");  
    makeXR(ax, XR);
    prtmat("XR", XR);
    matvec(XR, p1234, p);
    prtpt("XR*p1234", p);
    System.out.println("reverse rotate using -ax XR");
    makeXR(-ax, XR);
    prtmat("XR", XR);
    matvec(XR, p, p1);
    prtpt("XR*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("rotate around Y axis");
    makeYR(ay, YR);
    System.out.println("ay="+ay+" rad,  ay="+((ay/Pi)*180.0)+" deg");  
    prtmat("YR", YR);
    matvec(YR, p1234, p);
    prtpt("YR*p1234", p);
    System.out.println("reverse rotate YR");
    makeYR(-ay, YR);
    prtmat("YR", YR);
    matvec(YR, p, p1);
    prtpt("YR*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("rotate around Z axis");
    makeZR(az, ZR);
    System.out.println("az="+az+" rad,  az="+((az/Pi)*180.0)+" deg");  
    prtmat("ZR", ZR);
    matvec(ZR, p1234, p);
    prtpt("ZR*p1234", p);
    System.out.println("reverse rotate ZR");
    makeZR(-az, ZR);
    prtmat("ZR", ZR);
    matvec(ZR, p, p1);
    prtpt("ZR*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    System.out.println("rotate around W axis");
    makeWR(aw, WR);
    System.out.println("aw="+aw+" rad,  aw="+((aw/Pi)*180.0)+" deg");  
    prtmat("WR", WR);
    matvec(WR, p1234, p);
    prtpt("WR*p1234", p);
    System.out.println("reverse rotate WR");
    makeWR(-aw, WR);
    prtmat("WR", WR);
    matvec(WR, p, p1);
    prtpt("WR*p", p1);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ay; // 30 deg
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR30",p);
    ang = ang + ay;
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR60",p);
    ang = ang + ay;
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR90",p);
    ang = ang + ay;
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR120",p);
    ang = ang + ay;
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR150",p);
    ang = ang + ay;
    makeXR(ang,XR);
    matvec(XR,p1234,p);
    prtpt("XR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ay; // 30 deg
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR30",p);
    ang = ang + ay;
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR60",p);
    ang = ang + ay;
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR90",p);
    ang = ang + ay;
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR120",p);
    ang = ang + ay;
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR150",p);
    ang = ang + ay;
    makeYR(ang,YR);
    matvec(YR,p1234,p);
    prtpt("YR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ay; // 30 deg
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR30",p);
    ang = ang + ay;
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR60",p);
    ang = ang + ay;
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR90",p);
    ang = ang + ay;
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR120",p);
    ang = ang + ay;
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR150",p);
    ang = ang + ay;
    makeZR(ang,ZR);
    matvec(ZR,p1234,p);
    prtpt("ZR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ay; // 30 deg
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR30",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR60",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR90",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR120",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR150",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR210",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR240",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR270",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR300",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR330",p);
    ang = ang + ay;
    makeWR(ang,WR);
    matvec(WR,p1234,p);
    prtpt("WR360",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ay; // 30 deg
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR30",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR60",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR90",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR120",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR150",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR210",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR240",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR270",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR300",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR330",p);
    ang = ang + ay;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR360",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ax; // 15 deg
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR15",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR30",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR45",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR60",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR75",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR90",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR105",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR120",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR135",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR150",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR165",p);
    ang = ang + ax;
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");

    ang = ax; // 15 deg
    makeXYZWR(ang, ang, ang, ang, XYZWR);
    matvec(XYZWR,p1234,p);
    prtpt("XYZWR15",p);
    matvec(XYZWR,p,p1);
    prtpt("XYZWR30",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR45",p);
    matvec(XYZWR,p,p1);
    prtpt("XYZWR60",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR75",p);
    matvec(XYZWR,p,p1);
    prtpt("XYZWR90",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR105",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR120",p);
    matvec(XYZWR,p,p1);
    prtpt("XYZWR135",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR150",p);
    matvec(XYZWR,p,p1);
    prtpt("XYZWR165",p1);
    matvec(XYZWR,p1,p);
    prtpt("XYZWR180",p);
    System.out.println("check == p1234");
    System.out.println(" ");
    System.out.println(" ");    
  } // end test_matrix
	
  class MyThread extends Thread
  {
    MyThread()
    {
    } // end MyThread constructor
    
    public void run()
    {
      repaint(); // dummy in 2D
    } // end Run
  } // end class MyThread

  void sort(double val[], int n, int index[])
  {
    int k;
    double t;
    
    for(int i=0; i<n; i++) index[i] = i;
    for(int i=0; i<n-1; i++)
    {
      for(int j=i+1; j<n; j++)
      {
	if(val[j]<val[i])
	{
	  t = val[i];
	  val[i] = val[j];
	  val[j] = t;
	  k = index[i];
	  index[i] = index[j];
	  index[j] = k;
	} // end if
      } // end j
    } // end i
  } // end sort
    
  public static void main(String[] args)
  {
    String filename = args[0];
    new plot44(filename);
  }
} // end class  plot44.java