// light_stl.java  3D projection of a  .stl  file
//                                needs  center
// needs heapsort.java compiled

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


// solid rect
//   facet normal 0.0 0.0 -1.0
//     outer loop
//       vertex 1.0 1.0 1.0
//       vertex 5.0 4.0 1.0
//       vertex 1.0 4.0 1.0
//     endloop
//   endfacet
//   facet normal 0.0 0.0 -1.0
//     outer loop
//       vertex 1.0 1.0 1.0
//       vertex 5.0 1.0 1.0
//       vertex 5.0 4.0 1.0
//    endloop
//   endfacet
//...
// endsolid

public class light_stl  extends Frame
{
  // input .slt  convert to .dat format internally 
  int num_points;
  int num_polys;
  double x[] = new double[60000];
  double y[] = new double[60000];
  double z[] = new double[60000];
  int  data_polys[]  = new int[1000000];
  int pripoly[] = new int[40000];
  int priindex[] = new int[40000];
  double priz[] = new double[40000]; 
  double winxyz[] = new double[3];
  double xmin, xmax, ymin, ymax, zmin, zmax;
  double dxmin, dxmax, dymin, dymax, dzmin, dzmax;
  double dxinc, dyinc, dzinc;
  String fname;
  double size = 1.0;
  double gsize = 1.0;
  double nsize = -1.0;
  double scale = 200.0;
  boolean debug = false;
  int xm, ym, bm; // mouse
  double pos[] = {0.0, 0.0, 0.0};
  double roll = 0.0;  // rotate object 
  double pitch = 0.0;
  double heading = 0.0;
  int wire = 1; // ??? set to zero when debugged 
  int vert = 0;
  int colorz = 0;
  int noden = 1;
  double mmatrix[][] = new double[4][4];// 4 by 4 model transformation matrix 
  double pmatrix[][] = new double[4][4];// 4 by 4 persp transformation matrix 
  int viewport[] = new int[4];
  double pxyzw[] = new double[4];
  int win_w = 600;
  int win_h = 600;
  double M_PI = Math.PI;
  heapsort H = new heapsort();

  light_stl(String filename)
  {
    System.out.println("light_stl.java reading "+filename);
    read_data(filename);
    System.out.println("num_points="+num_points+", num_polys="+num_polys);  

    fname = filename;
    setTitle("light_stl");
    setSize(win_w,win_h);
    setBackground(Color.white);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    this.addMouseListener (new mousePressHandler());
  } // end light_stl

  void read_data(String filein)
  {

    int npt = 0;
    int pts[] = new int[3];
    int ipt; // unique     
    StringTokenizer T;
    String token;
    String solid = "solid";
    String facet = "facet";
    String normal = "normal";
    String outer = "outer";
    String loop = "loop";
    String vertex = "vertex";
    String endloop = "endloop";
    String endfacet = "endfacet";
    String endsolid = "endsolid";
      
    if(debug) System.out.println("light_stl.java reading "+filein);
    try
    {
      FileReader fin = new FileReader(filein);
      BufferedReader in = new BufferedReader(fin);
      String input_line;
      int len;

      while(true) // to  endsolid
      {
        input_line = in.readLine(); // 
        len = input_line.length();
        if(debug) System.out.println("input:"+input_line+"  len="+len);
        T = new StringTokenizer(input_line," ");
        token = T.nextToken(); // ?
        if(debug) System.out.println("token:"+token);
	if(token.equals(endsolid)) break;
	if(token.equals(endloop))
	{
	  data_polys[num_polys] = 3;
	  data_polys[num_polys+1] = pts[0];
	  data_polys[num_polys+2] = pts[1];
	  data_polys[num_polys+3] = pts[2];
	  num_polys += 4;
	  npt = 0;
	  continue;
	}
	if(token.equals(endfacet)) continue;
	if(token.equals(vertex))
	{
          token = T.nextToken(); // x
          if(debug) System.out.println("token:"+token);
          x[num_points] = Double.parseDouble(token);
          token = T.nextToken(); // y
          if(debug) System.out.println("token:"+token);
          y[num_points] = Double.parseDouble(token);
          token = T.nextToken(); // z
          if(debug) System.out.println("token:"+token);
          z[num_points] = Double.parseDouble(token);
	  num_points++;
	  ipt = unique();
	  pts[npt] = ipt; // one big per  .dat
	  npt++;
	}
      } // end while      
    } // end try
    catch(IOException exception)
    {
      System.out.println("datread.java "+exception);
    }
    System.out.println("unique data points");
    xmin = x[0];
    xmax = x[0];
    ymin = y[0];
    ymax = y[0];
    zmin = z[0];
    zmax = z[0];
    for(int i=0; i<num_points; i++)
    {
      xmin = Math.min(xmin, x[i]);
      xmax = Math.max(xmax, x[i]);
      ymin = Math.min(ymin, y[i]);
      ymax = Math.max(ymax, y[i]);
      zmin = Math.min(zmin, z[i]);
      zmax = Math.max(zmax, z[i]);
      if(debug) System.out.println("x="+x[i]+", y="+y[i]+", z="+z[i]);
    }
    size = Math.max(xmax, ymax);
    size = Math.max(size, zmax);
    nsize = Math.min(xmin, ymin);
    nsize = Math.min(nsize, zmin);
    size = Math.max(size, Math.abs(-nsize));
    size = size/2.0;
    dxmin = xmin*scale;
    dxmax = xmax*scale;
    dymin = ymin*scale;
    dymax = ymax*scale;
    dzmin = zmin*scale;
    dzmax = zmax*scale;
    dxinc = 1.0;
    dyinc = 1.0;
    dzinc = 1.0;
    gsize = Math.max((dxmax-dxmin),(dymax-dymin));
    gsize = Math.max(gsize,(dzmax-dzmin));
    gsize = gsize / 100.0;
    System.out.println("size="+size+", gsize="+gsize);

    if(debug) System.out.println("polys");
    for(int i=0; i<num_polys; i=i+4)
    {
	if(debug) System.out.println("pt "+(i/4)+"  "+data_polys[i]+"  "+
				     data_polys[i+1]+"  "+data_polys[i+2]+
				     "  "+data_polys[i+3]);
    } // end i
    num_polys = num_polys/4;
    System.out.println("num_points="+num_points+", num_polys="+num_polys);
    System.out.println("xmin="+xmin+", xmax="+xmax+
		       "\nymin="+ymin+", ymax="+ymax+
		       "\nzmin="+zmin+", zmax="+zmax);
  } // end read_data

  int unique() // keep unique x,y,z
  {
    if(num_points<2) return 1;
    for(int i=0; i<num_points-1; i++)
    {
      if(Math.abs(x[num_points-1]-x[i])<0.001 &&
	 Math.abs(y[num_points-1]-y[i])<0.001 &&
	 Math.abs(z[num_points-1]-z[i])<0.001)
      {
	num_points--;
	return i+1; // one too big
      }
    }
    return num_points;
  } // end unique
    
  public void paint(Graphics g)
  {
    int x1=0, x2=0, y1=0, y2=0;

    g.setColor(Color.black);
    g.drawString("xmin="+(dxmin/scale),  50,50);
    g.drawString("xmax="+(dxmax/scale),  50,65);
    g.drawString("ymin="+(dymin/scale), 250,50);
    g.drawString("ymax="+(dymax/scale), 250,65);
    g.drawString("zmin="+(dzmin/scale), 450,50);
    g.drawString("zmax="+(dzmax/scale), 450,65);
    g.drawString("light_stl for 3D surface on "+fname, 20, 20);
    g.drawString("C colored faces, W wireframe, V vertices, Nodes sz SZ size ", 20, 35);
    g.drawString("< X >  v Y ^  < Z >  < R >  ^ P v  v H ^    W    V    C    N   sz  SZ",
                 20, win_h-25);
    g.drawString("x="+pos[0]+", y="+pos[1]+", z-"+pos[2]+", R="+roll+", P="+
                 pitch+", H="+heading+", S="+size, 20, win_h-10);
    order_polys();

    if(colorz==1)    draw_face_dat(g);
    if(wire==1)      draw_wire_dat(g);
    if(noden==1)     draw_node_dat(g);
    else if(vert==1) draw_vert_dat(g);
    else             draw_face_dat(g);

  } // end paint

  class mousePressHandler extends MouseAdapter
  {
    public void mousePressed (MouseEvent e)
    {
      char key=' ';
      xm = e.getX();
      ym = e.getY();
      bm = e.getButton();
      System.out.println("mouse xm="+xm+", ym="+ym+", bm="+bm);
      if(ym>win_h-21 || ym<win_h-36) return;
      if(xm>18 && xm<28) key = 'x';
      if(xm>45 && xm<55) key = 'X';
      if(xm>62 && xm<72) key = 'y';
      if(xm>85 && xm<95) key = 'Y';
      if(xm>102 && xm<112) key = 'z';
      if(xm>125 && xm<135) key = 'Z';
      if(xm>148 && xm<158) key = 'r';
      if(xm>172 && xm<182) key = 'R';
      if(xm>191 && xm<201) key = 'p';
      if(xm>214 && xm<224) key = 'P';
      if(xm>230 && xm<240) key = 'h';
      if(xm>255 && xm<265) key = 'H';
      if(xm>275 && xm<295) key = 'W';
      if(xm>300 && xm<320) key = 'V';
      if(xm>325 && xm<345) key = 'C';
      if(xm>345 && xm<370) key = 'N';
      if(xm>370 && xm<400) key = 's';
      if(xm>400 && xm<430) key = 'S';
      if(debug)
      {
        System.out.println("at xm="+xm+", ym="+ym+", bm="+bm); // debug print
        System.out.println("key="+key);
      }
      if(key == 'x') {pos[0] -= dxinc;} // position object 
      if(key == 'X') {pos[0] += dxinc;}
      if(key == 'y') {pos[1] -= dyinc;} // fix Y down
      if(key == 'Y') {pos[1] += dyinc;}
      if(key == 'z') {pos[2] -= dzinc;}
      if(key == 'Z') {pos[2] += dzinc;}
      if(key == 'r') {roll   -= 15.0;  // rotate object 
	              if(roll<=-360) roll += 360.0;}
      if(key == 'R') {roll   += 15.0;
	              if(roll>=360) roll -= 360.0;}
      if(key == 'p') {pitch  -= 15.0;
	              if(pitch<=-360) pitch += 360.0;}
      if(key == 'P') {pitch  += 15.0;
	              if(pitch>=360) pitch -= 360.0;}
      if(key == 'h') {heading -= 15.0;
	              if(heading<=-360) heading += 360.0;}
      if(key == 'H') {heading += 15.0;
	              if(heading>=360) heading -= 360.0;}
      if(key == 'w' || key == 'W') {wire = 1-wire;} // 0,0,0 is colorz face
      if(key == 'v' || key == 'V') {vert = 1-vert; wire=0; colorz=0;}
      if(key == 'c' || key == 'C') {colorz = 1-colorz; wire=0; vert=0;}
      if(key == 'n' || key == 'N') {noden = 1-noden;}
      if(key == 's') {gsize = gsize*2.0;} // size smaller
      if(key == 'S') {gsize = gsize*0.5;} // size larger

      System.out.println("roll="+roll+", pitch="+pitch+", head="+heading+
                         ", size="+gsize); // debug print
      requestFocus();
      repaint();
    }
  } // end mousePressHandler


  void draw_face_dat(Graphics g)
  {
    int i, j, k, pts, pt, pt1;
    int skipz=0;
    double nxyz[] = new double[3];
    int xpts[] = new int[10];
    int ypts[] = new int[10];
    double xsav[] = new double[3];
    double ysav[] = new double[3];
    double zsav[] = new double[3];
    double zmin=0.0, zmax=0.0, zavg=0.0;
    double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4 ;

    build_mmat(); // from pitch, heading, roll, pos
    //glGetDoublev(GL_MODELVIEW_MATRIX, mmatrix); // save transformation 
    //glGetDoublev(GL_PROJECTION_MATRIX, pmatrix);

    if(debug) System.out.println("in draw_face_dat");

    // use num_polys to do that many faces 
    k = 0;
    for(i=0; i<num_polys; i++)
    {
      k = pripoly[priindex[i]]; // biggest z first
      pts = data_polys[k++];
      if(debug) System.out.println("i="+i+", wire pts="+pts);
      // draw
      pt1 = data_polys[k]; // save first point
      for(j=0; j<pts; j++)
      {
        pt = data_polys[k++];
        x1 = x[pt-1];
        y1 = y[pt-1];
        z1 = z[pt-1];
        if(j<3)
	{
	  xsav[j] = x1; // for normal
	  ysav[j] = y1;
	  zsav[j] = z1;
	}
        project(x1, y1, z1, winxyz);
        x1 = winxyz[0];
        y1 = winxyz[1];
        z1 = winxyz[2];
        xpts[j] = (int)winxyz[0];
        ypts[j] = win_h-(int)winxyz[1];
        if(debug)
        {
          System.out.println("face . . . . . . . . . . . . . . . .");
          System.out.println("fr point="+pt+", x1="+x[pt-1]+", y1="+y[pt-1]+
                                          ", z1="+z[pt-1]);
          System.out.println("fr winxyz"+j+", xw="+winxyz[0]+", yw="+winxyz[1]+
                             ", zw="+winxyz[2]);
	}
      } // end j

      // normal sets color
      x1 = xsav[0];
      y1 = ysav[0];
      z1 = zsav[0];
      x2 = xsav[1];
      y2 = ysav[1];
      z2 = zsav[1];
      x3 = xsav[2];
      y3 = ysav[2];
      z3 = zsav[2];
      normal(x1,y1,z1, x2,y2,z2, x3,y3,z3, g);
      g.fillPolygon(xpts, ypts, pts);

    } // end loop on num_polys 
  } // end draw_face_dat 

  void draw_wire_dat(Graphics g)
  {
    int i, j, k, pts, pt, pt1;
    int skipz=0;
    double nxyz[] = new double[3];
    double zmin=0.0, zmax=0.0, zavg=0.0;
    double x1, x2, y1, y2, z1, z2;

    build_mmat(); // from pitch, heading, roll, pos
    //glGetDoublev(GL_MODELVIEW_MATRIX, mmatrix); // save transformation 
    //glGetDoublev(GL_PROJECTION_MATRIX, pmatrix);

    if(debug) System.out.println("in draw_wire_dat");

    if(colorz==1)
    {
      zavg = 0.0;
      for(i=0; i<num_points; i++)
      {
        project(x[i], y[i], z[i], winxyz);
        zavg = zavg + winxyz[2];
	if(i==0) { zmin = winxyz[2]; zmax = winxyz[2]; }
	if(winxyz[2] < zmin) zmin = winxyz[2];
        if(winxyz[2] > zmax) zmax = winxyz[2];		
      }
      zavg = zavg/ (double)num_points;
    }

    // use num_polys to do wire frame 
    k = 0;
    for(i=0; i<num_polys; i++)
    {
      pts = data_polys[k++];
      if(debug) System.out.println("i="+i+", wire pts="+pts);
      // draw
      pt1 = data_polys[k]; // save first point
      for(j=0; j<pts; j++)
      {
        pt = data_polys[k++];
        x1 = x[pt-1];
        y1 = y[pt-1];
        z1 = z[pt-1];
        if(debug) System.out.println("j="+j+", pt="+pt);
        project(x1, y1, z1, winxyz);
        x1 = winxyz[0];
        y1 = winxyz[1];
        z1 = winxyz[2];
        if(debug)
        {
          System.out.println("wire . . . . . . . . . . . . . . . .");
          System.out.println("fr point="+pt+", x1="+x[pt-1]+", y1="+y[pt-1]+
                                          ", z1="+z[pt-1]);
          System.out.println("fr winxyz"+j+", xw="+winxyz[0]+", yw="+winxyz[1]+
                             ", zw="+winxyz[2]);
	}
        if(j==pts-1) // use first point to close
	{
          x2 = x[pt1-1];
          y2 = y[pt1-1];
          z2 = z[pt1-1];
          if(debug) System.out.println("tp pointcl="+pt1+", x2="+x2+", y2="+y2+
                                       ", z2="+z2);
	}
	else
	{
	  pt = data_polys[k]; // was incremented
          x2 = x[pt-1];
          y2 = y[pt-1];
          z2 = z[pt-1];
          if(debug) System.out.println("to point="+pt+", x2="+x2+", y2="+y2+
                                       ", z2="+z2);
	}
        project(x2, y2, z2, winxyz);
        x2 = winxyz[0];
        y2 = winxyz[1];
        z2 = winxyz[2];
        if(debug) System.out.println("to winxyz"+j+", xw="+winxyz[0]+
                                     ", yw="+winxyz[1]+ ", zw="+winxyz[2]);
        g.setColor(Color.black);
        g.drawLine((int)x1, win_h-(int)y1, (int)x2, win_h-(int)y2);
      }

    } // end loop on num_polys 
  } // end draw_wire_dat 


  void normal(double x0, double y0, double z0,
              double x1, double y1, double z1,
              double x2, double y2, double z2,
              Graphics g)
  {
    double ax, ay, az, bx, by, bz, nx, ny, nz, denom;
    float red, grn, blu; 
    ax = x2 - x0;
    ay = y2 - y0;
    az = z2 - z0;
    bx = x1 - x0;
    by = y1 - y0;
    bz = z1 - z0;
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;

    denom = Math.sqrt(nx*nx+ny*ny+nz*nz);
    nx = nx/denom;
    ny = ny/denom;
    nz = nz/denom;
    // System.out.println(x0+y0+z0);
    // System.out.println(x1+y1+z1);
    // System.out.println(x2+y2+z2);
    // System.out.println("nx="+nx+", ny="+ny+", nz="+nz);
    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));
  } // end normal


  void draw_vert_dat(Graphics g)
  {
    int i, j, k, pts, pt;
    double RGB[] = new double[3];
    double zc;
    double zavg, zmin=0.0, zmax=0.0;

    if(colorz==1)
    {
      zavg = 0.0;
      for(i=0; i<num_points; i++)
      {
        project(x[i], y[i], z[i], winxyz);
	zavg = zavg + winxyz[2];
	if(i==0) { zmin = winxyz[2]; zmax = winxyz[2]; }
	if(winxyz[2] < zmin) zmin = winxyz[2];
        if(winxyz[2] > zmax) zmax = winxyz[2];		
      }
      zavg = zavg/ (double)num_points;
    }

    // use num_polys to do vertices 
    k = 0;
    for(i=0; i<num_polys; i++)
    {
      pts = data_polys[k++];
      if(debug) System.out.println("vert pts="+pts);
      for(j=0; j<pts; j++)
      {
        pt = data_polys[k++];
        if(debug) System.out.println("j="+j+", wire look pt="+pt);
        project(x[pt-1], y[pt-1], z[pt-1], winxyz);
        if(colorz==1)
	{
	  zc =  winxyz[2];
          // colorf((zc-zmin)/(zmax-zmin), RGB);
          // setColorf(RGB);
        }
        g.setColor(Color.green);
        g.fillOval((int)winxyz[0]-2, win_h-(int)winxyz[1]-2, 4, 4);
      }
    } // end loop on num_polys 
  } // end draw_vert_dat 


  void draw_node_dat(Graphics g)
  {
    int i, j, k, pts, pt;
    // use num_polys to do vertices 
    k = 0;
    for(i=0; i<num_polys; i++)
    {
      pts = data_polys[k++];
      if(debug) System.out.println("vert pts="+pts);
      for(j=0; j<pts; j++)
      {
        pt = data_polys[k++];
        if(debug) System.out.println("j="+j+", wire look pt="+pt);
        project(x[pt-1], y[pt-1], z[pt-1], winxyz);
        g.setColor(Color.blue);
        g.drawString(" "+pt, (int)winxyz[0]+1, win_h-(int)winxyz[1]+1);
      }
    } // end loop on num_polys 
  } // end draw_node_dat 

  void build_mmat() // pitch, heading, roll, pos);
  {
    double rmatrix[][] = new double[4][4];
    double tmatrix[][] = new double[4][4];
    double c[][] = new double[4][4];
    double p[] = {0.0, 0.0, 0.0, 1.0};
    double pp[] = new double[4];

    // generate mmatrix and pmatrix
    final double imatrix[][] =  {{1.0, 0.0, 0.0, 0.0},
                                 {0.0, 1.0, 0.0, 0.0},
  	       	                 {0.0, 0.0, 1.0, 0.0},
		                 {0.0, 0.0, 0.0, 1.0}};
    double pll[] = new double[4];
    double pur[] = new double[4];
    double pllf[] = new double[4];
    double purf[] = new double[4];
    double xmin, xmax, ymin, ymax, near, far;
    double w=win_w;
    double h=win_h;
    double xt, yt, zt;

    xmin=-gsize;
    xmax= gsize;
    ymin=-gsize*h/w;
    ymax= gsize*h/w;
    near= 4.0*gsize;
    far = 20.0*gsize;

    if(debug)
    {
      System.out.println("in build_mmat");
      System.out.println("xmin="+xmin+", xmax="+xmax);
      System.out.println("ymin="+ymin+", ymax="+ymax);
      System.out.println("near="+near+", far ="+far);
    }
    pll[0]=xmin;        pll[1]=ymin;        pll[2]=-near; pll[3]=1.0;
    pur[0]=xmax;        pur[1]=ymax;        pur[2]=-near; pur[3]=1.0;
    pllf[0]=100.0*xmax; pllf[1]=100.0*ymin; pllf[2]=-far; pllf[3]=1.0;
    purf[0]=100.0*xmin; purf[1]=100.0*ymax; purf[2]=-far; purf[3]=1.0;
  
    pmatrix[0][0]=2.0*near/(xmax-xmin);
    pmatrix[0][1]=0.0;
    pmatrix[0][2]=(xmax+xmin)/(xmax-xmin);
    pmatrix[0][3]=0.0;
    pmatrix[1][0]=0.0;
    pmatrix[1][1]=2.0*near/(ymax-ymin);
    pmatrix[1][2]=(ymax+ymin)/(ymax-ymin);
    pmatrix[1][3]=0.0;
    pmatrix[2][0]=0.0;
    pmatrix[2][1]=0.0;
    pmatrix[2][2]=-(far+near)/(far-near);
    pmatrix[2][3]=-2.0*far*near/(far-near);
    pmatrix[3][0]=0.0;
    pmatrix[3][1]=0.0;
    pmatrix[3][2]=-1.0;
    pmatrix[3][3]=0.0;

    matcpy(imatrix, rmatrix);
    rmatrix[1][1]= Math.cos(M_PI*pitch/180.0);
    rmatrix[2][2]= Math.cos(M_PI*pitch/180.0);
    rmatrix[1][2]=-Math.sin(M_PI*pitch/180.0);
    rmatrix[2][1]= Math.sin(M_PI*pitch/180.0);
    matmul(pmatrix, rmatrix, c);
    matcpy(c, mmatrix);

    matcpy(imatrix, rmatrix);
    rmatrix[0][0]= Math.cos(M_PI*heading/180.0);
    rmatrix[2][2]= Math.cos(M_PI*heading/180.0);
    rmatrix[0][2]=-Math.sin(M_PI*heading/180.0);
    rmatrix[2][0]= Math.sin(M_PI*heading/180.0);
    matmul(mmatrix, rmatrix, c);
    matcpy(c, mmatrix);

    matcpy(imatrix, rmatrix);
    rmatrix[0][0]= Math.cos(M_PI*roll/180.0);
    rmatrix[1][1]= Math.cos(M_PI*roll/180.0);
    rmatrix[0][1]=-Math.sin(M_PI*roll/180.0);
    rmatrix[1][0]= Math.sin(M_PI*roll/180.0);
    matmul(mmatrix, rmatrix, c);
    matcpy(c, mmatrix);

    matcpy(imatrix, tmatrix);
    tmatrix[0][3] = pos[0];
    tmatrix[1][3] = pos[1];
    tmatrix[2][3] = pos[2];
    if(debug)
    {
      System.out.println("\n tmatrix");
      prmat(tmatrix);
    }
    matmul(mmatrix, tmatrix, c);
    matcpy(c, mmatrix);
    if(debug)
    {
      System.out.println("\n mmatrix rotated and translated \n");
      System.out.println("\n model matrix");
      prmat(mmatrix);
      System.out.println("\n projection matrix");
      prmat(pmatrix);
    }
  }

  void project(double xx, double yy, double zz, double wxyz[])
  {
    double c[][] = new double[4][4];
    double p[] = new double[4];
    double pp[] = new double[4];
    double ax, ay, az;
    // use mmatrix pmatrix

    p[0] = xx;
    p[1] = yy;
    p[2] = zz;
    p[3] = 1.0;
    matmul(pmatrix, mmatrix, c);
    vecmul(c, p, pp);
    ax = (double)win_w*(1.0+pp[0]/pp[3])/2.0;
    ay = (double)win_h*(1.0+pp[1]/pp[3])/2.0;
    az = (1.0+pp[2]/pp[3])/2.0;
    if(debug) System.out.println("x scr="+ax+", y scr="+ay+", at z="+az);
    wxyz[0] = ax;
    wxyz[1] = ay;
    wxyz[2] = az;
  } // end project

  void prmat(double matrix[][])
  {
    System.out.println("matrix[0][0..3]="+matrix[0][0]+", "+
                        matrix[0][1]+", "+matrix[0][2]+", "+matrix[0][3]);
    System.out.println("matrix[1][0..3]="+matrix[1][0]+", "+
                        matrix[1][1]+", "+matrix[1][2]+", "+matrix[1][3]);
    System.out.println("matrix[2][0..3]="+matrix[2][0]+", "+
                        matrix[2][1]+", "+matrix[2][2]+", "+matrix[2][3]);
    System.out.println("matrix[3][0..3]="+matrix[3][0]+", "+
                        matrix[3][1]+", "+matrix[3][2]+", "+matrix[3][3]);
  } // end prmat

  void matmul(double a[][], double b[][], double c[][])
  {
    for(int i=0; i<4; i++)
    {
      for(int j=0; j<4; j++)
      {
        c[i][j]=0.0;
        for(int k=0; k<4; k++)
        {
          c[i][j] = c[i][j] + a[i][k] * b[k][j];
        }
      }
    }
  } // end matmul

  void matcpy(double a[][], double c[][])
  {
    for(int i=0; i<4; i++)
    {
      for(int j=0; j<4; j++)
      {
        c[i][j] = a[i][j];
      }
    }
  } // end matcpy


  void vecmul(double a[][], double p[], double pp[])
  {
    for(int i=0; i<4; i++)
    {
      pp[i] = 0.0;
      for(int j=0; j<4; j++)
      {
        pp[i]= pp[i] + a[i][j] * p[j];
      }
    }
  } // end vecmul

  void order_polys()  // compute and sort pripoly[], priz[]
  {
    int k, pts, pt;
    double x1, y1, z1, zavg;

    k = 0;
    for(int i=0; i<num_polys; i++)
    {
      pripoly[i] = k;
      pts = data_polys[k++];
      zavg = 0.0;
      for(int j=0; j<pts; j++)
      {
        pt = data_polys[k++];
        x1 = x[pt-1];
        y1 = y[pt-1];
        z1 = z[pt-1];
        project(x1, y1, z1, winxyz);
        x1 = winxyz[0];
        y1 = winxyz[1];
        z1 = winxyz[2];
        zavg = zavg + z1;
      } // end j
      zavg = zavg/(double)pts;
      priz[i] = zavg;
    } // end i
    H.sortdindex(priz, num_polys, priindex);
    for(int i=0; i<num_polys; i++)
	if(debug) System.out.println("poly="+priindex[i]+", k="+
				     pripoly[priindex[i]]+", z="+priz[i]);
  } // end order_polys

  public static void main(String[] args)
  {
      new light_stl(args[0]); // e.g. rect.stl
  }  
} // end class light_stl