// light_normal_stl.java  projection may not be best
//                        needs to size and center
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_normal_stl  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 maxpolys;
  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 = true;
  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;
  int vert = 1;
  int colorz = 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;
  String stlfile;

  light_normal_stl(String filename)
  {
    System.out.println("light_normal_stl.java reading "+filename);
    stlfile = 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_normal_stl");
    setSize(win_w,win_h);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    this.addMouseListener (new mousePressHandler());
  } // end light_normal_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.drawString("light_normal_stl for 3D surface on "+fname, 40, 40);
    g.drawString("normal from blue dot on surface, to red dat", 40, 55);
    // g.drawString("stl file  "+stlfile, 50,50);
    g.drawString("xmin="+xmin,  50,70);
    g.drawString("xmax="+xmax,  50,85);
    g.drawString("ymin="+ymin, 250,70);
    g.drawString("ymax="+ymax, 250,85);
    g.drawString("zmin="+zmin, 450,70);
    g.drawString("zmax="+zmax, 450,85);
    g.drawString("< X >  v Y ^  < Z >  < R >  ^ P v  v H ^    W    V    C   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);
    if(wire==1) draw_wire_dat(g);
    if(vert==1) draw_vert_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 = 's';
      if(xm>370 && xm<400) key = 'S';
      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;}
      if(key == 'v' || key == 'V') {vert = 1-vert;}
      if(key == 'c' || key == 'C') {colorz = 1-colorz;}
      if(key == 's') {gsize = gsize*1.15;}
      if(key == 'S') {gsize = gsize*0.85;}

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


  void draw_wire_dat(Graphics g)
  {
    int i, j, k, pts, pt, pt1;
    int skipz=0;
    double winxyz[] = new double[3];
    double nxyz[] = new double[3];
    double avgx, avgy, avgz, zmin=0.0, zmax=0.0, zavg=0.0;
    double lenn = size/15.0; // for normal vector
    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);

    System.out.println("in draw_wire_dat  lenn="+lenn);

    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;
    }
    setForeground(Color.black);

    // 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);
      avgx = 0.0;
      avgy = 0.0;
      avgz = 0.0;
      compute_normal(k, nxyz);
      // 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];
        avgx += x1;
        avgy += y1;
        avgz += z1;
        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);
      }
      avgx /= (double)pts;
      avgy /= (double)pts;
      avgz /= (double)pts;;
      if(debug) System.out.println("avgx="+avgx+", avgy="+avgy+", avgz="+avgz);
      x1 = avgx;
      y1 = avgy;
      z1 = avgz;
      project(x1, y1, z1, winxyz);
      x1 = winxyz[0];
      y1 = winxyz[1];
      z1 = winxyz[2];
      x2 = avgx+lenn*nxyz[0];
      y2 = avgy+lenn*nxyz[1];
      z2 = avgz+lenn*nxyz[2];
      project(x2, y2, z2, winxyz);
      x2 = winxyz[0];
      y2 = winxyz[1];
      z2 = winxyz[2];
      // normal line blue to red
      g.setColor(Color.black);
      g.drawLine((int)x1, win_h-(int)y1, (int)x2, win_h-(int)y2); 
      g.setColor(Color.blue);
      g.fillOval((int)x1-2, win_h-(int)y1-2, 4, 4);
      g.setColor(Color.red);
      g.fillOval((int)x2-2, win_h-(int)y2-1, 4, 4);
      g.setColor(Color.black);
    } // end loop on num_polys 
  } // end draw_wire_dat 

  void draw_vert_dat(Graphics g)
  {
    int i, j, k, pts, pt;
    double winxyz[] = new double[3];
    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);
        g.drawString(" "+pt,(int)winxyz[0]-2, win_h-(int)winxyz[1]-6);
      }
    } // end loop on num_polys 
  } // end draw_vert_dat 

  void compute_normal(int ivert, double nxyzout[])
  {
    int j, k, pt;
    double ax=0.0, ay=0.0, az=0.0, bx=0.0, by=0.0, bz=0.0;
    double nx, ny, nz, s;
    double v[][] = new double[3][3];

    k = ivert;
    for(j=0; j<3; j++) // use first three points, assume planar 
    {
      pt = data_polys[k++];
      if(debug) System.out.println("j="+j+", pt="+pt);
      v[j][0] = x[pt-1];
      v[j][1] = y[pt-1];
      v[j][2] = z[pt-1];
    }
    // compute normals, assuming planar, normalize 
    ax = v[2][0] - v[1][0];
    ay = v[2][1] - v[1][1];
    az = v[2][2] - v[1][2];
    bx = v[1][0] - v[0][0];
    by = v[1][1] - v[0][1];
    bz = v[1][2] - v[0][2];
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;
    s = Math.sqrt(nx*nx+ny*ny+nz*nz);
    nx = nx / s;
    ny = ny / s;
    nz = nz / s;
    if(debug)
    {
      System.out.println("\ncompute_normal for ivert="+ivert);
      System.out.println("nx="+nx+", ny="+ny+", nz="+nz); 
    }
    nxyzout[0] = nx;
    nxyzout[1] = ny;
    nxyzout[2] = nz;
  } // end compute normal 

  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;

  
    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];
    System.out.println("\n tmatrix");
    prmat(tmatrix);
    matmul(mmatrix, tmatrix, c);
    matcpy(c, mmatrix);
    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 winxyz[])
  {
    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);
    winxyz[0] = ax;
    winxyz[1] = ay;
    winxyz[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

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