/* volume_dat2.c   compute the volume of a single closed .dat graphics figure*/
/*                 compile with datread.c                                    */
/*                 volume_dat  bull.dat > volume_dat_bull.out                */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "datread.h"

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef max
#define max(x,y) ((x)>(y)?(x):(y))
#undef min
#define min(x,y) ((x)<(y)?(x):(y))

float znormal(int n, float x[], float y[], float z[], float *area);

static int num_points;            /* for datread  */
static int num_polys;
static dpts * data_points;
static int *  data_verts;
static float size = 1.0;     /* auto scale */
static char filename[100] = "cube1.dat";
static int status = -1;
static int debug = 0;
static float xmin, xmax, ymin, ymax, zmin, zmax;

int main(int argc, char * argv[])
{
  float area, avgz, vol, dvol, minz, znorm, tarea;
  float offs;
  int i, j, k, npts, iv;
  float x[10], y[10], z[10]; /* max 10 points in polygon */

  if(argc>1)
  {
    strcpy(filename,argv[1]);
  }
  printf("volume_dat2.c reading %s \n", filename);

  status = datread(filename, &data_points, &num_points,
                             &data_verts, &num_polys, &size);
  printf("status=%d, zmax=%g, points=%d, polys=%d \n",
	 status, size, num_points, num_polys);

  xmin = data_points[0].x;
  xmax = xmin;
  ymin = data_points[0].y;
  ymax = ymin;
  zmin = data_points[0].z;
  zmax = zmin;
  
  /* 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 */
  offs = size; /* z's must be positive */
  k = 0; /* npts, then npts node numbers, etc */ 
  for(i=0; i<num_polys; i++)
  {
    npts = data_verts[k];
    k++;
    for(j=0; j<npts; j++)
    {
      iv = data_verts[k];
      k++;
      x[j] = data_points[iv-1].x;
      y[j] = data_points[iv-1].y;
      z[j] = data_points[iv-1].z;
      xmin = min(xmin,x[j]);
      xmax = max(xmax,x[j]);
      ymin = min(ymin,y[j]);
      ymax = max(ymax,y[j]);
      zmin = min(zmin,z[j]);
      zmax = max(zmax,z[j]);
    }
    avgz = 0.0;
    for(j=0; j<npts; j++)
      avgz = avgz + (z[j]+offs); /* z's must be positive */
    avgz = avgz/(float)npts; /* accurate for triangles, approx for other */

    x[npts] = x[0];
    y[npts] = y[0];
    znorm  = znormal(npts, x, y, z, &area);
    tarea = tarea + area;
    dvol = area * znorm * avgz;
    vol = vol + dvol;
    if(debug) printf("area=%g, avgz=%g, znorm=%g, dvol=%g, vol=%g \n",
                     area, avgz, znorm, dvol, vol );
    if(debug) if(avgz>offs && znorm<0.0) printf("upper normal is negative \n");
    if(debug) if(avgz<offs && znorm>0.0) printf("lower normal is positive \n");
  }
  printf("\nxmin=%f, xmax=%f, ymin=%f, ymax=%f \n",
         xmin, xmax, ymin, ymax);
  printf("zmin=%f, zmax=%f \n",
          zmin, zmax);
  printf("enclosing area= %g, enclosing volume= %g \n",
         2.0*((xmax-xmin)*(ymax-ymin)+(xmax-xmin)*(zmax-zmin)+(ymax-ymin)*(zmax-zmin)),
         (xmax-xmin)*(ymax-ymin)*(zmax-zmin));
         
  printf("\nfinal total area = %g, total volume = %g \n\n", tarea, abs(vol));
  return 0;
} /* end main of volume_dat.c */

float znormal(int n, float x[], float y[], float z[], float *area)
{
  float ax, ay, az, bx, by, bz, nx, ny, nz, ss;
  float cx, cy, cz, sss, sstot, ssstot;
  int i;

  sstot = 0.0;
  ssstot = 0.0;
  for(i=2; i<n; i++) /* once for three point triangle */
  {
    ax = x[i] - x[0]; /* |ax ay az| is a direction vector */
    ay = y[i] - y[0];
    az = z[i] - z[0];
    bx = x[i-1] - x[0]; /* |bx by bz| is a direction vector */
    by = y[i-1] - y[0];
    bz = z[i-1] - z[0];
    nx = ay*bz-az*by;   /* |nx ny nz| is a direction vector */
    ny = az*bx-ax*bz;   /* n is vector cross product  a x b */
    nz = ax*by-ay*bx;
    ss = sqrt(nx*nx+ny*ny+nz*nz); /* normalize length to 1.0 */
    if(i==2) /* use only first three points for normal */
    {        /* inaccurate for non planar faces        */
      nx = nx / ss;
      ny = ny / ss;
      nz = nz / ss;
    }
    ss = ss/2.0;
    sstot = sstot + ss;

    /* check */
    cx = x[i] - x[i-1];
    cy = y[i] - y[i-1];
    cz = z[i] - z[i-1];
    ax = sqrt(ax*ax+ay*ay+az*az);
    bx = sqrt(bx*bx+by*by+bz*bz);
    cx = sqrt(cx*cx+cy*cy+cz*cz);
    sss =(ax+bx+cx)/2.0;
    sss = sqrt(sss*(sss-ax)*(sss-bx)*(sss-cx));
    ssstot = ssstot + sss;
  }
  if(abs(ssstot-sstot)>0.1)
     printf("bad area ss=%g, sss=%g \n", sstot, ssstot);
  *area = sstot;
  return nz;
} /* end znormal */