// light_dat.java  projection may not be best
//                                needs  center
// needs heapsort.java compiled

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

public class light_dat  extends Frame
{
  // .dat format file input
  int num_points;
  int num_polys;
  double x[] = new double[20000];
  double y[] = new double[20000];
  double z[] = new double[20000];
  int  data_polys[]  = new int[256000];
  int pripoly[] = new int[18000];
  int priindex[] = new int[18000];
  double priz[] = new double[18000]; 
  double winxyz[] = new double[3];
  int maxpolys;
  double dxmin, dxmax, dymin, dymax, dzmin, dzmax;
  double dxctr, dyctr, dzctr, dxinc, dyinc, dzinc;
  String fname;
  double size = 1.0;
  double gsize = 1.0;
  double nsize = -1.0;
  double scale = 0.15;
  double  offs = 1.0; // not used
  boolean debug = true;
  int xm, ym, bm; // mouse
  double pos[] = {0.0, 0.0, 0.0};
  // keys or mouse clicks
  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;
  int nonum = 0;
  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_dat(String filename)
  {
    System.out.println("light_dat.java reading "+filename);
    read_data(filename);
    System.out.println("num_points="+num_points+", num_polys="+num_polys);  
    System.out.println("size="+size+", maxpolys="+maxpolys);

    fname = filename;
    setTitle("light_dat");
    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_dat

  void read_data(String filename)
  {
    String input_line;
    int k=0;
    int ntok, pts;
    StringTokenizer T;
    String token[] = new String[255];

    try
    {
      FileReader fin  = new FileReader(filename);
      BufferedReader file_in = new BufferedReader(fin);

      if(debug) System.out.println("read_data reading "+filename);

      input_line = file_in.readLine();
      if(debug) System.out.println(input_line);
      T = new StringTokenizer(input_line," ");
      ntok=0;
      while(T.hasMoreTokens())
      {
        token[ntok]=T.nextToken();
        if(debug) System.out.println("token["+ntok+"]=|"+token[ntok]+"|");
        ntok++;
      }
      num_points = Integer.parseInt(token[0]);
      num_polys  = Integer.parseInt(token[1]);
      System.out.println("num_points="+num_points+", num_polys="+num_polys);

      // read points
      for(int i=0; i<num_points; i++)
      {
	// System.out.println("num_points i="+i);
        input_line = file_in.readLine();
        if(i<10 && debug) System.out.println(input_line);
        T = new StringTokenizer(input_line," ");
        ntok=0;
        while(T.hasMoreTokens())
        {
          token[ntok]=T.nextToken();
          if(i<5 && debug) System.out.println("token["+ntok+"]=|"+token[ntok]+"|");
          ntok++;
	}
        x[i] = Double.parseDouble(token[0])*scale;
        y[i] = Double.parseDouble(token[1])*scale;
        z[i] = Double.parseDouble(token[2])*scale;
        size = Math.max(size, x[i]);
        size = Math.max(size, y[i]);
        size = Math.max(size, z[i]);
        nsize = Math.min(nsize, x[i]);
        nsize = Math.min(nsize, y[i]);
        nsize = Math.min(nsize, z[i]);
      }
      // System.out.println("points read, read poly");
      size = Math.max(size, Math.abs(-nsize));
      size = size/2.0;
      System.out.println("size/2 ="+size);
      
      // read polygons
      for(int i=0; i<num_polys; i++)
      {
        input_line = file_in.readLine();
        if(debug) System.out.println(input_line);
        T = new StringTokenizer(input_line," ");
        ntok=0;
        while(T.hasMoreTokens())
        {
          token[ntok]=T.nextToken();
          if(debug) System.out.println("token["+ntok+"]=|"+token[ntok]+"|");
          ntok++;
	}
        pts = Integer.parseInt(token[0]);
        if(ntok != pts+1 && debug) System.out.println("poly="+i+", pts="+pts);
        data_polys[k] = pts;
        if(i<5 && debug) System.out.println("i="+i+", pts="+pts);
        k++;
        for(int j=0; j<pts; j++)
	{
          data_polys[k] = Integer.parseInt(token[j+1]);
          if(i<5 && debug) System.out.println("j="+j+", pt="+data_polys[k]);
          k++;
	}
        maxpolys = k;
	// System.out.println("k="+k);
      }
      file_in.close();
      if(debug)System.out.println("finished read_data");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println(filename+" not found");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("java exception in read_data"+exception);
    }
    dxmin = 10000.0;
    dxmax = -10000.0;
    dymin = 10000.0;
    dymax = -10000.0;
    dzmin = 10000.0;
    dzmax = -10000.0;
    for(int i=0; i<num_points; i++)
    {
      dxmin = Math.min(dxmin,x[i]); // allready scaled
      dxmax = Math.max(dxmax,x[i]); // size, nsize allready
      dymin = Math.min(dymin,y[i]);
      dymax = Math.max(dymax,y[i]);
      dzmin = Math.min(dzmin,z[i]);
      dzmax = Math.max(dzmax,z[i]);
    }
    System.out.println("dxmin="+(dxmin/scale)+", dxmax="+(dxmax/scale));
    System.out.println("dymin="+(dymin/scale)+", dymax="+(dymax/scale));
    System.out.println("dzmin="+(dzmin/scale)+", dzmax="+(dzmax/scale));
    dxctr = (dxmax+dxmin)/2.0;
    dyctr = (dymax+dymin)/2.0;
    dzctr = (dzmax+dzmin)/2.0;
    System.out.println("dxctr="+dxctr+", dyctr="+dyctr+", dzctr="+dzctr);
    dxinc = (dxmax-dxmin)/10.0;
    dyinc = (dymax-dymin)/10.0;
    dzinc = (dzmax-dzmin)/10.0;

    gsize = Math.max((dxmax-dxmin),(dymax-dymin));
    gsize = Math.max(gsize,(dzmax-dzmin));
    gsize = gsize * 1.25;

    if(debug) System.out.println("dxinc="+dxinc+", dyinc="+dyinc+
                                 ", dzinc="+dzinc);
  } // end read_data


  public void paint(Graphics g)
  {
    int x1=0, x2=0, y1=0, y2=0;

    g.setColor(Color.black);
    g.drawString("xmin="+(dxmin/scale),  50,80);
    g.drawString("xmax="+(dxmax/scale),  50,95);
    g.drawString("ymin="+(dymin/scale), 200,80);
    g.drawString("ymax="+(dymax/scale), 200,95);
    g.drawString("zmin="+(dzmin/scale), 350,80);
    g.drawString("zmax="+(dzmax/scale), 350,95);
    g.drawString("vertex="+num_points,  450,80);
    g.drawString("polygons="+num_polys, 450,95);
    g.drawString("light_dat.java for 3D surface on "+fname, 20, 50);
    g.drawString("C colored faces, W wireframe, V vertices, Nodes sz SZ size ", 20, 65);
    g.drawString("< X >  v Y ^  < Z >  < R >  ^ P v  v H ^    W    V    C    N   sz  SZ    Nn",
                 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(xm>430 && xm<460) key = 'n';
      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;}
      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*1.15;} // size smaller
      if(key == 'S') {gsize = gsize*0.85;} // size larger
      if(key == 'n') {nonum = 1-nonum;}
      
      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++)
    //   System.out.println("poly="+priindex[i]+", k="+pripoly[priindex[i]]+
    //                      ", z="+priz[i]);
  } // end order_polys

  public static void main(String[] args)
  {
      new light_dat(args[0]); // e.g. blivet.dat
  }  
} // end class light_dat