/* volume_dat.c  compute the volume of a single closed .dat graphics figure */
/*               compile with datread.c                                     */

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

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))

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;

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_dat.c reading %s \n", filename);

  status = datread(filename, &data_points, &num_points,
                             &data_verts, &num_polys, &size);
  printf("status=%d, size=%g\n", status, size);
  
  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 = 10.0; /* z's must be positive */
  k = 0;
  for(i=0; i<num_polys; i++)
  {
    npts = data_verts[k];
    k++;
    for(j=0; j<npts; j++)
    {
      iv = data_verts[k];
      x[j] = data_points[iv-1].x;
      y[j] = data_points[iv-1].y;
      z[j] = data_points[iv-1].z;
      k++;
    }
    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;
    if(debug) printf("area=%g, avgz=%g, znorm=%g, dvol=%g \n",
                     area, avgz, znorm, dvol );
    vol = vol + dvol;
    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("\nfinal total area = %g, total volume = %g \n\n", tarea, 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];
    ay = y[i] - y[0];
    az = z[i] - z[0];
    bx = x[i-1] - x[0];
    by = y[i-1] - y[0];
    bz = z[i-1] - z[0];
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;
    ss = sqrt(nx*nx+ny*ny+nz*nz);
    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.0001)
     printf("bad area ss=%g, sss=%g \n", sstot, ssstot);
  *area = sstot;
  return nz;
} /* end znormal */