// volume_dat2.java  compute the volume of a single closed .dat graphics figure
//          compile with javac -cp . datread.java
//          compile with javac -cp . volume_dat2.java
//          run with java -cp . volume_dat  bull.dat > volume_dat_bull.out
//          sailors algorithm - normal vector must be correct - best triangles

public class volume_dat2
{
  boolean debug = false;
    
  volume_dat2(String filename)
  {
    datread DR; // = new datread(filename); // defines data structures
    int num_points = 0; // set from datread
    int num_polys = 0; // set from datread
    double xmin, xmax, ymin, ymax, zmin, zmax;
    double x1=0.0, x2=0.0, x3=0.0, y1=0.0, y2=0.0, y3=0.0,
	   z1=0.0, z2=0.0, z3=0.0;
    double area=0.0, avgz, vol, dvol, minz, znorm, tarea;
    double earea, evol;
    double a, b, c;
    int i, j, k, npts, iv, pt, pts;
    double x[] = new double[32]; // max 32 points in polygon 
    double y[] = new double[32];
    double z[] = new double[32];
    double ta[] = new double[1024];

    System.out.println("volume_dat2.java reading "+filename);
    DR = new datread(filename); // defines data structures
    num_points = DR.num_points;
    num_polys = DR.num_polys;
    System.out.println("num_points="+num_points+
		       ", num_polys="+num_polys);

    xmin = DR.data_points[0].x;
    xmax = xmin;
    ymin = DR.data_points[0].y;
    ymax = ymin;
    zmin = DR.data_points[0].z;
    zmax = zmin;

    System.out.println("print vertices");
    for(i=0; i<num_points; i++) // vertices
    {
      System.out.println(+(i+1)+"  "+DR.data_points[i].x+
                         "  "+DR.data_points[i].y+"  "+DR.data_points[i].z);
      xmin = Math.min(xmin,DR.data_points[i].x);
      xmax = Math.max(xmax,DR.data_points[i].x);
      ymin = Math.min(ymin,DR.data_points[i].y);
      ymax = Math.max(ymax,DR.data_points[i].y);
      zmin = Math.min(zmin,DR.data_points[i].z);
      zmax = Math.max(zmax,DR.data_points[i].z);
    } // end print vertices

    System.out.println(" ");
    System.out.println("xmin="+xmin+", xmax="+xmax+", ymin="+ymin+
		       ", ymax="+ymax);
    System.out.println("zmin="+zmin+", zmax="+zmax);
    earea = 2.0*((xmax-xmin)*(ymax-ymin)+(xmax-xmin)*(zmax-zmin)+
		 (ymax-ymin)*(zmax-zmin));
    evol = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
    System.out.println("enclosing area="+earea+
		       ", enclosing volume="+evol);
    System.out.println(" ");


    System.out.println("move all vertices positive, 0 to max");
    for(i=0; i<num_points; i++) // vertices
    {
      DR.data_points[i].x -= xmin;
      DR.data_points[i].y -= ymin;
      DR.data_points[i].z -= zmin;
    } // end move all vertices

    area = 0.0;
    k = 0;
    for(i=0; i<num_polys; i++)
    {
      pts = DR.data_polys[k++];
      System.out.print(+pts); // no cr
      for(j=0; j<pts; j++)
      {
	pt = DR.data_polys[k++];
        System.out.print(" "+ pt);
        if(j<3)
	{ 
	  x[j] = DR.data_points[pt-1].x;
	  y[j] = DR.data_points[pt-1].y;
	  z[j] = DR.data_points[pt-1].z;
	} // end j<3
      } // end j
      System.out.println(" "); // point numbers
      ta[i] = tarea(x, y, z);
      area += ta[i];
      if(debug) System.out.println("x0="+x[0]+" y0="+y[0]+" z0="+z[0]);
      if(debug) System.out.println("x1="+x[1]+" y1="+y[1]+" z1="+z[1]);
      if(debug) System.out.println("x2="+x[2]+" y2="+y[2]+" z2="+z[2]);
      if(debug) System.out.println("tarea["+i+"]="+ta[i]);
      if(debug) System.out.println(" ");
    } // end i  print polys
    System.out.println("area of triangles ="+area);
	
  
    // uncomment if data is clockwise looking from outside 
    // datccw(data_points[], num_points, data_verts[], num_polys); 
    // if(debug) datprint(data_points[], num_points, data_verts[], num_polys);

    vol = 0.0; // total volume 
    tarea = 0.0; // total area  
    k = 0; // npts, then npts node numbers, etc  
    for(i=0; i<num_polys; i++)
    {
      npts = DR.data_polys[k];
      k++;
      for(j=0; j<npts; j++)
      {
        iv = DR.data_polys[k];
        k++;
        x[j] = DR.data_points[iv-1].x;
        y[j] = DR.data_points[iv-1].y;
        z[j] = DR.data_points[iv-1].z;
      }
      area = tarea(x, y, z);
      avgz = 0.0;
      for(j=0; j<npts; j++)
        avgz = avgz + z[j]; // all z's are positive 
      avgz = avgz/(double)npts; // accurate for triangles, approx for other 

      znorm  = znormal(x, y, z);
      tarea = tarea + area;
      dvol = area * znorm * avgz;
      vol = vol + dvol;
      if(debug) System.out.println("area="+area+", avgz="+avgz+
		", znorm="+znorm+ ", dvol="+dvol+", vol="+vol);
    }
    System.out.println("final total area="+tarea+
		     ", total volume="+Math.abs(vol));
    System.out.println("volume_dat2.java  finished ");
  } // end volume_dat2

  double tarea(double x[], double y[], double z[]) // triangle area
  {
    double ax, ay, az, bx, by, bz, cx, cy, cz;
    double la, lb, lc, p;
    
    ax = x[2] - x[0];
    ay = y[2] - y[0];
    az = z[2] - z[0];
    bx = x[1] - x[0];
    by = y[1] - y[0];
    bz = z[1] - z[0];
    cx = x[2] - x[1];
    cy = y[2] - y[1];
    cz = z[2] - z[1];
    la = Math.sqrt(ax*ax + ay*ay + az*az);
    lb = Math.sqrt(bx*bx + by*by + bz*bz);
    lc = Math.sqrt(cx*cx + cy*cy + cz*cz);
    p = (la+lb+lc)/2.0;
    return Math.sqrt(p*(p-la)*(p-lb)*(p-lc));
  } // end tarea

  double znormal(double x[], double y[], double z[]) // of triangle
  {
    double ax, ay, az, bx, by, bz, nx, ny, nz=0.0, ss;

    ax = x[2] - x[0];
    ay = y[2] - y[0];
    az = z[2] - z[0];
    bx = x[1] - x[0];
    by = y[1] - y[0];
    bz = z[1] - z[0];
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;
    ss = Math.sqrt(nx*nx+ny*ny+nz*nz);
    nx = nx / ss;
    ny = ny / ss;
    nz = nz / ss;
    return nz;
  } // end znormal 

  public static void main(String[] args)
  {
    String filename = args[0];
    new volume_dat2(filename);
  }
} // end class  volume_dat2.java