// Copyright 2001 Ken Perlin

// additions for VRML conversion- jeremi : 
// public Geometry cylinder(int k, double radius, double height)
// public Geometry cube( double width, double height, double depth)
// public Geometry ball(int n, double radius)
// public void computeSurfaceNormals()

// additions for Picking - C. Robbins
// public void calculateBoundingBox()
//

package render;

import java.awt.*;

// A VERY SIMPLE 3D RENDERER BUILT IN JAVA 1.0 - KEN PERLIN

public class Geometry
{
	
   private String notice = "Copyright 2001 Ken Perlin. All rights reserved.";

   public Material material;
   public Geometry child[], refGeometry;
   public double[] color = { 1,1,1, 0,0,0,1 };
   public double[][] matrix = new double[4][4];
   public double[][] globalMatrix = new double[4][4];
   public double[][] refMatrix = new double[4][4];
   public int[][] faces;
   public double[][] vertices, refVertices;
   public double transparency = 0, glow[] = {0,0,0};
   public int meshRowSize = -1;
   public boolean computedMeshNormals = false;
   public static double[][] mat = new double[4][4];
   public boolean modified = false;
   public double[][] boundingBox;              // Two points for bounding box,
                                               //    boundingBox[0] = lower left corner
                                               //    boundingBox[1] = upper right corner
   public Geometry boundingBoxGeometry;        // Cube geometry representing bounding box

   public boolean boundingBoxVisible = false;  // If true geometry will include a transparent
                                               //  bounding box cube when rendered.
   public boolean isBoundingBox = false;       // If true then this geometry is actually the
                                               //  geometry represenation of a boundingBox
	                                                                           
   public static double pullWeight = 1;
   public static int pullMask = 0xffffffff;

   public Geometry() {
      Matrix.identity(matrix);
      Matrix.identity(globalMatrix);
      Matrix.identity(m());
   }

   public void setRefGeometry(Geometry s) { refGeometry = s; }

   /**
      Add a new child shape to a shape.
   */
   public Geometry add() {
      Geometry s = add(new Geometry());
      s.material = material;
      return s;
   }
   /**
      Add an existing shape, as a child to this shape.
   */
   public Geometry add(Geometry s) {
      if (child == null)
         child = new Geometry[16];
      else if (child[child.length-1] != null) {
         Geometry c[] = child;
         child = new Geometry[2*c.length];
         for (int i = 0 ; i < c.length ; i++)
            child[i] = c[i];
      }

      for (int i = 0 ; i < child.length ; i++)
         if (child[i] == null) {
            child[i] = s;
            break;
         }

     return s;
   }
   /**
      Delete a child of a shape.
   */
   public Geometry delete(Geometry s) {
      if (child != null)
         for (int i = 0 ; i < child.length ; i++)
            if (child[i] == s) {
	       delete(i);
	       break;
	    }
      return this;
   }
   /**
      Delete the nth child of a shape.
   */
   public Geometry delete(int n) {
      if (child != null && n >= 0 && n < child.length) {
         for ( ; n < child.length-1 && child[n+1] != null ; n++)
            child[n] = child[n+1];
         child[n] = null;
      }
      return this;
   }
   /**
      Set the Material for this shape to be Material m.
   */
   public Geometry setMaterial(Material m) {
      material = m;
      if (child != null)
         for (int i = 0 ; i < child.length ; i++)
            if (child[i] != null)
               child[i].setMaterial(m);
      return this;
   }
   /**
      Set the transformation matrix for this shape.
   */
   public void setMatrix(double m[][]) {
      Matrix.copy(m, matrix);
   }
   /**
   */
   public double[][] copyVertices(double src[][]) {
      return copyVertices(src, new double[src.length][6]);
   }

   public double[][] copyVertices(double src[][], double dst[][]) {
      for (int i = 0 ; i < src.length ; i++)
      for (int j = 0 ; j < 6 ; j++)
         dst[i][j] = src[i][j];
      return dst;
   }

   public void setVertex(int i,double x , double y , double z ,
                               double nx, double ny, double nz) {
      setVertex(vertices[i], x,y,z, nx,ny,nz);
   }

   public static void setVertex(double v[],double x , double y , double z ,
                                           double nx, double ny, double nz) {
      v[0] = x;
      v[1] = y;
      v[2] = z;
      v[3] = nx;
      v[4] = ny;
      v[5] = nz;
   }

// CUBE

   public Geometry cube() {
      faces = cubeFaces;
      vertices = cubeVertices;
      return this;
   }

   private static int[][] cubeFaces = {
      { 0, 1, 2, 3}, { 4, 5, 6, 7},
      { 8, 9,10,11}, {12,13,14,15},
      {16,17,18,19}, {20,21,22,23},
   };

   private static double N=-1, P=1;
   private static double[][] cubeVertices = {

      {N,N,N, N,0,0},{N,N,P, N,0,0},{N,P,P, N,0,0},{N,P,N, N,0,0},
      {P,N,N, P,0,0},{P,P,N, P,0,0},{P,P,P, P,0,0},{P,N,P, P,0,0},

      {N,N,N, 0,N,0},{P,N,N, 0,N,0},{P,N,P, 0,N,0},{N,N,P, 0,N,0},
      {N,P,N, 0,P,0},{N,P,P, 0,P,0},{P,P,P, 0,P,0},{P,P,N, 0,P,0},

      {N,N,N, 0,0,N},{N,P,N, 0,0,N},{P,P,N, 0,0,N},{P,N,N, 0,0,N},
      {N,N,P, 0,0,P},{P,N,P, 0,0,P},{P,P,P, 0,0,P},{N,P,P, 0,0,P},
   };


// CUBE FOR VRML (includes size components)

    public Geometry cube( double width, double height, double depth) {
	
	add().cube();
	
	Matrix.identity(m);
	Matrix.scale (m,  width/2.0, height/2.0, depth/2.0 );
	child[0].setMatrix(m);
	//scale ( width/2.0, height/2.0, depth/2.0 );
	//transform(child[0]);
	return this;
    }

// BEZELED CUBE

   public Geometry bezeledCube(double r) {
      faces = bezeledCubeFaces;
      vertices = copyVertices(bezeledCubeVertices);
      for (int i = 0 ; i < vertices.length ; i++)
      for (int j = 0 ; j < 3 ; j++) {
         if (vertices[i][j] == n) vertices[i][j] = r-1;
         if (vertices[i][j] == p) vertices[i][j] = 1-r;
      }
      return this;
   }

   private static int[][] bezeledCubeFaces = {
      { 0, 1, 2, 3}, { 4, 5, 6, 7},
      { 8, 9,10,11}, {12,13,14,15},
      {16,17,18,19}, {20,21,22,23},

      { 8,11, 1, 0}, {20,23, 2, 1},
      {13,12, 3, 2}, {17,16, 0, 3},

      {16,19, 9, 8}, {11,10,21,20},
      {23,22,14,13}, {12,15,18,17},

      { 4, 7,10, 9}, { 7, 6,22,21},
      { 6, 5,15,14}, { 5, 4,19,18},

      {16, 8, 0}, {11,20, 1}, {23,13, 2}, {12,17, 3},
      { 9,19, 4}, {21,10, 7}, {14,22, 6}, {18,15, 5},
   };

   private static double n=-.9,p=.9;
   private static double[][] bezeledCubeVertices = {

      {N,n,n, N,0,0},{N,n,p, N,0,0},{N,p,p, N,0,0},{N,p,n, N,0,0},
      {P,n,n, P,0,0},{P,p,n, P,0,0},{P,p,p, P,0,0},{P,n,p, P,0,0},

      {n,N,n, 0,N,0},{p,N,n, 0,N,0},{p,N,p, 0,N,0},{n,N,p, 0,N,0},
      {n,P,n, 0,P,0},{n,P,p, 0,P,0},{p,P,p, 0,P,0},{p,P,n, 0,P,0},

      {n,n,N, 0,0,N},{n,p,N, 0,0,N},{p,p,N, 0,0,N},{p,n,N, 0,0,N},
      {n,n,P, 0,0,P},{p,n,P, 0,0,P},{p,p,P, 0,0,P},{n,p,P, 0,0,P},
   };

// POLYGON

   public Geometry polygon(double X[], double Y[]) {

      int k = X.length;
      faces = new int[1][k];
      vertices = new double[k][6];

      for (int i = 0 ; i < k ; i++) {
         faces[0][i] = i;
         setVertex(i, X[i],Y[i],0, 0,0,1);
      } 

      return this;
   }


// DISK

   public Geometry disk(int k) {

      faces = new int[k][3];
      vertices = new double[k+1][6];

      for (int i = 0 ; i < k ; i++) {
	      faces[i][0] = i;
	      faces[i][1] = (i+1) % k;
	      faces[i][2] = k;
         double theta = 2 * Math.PI * i / k;
         setVertex(i, Math.cos(theta),Math.sin(theta),0, 0,0,1);
      }
      setVertex(k, 0,0,0, 0,0,1);
/*
      faces = new int[1][k];
      vertices = new double[k][6];

      for (int i = 0 ; i < k ; i++) {
         faces[0][i] = i;
         double theta = 2 * Math.PI * i / k;
         setVertex(i, Math.cos(theta),Math.sin(theta),0, 0,0,1);
      } 
*/
      return this;
   }

   private double m[][] = new double[4][4];

// CYLINDER

   public Geometry cylinder(int k) {

      child = null;

      add().tube(k);
      add().disk(k);
      add().disk(k);

      Matrix.identity(m);
      Matrix.scale(m, 1,-1,-1);
      Matrix.translate(m, 0,0,1);
      child[1].setMatrix(m);

      Matrix.identity(m);
      Matrix.translate(m, 0,0,1);
      child[2].setMatrix(m);

      return this;
   }

// CYLINDER FOR VRML DEFAULTS with size accomodation

    public Geometry cylinder(int k, double radius, double height) {

      child = null;

      add().tube(k);
      add().disk(k);
      add().disk(k);

      Matrix.identity(m);
      Matrix.scale(m, 1,-1,-1);
      Matrix.rotateX (m,  Math.PI/2.0 );
      Matrix.translate(m, 0,0,height/2.0);
      Matrix.scale(m, radius, radius, radius);
      child[1].setMatrix(m);

      Matrix.identity(m);
      Matrix.rotateX (m,  Math.PI/2.0 );
      Matrix.translate(m, 0,0,height/2.0);
      Matrix.scale(m, radius, radius, radius);
      child[2].setMatrix(m);

      rotateX( Math.PI / 2.0 );
      scale( (2.0 * radius)/2.0, (2.0 * radius) / 2.0, height/2.0 );
      transform( child[0] );
      
      return this;
    }

// PILL (TUBE WITH ROUNDED ENDS)

   public Geometry pill(int k, double len, double bulge) {
      return pill(k, len, bulge, 1);
   }

// TAPERED PILL (TAPERED TUBE WITH ROUNDED ENDS)

   public Geometry pill(int k, double len, double bulge, double taper) {

      child = null;

      add().tube(k, taper);
      add().globe(k,k/2,0,1,0,.5);
      add().globe(k,k/2,0,1,.5,1);

      Matrix.identity(m);
      Matrix.scale(m, 1,1,len);
      child[0].setMatrix(m);

      Matrix.identity(m);
      Matrix.translate(m, 0,0,-len);
      Matrix.scale(m, 1,1,bulge);
      child[1].setMatrix(m);

      Matrix.identity(m);
      Matrix.translate(m, 0,0,len);
      Matrix.scale(m, taper,taper,bulge);
      child[2].setMatrix(m);

      return this;
   }

// BALL

   public Geometry ball(int n) {

      child = null;

      for (int i = 0 ; i < 6 ; i++) {
         Geometry s = add().newBallFace(n);
         switch (i) {
         case 1: Matrix.rotateX(s.matrix, Math.PI/2); break;
         case 2: Matrix.rotateX(s.matrix, Math.PI  ); break;
         case 3: Matrix.rotateX(s.matrix,-Math.PI/2); break;
         case 4: Matrix.rotateY(s.matrix, Math.PI/2); break;
         case 5: Matrix.rotateY(s.matrix,-Math.PI/2); break;
         }
      }
      return this;
   }

    // BALL FOR VRML (with radius component)
    public Geometry ball(int n, double radius) {

      child = null;
      add().ball(n);
      scale( radius, radius, radius );
      transform ( child[0] );
      return this;
   }

   private Geometry newBallFace(int n) {
      newRectangularMesh(n, n);
      int N = 0;
      for (int j = 0 ; j <= n ; j++)
      for (int i = 0 ; i <= n ; i++) {
         double x = Math.tan(Math.PI/4 * (j-n/2) / (n/2));
         double y = Math.tan(Math.PI/4 * (i-n/2) / (n/2));
         double r = Math.sqrt(x*x + y*y + 1);
         x /= r;
         y /= r;
         double z = -1 / r;
         setVertex(N++, x, y, z, x, y, z);
      }
      computedMeshNormals = true;
      return this;
   }

// GLOBE (LONGITUDE/LATITUDE SPHERE)

   public Geometry globe(int m,int n) {
      return globe(m, n, 0, 1, 0, 1);
   }

// PARAMETRIC SUBSECTION OF A GLOBE

   public Geometry globe(int m,int n,double uLo,double uHi,double vLo,double vHi) {
      newRectangularMesh(m, n);
      int N = 0;
      for (int j = 0 ; j <= n ; j++)
      for (int i = 0 ; i <= m ; i++) {
         double u = uLo + i * (uHi - uLo) / m;
         double v = vLo + j * (vHi - vLo) / n;
         double theta = 2 * u * Math.PI;
         double phi = (v-.5) * Math.PI;
         double x = Math.cos(phi) * Math.cos(theta);
         double y = Math.cos(phi) * Math.sin(theta);
         double z = Math.sin(phi);
         setVertex(N++, x,y,z, x,y,z);
      }
      computedMeshNormals = true;
      return this;
   }

// TUBE

   public Geometry tube(int k) {
      double path[][] = {{0,0,-1, 1,0,0},{0,0,0, 1,0,0},{0,0,1, 1,0,0}};
      return wire(k, path, 1);
   }

// EXTRUDE A WIRE ALONG A PATH

   public Geometry wire(int m, int n, double key[][], double r) {
      return wire(m, makePath(n, key), r);
   }

   public Geometry wire(int m, double path[][], double r) {
      return extrusion(makeCircle(m, r), path);
   }

// EXTRUDE A PATH ALONG ANOTHER PATH

   public Geometry extrusion(double O[][], double P[][]) {

      int m = O.length-1;
      int n = P.length-1;

      newRectangularMesh(m, n);

      double U[] = new double[3];
      double V[] = new double[3];
      double W[] = new double[3];

      boolean loop = same(P[0],P[n]);

      int N = 0;
      for (int j = 0 ; j <= n ; j++) {
	 for (int k = 0 ; k < 3 ; k++) {
            U[k] = P[j][k+3];
	    W[k] = j == n ? ( loop ? P[1][k]-P[0][k]
                                   : P[n][k]-P[n-1][k] )
                          : P[j+1][k]-P[j][k];
         }
         double radius = Vec.norm(U);
         computeCrossVectors(W, U, V);
         for (int i = 0 ; i <= m ; i++) {
            double x = O[i][0];
            double y = O[i][1];
            double z = O[i][2];
	    for (int k = 0 ; k < 3 ; k++)
	       vertices[N][k] = P[j][k] + radius * (x*U[k] - y*V[k] + z*W[k]);
            N++;
         }
         double ux = U[0],uy = U[1],uz = U[2];
      }
      if (loop)
         for (int i = 0 ; i <= m ; i++)
            for (int k = 0 ; k < 3 ; k++)
               vertices[indx(m,n,i,n)][k] = vertices[indx(m,n,i,0)][k];

      return this;
   }

   public void computeCrossVectors(double W[], double U[], double V[]) {
      Vec.normalize(W);
      Vec.normalize(U);
      Vec.cross(W,U,V);
      Vec.normalize(V);
      Vec.cross(V,W,U);
      Vec.normalize(U);
   }

   private int _m,_n;
   private double[] V(int i,int j) {
      return vertices[indx(_m,_n,i,j)];
   }

    public void computeSurfaceNormals() {
	// 
	// needs to be inserted into the renderer
	// used for normals computation of indexedFaceSurfaces
	//
	// jeremi july2002
	//
	
	double faceNormals[][] = new double [ faces.length ][3] ;
	double A[] = new double[3];
	double B[] = new double[3];
	
	double vertNormals[][] = new double [ vertices.length ][4] ; 
	//each entry contains accumulated values + count of normals to compute avg. 

	// first compute normals of faces.
	for (int i=0; i< faces.length; i++) {
	    //for each face
	    
	    for (int k=0; k<3; k++) {
		//for each dimension
		
		A[k] = vertices[ faces[i][0] ][ k ] - vertices[ faces[i][1] ][ k ];
		B[k] = vertices[ faces[i][1] ][ k ] - vertices[ faces[i][2] ][ k ];
	    } 
	    Vec.cross( A, B, faceNormals[i] );
	    //Vec.normalize( faceNormals[i] );
	}

	for ( int i=0; i < vertNormals.length; i++ )
	    for ( int k=0; k<4; k++)
		vertNormals[i][k] = 0;

	for (int i=0; i< faces.length; i++) {
	    for (int j=0; j<faces[i].length; j++) {
		for (int k=0; k<3; k++) {
		    vertNormals[ faces[i][j] ][k] += faceNormals[i][k];
		}
		vertNormals[ faces[i][j] ][ 3 ] += 1;
	    }
	}

	for (int i=0; i<vertNormals.length; i++)
	    if ( vertNormals[i][3] != 0) {
		//if normals were contributed .. compute avg
		vertices[i][3] = -vertNormals[i][0] / vertNormals[i][3];
		vertices[i][4] = -vertNormals[i][1] / vertNormals[i][3];
		vertices[i][5] = -vertNormals[i][2] / vertNormals[i][3];
	       
	    }

	computedMeshNormals = true;
    }

    public void computeMeshNormals() {
      if (computedMeshNormals)
         return;

      if (meshRowSize < 0)
         return;

      int m = meshRowSize, n = vertices.length / (m+1) - 1;

      _m = m;
      _n = n;

      double nn[] = new double[3];
      double A[] = new double[3], B[] = new double[3];
      boolean loopI = same(vertices[indx(m,n,0,0)],vertices[indx(m,n,m,0)]);
      boolean loopJ = same(vertices[indx(m,n,0,0)],vertices[indx(m,n,0,n)]);
      int N = 0;
      for (int j = 0 ; j <= n ; j++)
      for (int i = 0 ; i <= m ; i++) {

         if (i == 0)
          for (int k = 0 ; k < 3 ; k++)
           A[k] = 1.5*(V(i+1,j)[k]-V(i,j)[k]) - .5*(V(i+2,j)[k]-V(i+1,j)[k]);
         else if (i == m)
          for (int k = 0 ; k < 3 ; k++)
           A[k] = 1.5*(V(i,j)[k]-V(i-1,j)[k]) - .5*(V(i-1,j)[k]-V(i-2,j)[k]);
         else
          for (int k = 0 ; k < 3 ; k++)
           A[k] = .5 * ( V(i+1,j)[k] - V(i-1,j)[k] );

         if (j == 0)
          for (int k = 0 ; k < 3 ; k++)
           B[k] = 1.5*(V(i,j+1)[k]-V(i,j)[k]) - .5*(V(i,j+2)[k]-V(i,j+1)[k]);
         else if (j == n)
          for (int k = 0 ; k < 3 ; k++)
           B[k] = 1.5*(V(i,j)[k]-V(i,j-1)[k]) - .5*(V(i,j-1)[k]-V(i,j-2)[k]);
         else
          for (int k = 0 ; k < 3 ; k++)
           B[k] = .5 * ( V(i,j+1)[k] - V(i,j-1)[k] );
/*
         int i0 = loopI ? (m+i-1)%m : Math.max(0,i-1);
         int i1 = loopI ? (m+i+1)%m : Math.min(m,i+1);
         int j0 = loopJ ? (n+j-1)%n : Math.max(0,j-1);
         int j1 = loopJ ? (n+j+1)%n : Math.min(n,j+1);
         double a[] = vertices[indx(m,n,i0,j0)];
         double b[] = vertices[indx(m,n,i1,j0)];
         double c[] = vertices[indx(m,n,i0,j1)];
         double d[] = vertices[indx(m,n,i1,j1)];
         for (int k = 0 ; k < 3 ; k++) {
            A[k] = d[k] - a[k];
            B[k] = c[k] - b[k];
         }
*/
         Vec.cross(B, A, nn);
         Vec.normalize(nn);
	 for (int k = 0 ; k < 3 ; k++)
            vertices[N][k+3] = nn[k];
         N++;
      }
      computedMeshNormals = true;
   }

   private int indx(int m, int n, int i, int j) {
      return j * (m+1) + i;
   }

   public Geometry tube(int m, double taper) {
      double T[] = {-1,1};
      double C[] = {1,taper};
      return latheGen(m, T, C, false);
   }

// LATHE (OBJECT FORMED ON A LATHE)

   /**
      Lathe:

         Extrude a curve of varying radius around a circle.

      Arguments:

         m is the number of steps around the circle
         n is the number of steps to sample along the path

         Z and R specify keypoints along the curve such that:

            Z is an array containing successive z coordinates
            R is an array containing successive radii

      Algorithm:

         The path is first expanded out to a smooth curve,
         so that it is sampled evenly in n steps.  The curve
         preserves local maxima and minima in radius (eg:
	 if a keypoint has a greater radius than both
	 its neighbors, then that point will be a local
	 maximum in the resulting curve).

         Then this expanded path is extruded around a circle.
   */
   public Geometry lathe(int m, int n, double Z[], double R[]) {

      double T[] = new double[n+1];
      double C[] = new double[n+1];
      makeCurve(Z, R, T, C);
      return latheGen(m, T, C, true);
   }

   public Geometry latheGen(int m, double T[], double C[], boolean round) {
      
      int n = T.length-1;

      newRectangularMesh(m, n);

      for (int i = 0 ; i <= n ; i++)
      for (int j = 0 ; j <= m ; j++) {
         double theta = 2 * Math.PI * j / m;
         double x = Math.cos(theta);
         double y = Math.sin(theta);
         double z = T[i];
         double r = C[i];
         double sign = T[0] < T[n] ? 1 : -1;
         double dr;
         if (round)
            dr = i==0 ? -sign : i==n ? sign : (C[i+1] - C[i-1])/(2*r);
         else
            dr = (C[1]-C[0])/r;
         double nn[] = {r*x,r*y,dr};
         Vec.normalize(nn);
         setVertex(i * (m+1) + j, r*x,r*y,z, nn[0],nn[1],nn[2]);
      }
      computedMeshNormals = true;

      return this;
   }

   private void newRectangularMesh(int m, int n) {

      meshRowSize = m;

      faces = new int[m*n][4];

      for (int k = 0 ; k < n ; k++)
      for (int j = 0 ; j < m ; j++) {
         int f = k *  m    + j;
         int v = k * (m+1) + j;
         faces[f][0] = v;
         faces[f][1] = v + 1;
         faces[f][2] = v + m+1 + 1;
         faces[f][3] = v + m+1;
      }

      vertices = new double[(m+1)*(n+1)][6];

      computedMeshNormals = false;
   }

   public static double[][] makePath(int n, double key[][]) {

      double P[][] = new double[n+1][6];

      for (int i = 0 ; i <= n ; i++) {
	 double t = i / (n+.001);
	 double f = t * (key.length-1);
	 int k = (int)f;
	 for (int j = 0 ; j < 3 ; j++)
	    P[i][j] = hermite(0,1, key[k][  j], key[k+1][  j],
                                   key[k][3+j], key[k+1][3+j], f % 1.0);
      }
      return P;
   }

   public static double[][] makeCircle(int n, double radius) {
      double P[][] = new double[n+1][6];
      for (int i = 0 ; i <= n ; i++) {
         double theta = 2 * Math.PI * i / n;
         double cos = Math.cos(theta);
         double sin = Math.sin(theta);

         P[i][0] = radius * cos; // LOCATION
         P[i][1] = radius * sin;
         P[i][2] = 0;

         P[i][3] = cos;          // NORMAL DIRECTION
         P[i][4] = sin;
         P[i][5] = 0;
      }
      return P;
   }

   public static double[][] makeRectangle(double w, double h) {
      double P[][] = new double[8][6];

      setVertex(P[0], -w,-h,0, -1,0,0);
      setVertex(P[1], -w, h,0, -1,0,0);

      setVertex(P[2], -w, h,0, 0,1,0);
      setVertex(P[3],  w, h,0, 0,1,0);

      setVertex(P[4],  w, h,0, 1,0,0);
      setVertex(P[5],  w,-h,0, 1,0,0);

      setVertex(P[6],  w,-h,0, 0,-1,0);
      setVertex(P[7], -w,-h,0, 0,-1,0);

      return P;
   }

   public static void makeCurve(double X[],double Y[], double T[],double C[]) {
      double S[] = new double[X.length]; // SLOPE
      int n = X.length;

      for (int i = 1 ; i < n-1 ; i++)
         S[i] = (Y[i]>=Y[i-1]) == (Y[i]>=Y[i+1])
              ? 0
              : ( (Y[i+1] - Y[i  ]) * (X[i  ] - X[i-1])+
                  (Y[i  ] - Y[i-1]) * (X[i+1] - X[i  ]))
                /((X[i+1] - X[i-1]) * (X[i+1] - X[i-1]));

      S[ 0 ] = 2 * (Y[ 1 ]-Y[ 0 ]) / (X[ 1 ]-X[ 0 ]) - S[ 1 ];
      S[n-1] = 2 * (Y[n-1]-Y[n-2]) / (X[n-1]-X[n-2]) - S[n-2];

      int k = C.length;
      for (int j = 0 ; j < k ; j++) {
         double t = j / (k-.99);
         double x = X[0] + t * (X[n-1] - X[0]);
         int i = 0;
         for ( ; i < n-2 ; i++)
            if (x >= X[i] && x < X[i+1])
               break;
         T[j] = x;
         C[j] = hermite(X[i], X[i+1], Y[i], Y[i+1], S[i], S[i+1], x);
      }
   }

   private static double hermite(double x0,double x1,double y0,double y1,
                          double s0,double s1,double x) {
      double t = (x - x0) / (x1 - x0);
      double s = 1-t;
      return y0 * s*s*(3 - 2*s)
           + s0 * (1-s)*s*s
           - s1 * (1-t)*t*t
           + y1 * t*t*(3 - 2*t);
   }

   private boolean same(double a[], double b[]) {
      return Math.abs(a[0]-b[0])<.01 &&
             Math.abs(a[1]-b[1])<.01 &&
             Math.abs(a[2]-b[2])<.01 ;
   }

   public void superquadric(double p, double h) {
      superquadric(globalMatrix, p, h);
   }

   public void superquadric(double mat[][], double p, double h) {

      if (child != null) {
         for (int i = 0 ; i < child.length ; i++)
            if (child[i] != null) {
               double tmp[][] = new double[4][4];
               Matrix.copy(mat, tmp);
               Matrix.preMultiply(tmp, child[i].matrix);
               child[i].superquadric(tmp, p, h);
            }
         return;
      }

      double v[][] = vertices, x, y, z, s;
      double w[] = new double[3];
      double inv[][] = new double[4][4];
      Matrix.invert(mat, inv);
      for (int k = 0 ; k < v.length ; k++) {
         transform(v[k], mat, w);
         x = Math.abs(w[0]);
         y = Math.abs(w[1]);
         z = Math.abs(w[2]);
         double t = Math.pow(Math.pow(x,p)+Math.pow(y,p)+Math.pow(z,p),1/p);
         w[0] /= t;
         w[1] /= t;
         w[2] /= t;

         if (h > 0) {
            t = Math.pow(w[0]*w[0] + w[1]*w[1] + w[2]*w[2], h);
            w[0] *= t;
            w[1] *= t;
            w[2] *= t;
         }

         transform(w, inv, v[k]);
      }
      computedMeshNormals = false;
   }

   public void addNoise(double freq, double ampl) {
      addNoise(globalMatrix, freq, ampl);
   }

   public void addNoise(double mat[][], double freq, double ampl) {

      if (child != null) {
         for (int i = 0 ; i < child.length ; i++)
            if (child[i] != null) {
               double tmp[][] = new double[4][4];
               Matrix.copy(mat, tmp);
               Matrix.preMultiply(tmp, child[i].matrix);
               child[i].addNoise(tmp, freq, ampl);
            }
         return;
      }

      double v[][] = vertices, x, y, z, s;
      double w[] = new double[3];
      double vn[] = new double[3];
      double wn[] = new double[3];
      double inv[][] = new double[4][4];
      double matn[][] = new double[4][4];
      Matrix.copy(mat, matn);
      matn[0][3] = matn[1][3] = matn[2][3] = 0;
      Matrix.invert(mat, inv);
      for (int k = 0 ; k < v.length ; k++) {
         transform(v[k], mat, w);
         for (int j = 0 ; j < 3 ; j++)
            vn[j] = v[k][3+j];
         transform(vn, matn, wn);
         x = freq*w[0];
         y = freq*w[1];
         z = freq*w[2] + 100;
         double t = ampl * Noise.noise(x,y,z);
         for (int j = 0 ; j < 3 ; j++)
            w[j] += t * wn[j];
         transform(w, inv, v[k]);
      }
      computedMeshNormals = false;
   }

   public void sew(Geometry s, int a, int b) {
      sew(s, a, b, null);
   }

   public void sew(Geometry s, int a, int b, double m[][]) {
      int am = meshRowSize, an = vertices.length / (am+1) - 1;
      int bm = s.meshRowSize, bn = s.vertices.length / (bm+1) - 1;
      int j1 = a==0 ? 0 : an*(am+1);
      int j2 = b==0 ? 0 : bn*(bm+1);
      if (am >= bm)
         for (int i = 0 ; i <= am ; i++)
            sewVertex(vertices[j1 + i], m, s.vertices[j2 + i*bm/am]);
      else
         for (int i = 0 ; i <= bm ; i++)
            sewVertex(vertices[j1 + i*am/bm], m, s.vertices[j2 + i]);
   }

   private void sewVertex(double v1[], double m[][], double v2[]) {
      if (m != null)
         transformVertex(v1, m, v2);
      else
         copyVertex(v1, v2);
   }

   public void copyVertex(double src[], double dst[]) {
      for (int j = 0 ; j < dst.length ; j++)
         dst[j] = src[j];
   }

   void transformVertex(double src[], double m[][], double dst[]) {
       for (int i = 0 ; i < 3 ; i++)
          dst[i] = src[0]*m[i][0] +
                   src[1]*m[i][1] +
                   src[2]*m[i][2] + m[i][3];
       dst[3] = src[3];
       dst[4] = src[4];
       dst[5] = src[5];
   }

   public int pull(double m[][], double x0,double x1,double x2,
                                  double y0,double y1,double y2,
                                  double z0,double z1,double z2) {
      double tmp1[][] = new double[4][4];
      Matrix.identity(tmp1);

      double tmp2[][] = new double[4][4];
      Matrix.invert(globalMatrix, tmp2);
      Matrix.preMultiply(tmp2, m);

      return pull(tmp1, tmp2, x0,x1,x2, y0,y1,y2, z0,z1,z2);
   }

   final double UNDEFINED = .001234;
   double vCache[][];

   private double A[] = {0,0,0}, B[] = {0,0,0};
   private double inverseOldRelMatrix[][] = new double[4][4];

   public int pull(double oldRelMatrix[][],
                    double newRelMatrix[][], double x0,double x1,double x2,
                                             double y0,double y1,double y2,
                                             double z0,double z1,double z2) {
      int mask = 0;

      if (child != null) {
         double v[] = new double[3];
         for (int i = 0 ; i < child.length ; i++)
            if (child[i] != null) {

               double tmp1[][] = new double[4][4];
               Matrix.copy(oldRelMatrix, tmp1);
               Matrix.preMultiply(tmp1, child[i].matrix);

               double tmp2[][] = new double[4][4];
               Matrix.copy(newRelMatrix, tmp2);
               Matrix.preMultiply(tmp2, child[i].matrix);

               if ( ( pullMask & (1<<i) ) != 0)
                  if (child[i].pull(tmp1,tmp2,x0,x1,x2,y0,y1,y2,z0,z1,z2)!=0)
                     mask |= (1 << i);
            }
         return mask;
      }

      if (refVertices == null) {
         refVertices = copyVertices(vertices);
         vCache = new double[vertices.length][3];
      }
      if (!modified) {
         copyVertices(refVertices, vertices);
         for (int i = 0 ; i < vCache.length ; i++)
            vCache[i][0] = UNDEFINED;
         modified = true;
      }

      Matrix.invert(oldRelMatrix, inverseOldRelMatrix);

      boolean pulled = false;

      double f, fx, fy, fz;
      for (int i = 0 ; i < vertices.length ; i++) {
         double v[] = vertices[i];
         double w[] = vCache[i];
         if (w[0] == UNDEFINED)
            transform(v, oldRelMatrix, w);

         if ((fx = lump(w[0],x0,x1,x2)) >= 1) continue;
         if ((fy = lump(w[1],y0,y1,y2)) >= 1) continue;
         if ((fz = lump(w[2],z0,z1,z2)) >= 1) continue;

         pulled = true;

         if (fx < 0) fx = 0;
         if (fy < 0) fy = 0;
         if (fz < 0) fz = 0;
         f = fx*fx + fy*fy + fz*fz;

         if (f < 1) {
            if (f == 0)
               transform(v, newRelMatrix, w);
            else {
               transform(v, newRelMatrix, B);
               f = dropoff[(int)(D*f)] * pullWeight;
               w[0] += f * (B[0] - w[0]);
               w[1] += f * (B[1] - w[1]);
               w[2] += f * (B[2] - w[2]);
            }
            transform(w, inverseOldRelMatrix, v);
         }
      }
      if (pulled)
         computedMeshNormals = false;

      return pulled ? 1 : 0;
   }

   private static final int D = 1000;
   private static double dropoff[] = new double[D];
   private static boolean initDropoff = computeDropoff();
   private static boolean computeDropoff() {
      for (int i = 0 ; i < D ; i++) {
         double f = (double)i / D;
         dropoff[i] = .5 + .5 * Math.cos(Math.PI * Math.sqrt(f));
      }
      return true;
   }
   private double lump(double t, double a, double b, double c) {
      if (a==c)
         return 0;
      t = a==b || (t<b && c<b || t>b && c>b) ? (t-b)/(c-b) : (t-b)/(a-b);
      return t>=1 ? 1 : t;
   }
   void transform(double src[], double m[][], double dst[]) {
      dst[0] = m[0][0]*src[0] + m[0][1]*src[1] + m[0][2]*src[2] + m[0][3];
      dst[1] = m[1][0]*src[0] + m[1][1]*src[1] + m[1][2]*src[2] + m[1][3];
      dst[2] = m[2][0]*src[0] + m[2][1]*src[1] + m[2][2]*src[2] + m[2][3];
   }


   // Robbins - PUBLIC METHOD TO RETURN NUMBER OF CHILD GEOMETRIES FOR THIS GEOMETRY
   public int numChildren() {
      
      if(this.child == null) return 0;
      
      int childCount = 0;
      while(childCount < this.child.length && this.child[childCount] != null) {
      		childCount++;
      }
      return childCount;

   }
   
            
   // Robbins - PUBLIC METHOD TO PRINT THIS GEOMETRY'S BOUNDING BOX TO STDOUT         
   public void printBoundingBox() {
      printBoundingBox(this.boundingBox);
   }

   // Robbins - PUBLIC METHOD TO PRINT A GIVEN BOUNDING BOX TO STDOUT   
   public void printBoundingBox(double[][] bb) {      
      if (bb != null) {         
         System.out.print("Bounding Box = { ");
         System.out.print(Vec.toString(bb[0]));
         System.out.print(Vec.toString(bb[1]));
         System.out.println();
         
      }      
   }
         
    // Robbins - PRIVATE METHOD TO CALCULATE THE LOWER LEFT AND UPPER RIGHT 
    //           BOUNDING BOX CORNERS FOR A GIVEN GEOMETRY 
    // Utility function for public method calculateBoundingBox()
    // calculates bounding box lower left and upper right corners from existing vertices
    // then recursively updates these corners from child geometries. 
    // returns lower left and upper right corners in provided locations
    private void calculateCorners(double[] lowerLeft, double[] upperRight,	Geometry g, double[][] curMatrix) {

	final boolean debugLoud = false;

	double[][] identityMatrix = new double[4][4];
	double[] currentVertex = new double[6];
	Matrix.identity(identityMatrix);

	if (g.vertices != null) {

	    if (debugLoud) {
		System.out.println("      " + Matrix.toString(identityMatrix));
		System.out.println("      " + Matrix.toString(curMatrix) + "  ");
	    }

	    for (int counter = 0; counter < g.vertices.length; counter++) {

		if (Matrix.equals(g.matrix, identityMatrix)) {
		    Vec.copy(g.vertices[counter], currentVertex);
		} else {
		    Vec.copy(
			     Vec.matrixVectorMultiply(g.vertices[counter], g.matrix),
			     currentVertex);
		}

		if (debugLoud) System.out.println(Vec.toString(currentVertex));

		if (lowerLeft[0] > currentVertex[0])
		    lowerLeft[0] = currentVertex[0];
		else if (upperRight[0] < currentVertex[0])
		    upperRight[0] = currentVertex[0];

		if (lowerLeft[1] > currentVertex[1])
		    lowerLeft[1] = currentVertex[1];
		else if (upperRight[1] < currentVertex[1])
		    upperRight[1] = currentVertex[1];

		if (lowerLeft[2] > currentVertex[2])
		    lowerLeft[2] = currentVertex[2];
		else if (upperRight[2] < currentVertex[2])
		    upperRight[2] = currentVertex[2];

	    }
	}
	if (g.child != null) {
	    for (int childIndex = 0;
		 childIndex < g.child.length && g.child[childIndex] != null;
		 childIndex++) {
		if (debugLoud)
		    System.out.println(
				       "  processing child # " + childIndex + ":");
		calculateCorners(lowerLeft, upperRight, g.child[childIndex], curMatrix);
	    }
	}
    }
    
    
    // Robbins - PUBLIC METHOD TO CALCULATE BOUNDING BOX 
    // calculates and sets bounding box this geometry
    // returns cube geometry representing Bounding Box
    public double[][] calculateBoundingBox() {

	final boolean debugLoud = false;

	boundingBox = new double[2][3];
      
	double lowerLeft[] = {Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY,	Double.POSITIVE_INFINITY },
	    upperRight[] = {	Double.NEGATIVE_INFINITY,	Double.NEGATIVE_INFINITY,	Double.NEGATIVE_INFINITY };


	    calculateCorners(lowerLeft, upperRight, this, this.matrix);

	    boundingBox[0] = lowerLeft;
	    boundingBox[1] = upperRight;
	    if (debugLoud) printBoundingBox();

	    double width = upperRight[0] - lowerLeft[0];
	    double height = upperRight[1] - lowerLeft[1];
	    double depth = upperRight[2] - lowerLeft[2];

	    if (this.child == null) {
		Geometry moveRootToChild = this.add();
		moveRootToChild.vertices = this.vertices;
		this.vertices = null;
		moveRootToChild.faces = this.faces;
		this.faces = null;
	    }

	    boundingBoxGeometry = add();
	    boundingBoxGeometry.cube(width, height, depth);
	    boundingBoxGeometry.isBoundingBox = true;
	    push();
  	       translate(lowerLeft[0] + (width / 2),	lowerLeft[1] + (height / 2),	lowerLeft[2] + (depth / 2));
	       transform(boundingBoxGeometry);
	    pop();
	    boundingBoxGeometry.setMaterial(new Material().setColor(0.8, 0.8, 0.8, 0.0, 0.0, 0.0, 0.0).setTransparency(1.0));
	    boundingBoxGeometry.transparency = 1;

	    return boundingBox;

    }
	
    // Robbins - PUBLIC METHOD TO MAKE BOUNDING BOX (IN)VISIBLE
    // sets boundingBoxVisible value for this Geometry and adds or deletes boundingBoxGeometry
    // as necessary
    public void setBoundingBoxVisible(boolean visibiltyRequest) {

	Material boundingBoxVisibleMaterial = new Material();  // transparent grey-ish glass
	boundingBoxVisibleMaterial.setColor( 0.8, 0.8, 0.8, 0.0, 0.0, 0.0, 0.0);

	if (this.boundingBox == null) return;
		
	if (visibiltyRequest == true) {
            
	    if (this.boundingBoxVisible == false) {
		boundingBoxVisibleMaterial.setTransparency(0.7);
		boundingBoxGeometry.setMaterial(boundingBoxVisibleMaterial);
		this.boundingBoxVisible = true;
	    }

	} else {
			
	    if (this.boundingBoxVisible == true) {
		boundingBoxVisibleMaterial.setTransparency(1.0);
		boundingBoxGeometry.setMaterial(boundingBoxVisibleMaterial);
		this.boundingBoxVisible = false;
	    }
			
	}		
    }
   
    // Robbins - PUBLIC METHOD TO RETURN NUMBER OF VERTICES IN THIS GEOMETRY
    public int getVertexCount() {
        
	int currentCount = 0;
	int numChildren = this.numChildren();
      
	if (numChildren > 0) {
	    for(int childIndex = 0; childIndex < numChildren; childIndex++) {
		currentCount += child[childIndex].getVertexCount();
	    }
	}
      
	if (vertices != null) {
	    currentCount += this.vertices.length;
	} 
	return currentCount;
      
    }
      
   // PUBLIC METHODS TO LET THE PROGRAMMER MANIPULATE A MATRIX STACK

   public void identity() { Matrix.identity(m()); }
   public double[][] m() { return mstack[top]; }
   public void pop() { top--; }
   public void push() { Matrix.copy(m(),mstack[top+1]); top++; }
   public void rotateX(double t) { Matrix.rotateX(m(), t); }
   public void rotateY(double t) { Matrix.rotateY(m(), t); }
   public void rotateZ(double t) { Matrix.rotateZ(m(), t); }
   public void scale(double x,double y,double z) { Matrix.scale(m(),x,y,z); }
   public void transform() { setMatrix(m()); }
   public void transform(Geometry s) { s.setMatrix(m()); }
   public void translate(double x,double y,double z) { Matrix.translate(m(), x, y, z); }

   private double mstack[][][] = new double[10][4][4]; // THE MATRIX STACK
   private int top = 0;                // MATRIX STACK POINTER
}