/* sphere_div.c  divide a tetrahedron^H^H^H^H^H^H^H^H^H   */
/*                        octahedron, write a .dat file   */
/*               executable link with datread             */
/* now has spike option for split == 10                   */

/* sphere_div out_file.dat 3  # last number is split, max=4 */

/*  tetrahedron   4 vertices     4 faces     6 edges  */
/*  split 1      10 vertices    16 faces    24 edges  */
/*  split 2      34 vertices    64 faces    96 edges  */
/*  split 3     130 vertices   256 faces   384 edges  */
/*  split 4     514 vertices  1024 faces  1536 edges  */

/*  octahedron    6 vertices     8 faces    12 edges  */
/*  split 1      18 vertices    32 faces    48 edges  */
/*  split 2      66 vertices   128 faces   192 edges  */
/*  split 3     258 vertices   512 faces   768 edges  */
/*  split 4    1026 vertices  2048 faces  3072 edges  */

/* num_points = vertices,  num_polys = faces */

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

#undef  abs
#define abs(a) ((a)>0.0?(a):(-(a)))

static int num_points = 0; /* for datread, datwrite and datclean */
static int new_point = 0;  /* new num_points */
static int num_polys = 0;
static int new_polys = 0;      /* new num_polys */
static int num_edges = 0;
static dpts data_points[2050]; /* each has .x .y .z */
static dpts new_points[6150]; /* each has .x .y .z spikes */
static int  data_verts[8200]; /* 4*num_polys is count=3 index */
static int  new_verts[24600];  /* spikes */
static int  link_verts[2050][3]; /* links to three edges */
static int  edge_list[3100][3]; /* old_vert, old_vert, mid_vert, poly */
static float size;
static int status = -1;
static int splits = 1;
static char out_name[128];
static double Pi = 3.141592653589793238462643383279502884197;
static int spike = 0;
static int npts = 0;

typedef double vert_typ[3];
static vert_typ vertices[4]; /* initial tetrahedron */

static void tetrahedron_compute();
static void octahedron_compute();
static void split(int splits);
static double d3len(vert_typ a, vert_typ b);
static void vertex_compute();
static void polygon(int a, int b, int c);
static int insert_edge(int a, int b);
static void split_poly(int i);
static void compute_normal(int ivert, 
                           double *nxout, double *nyout, double *nzout);

int main(int argc, char * argv[])
{
  if(argc>1)
  {
    strcpy(out_name, argv[1]);
    if(argc>2) splits = atoi(argv[2]);
    if(splits>9) { splits=1; spike=1;} // Medieval flail spiked ball
    if(splits<1 || splits>4) splits=1;
    printf("sphere_div about to split %d and write to %s \n",
            splits, argv[1]);
  }
  else
  {
    printf("sphere_div about to split %d and write to sphere_div.dat \n",
           splits);
    strcpy(out_name, "sphere_div.dat");
  }  
  
  split(splits);  /* do the work, set up for datwrite */
  if(spike==1)
  {
    status = datwrite(out_name, new_points, new_point,
                                new_verts,  new_polys);
  printf("sphere_div has written spike to %s \n", out_name);
  }  
  else
  {
    status = datwrite(out_name, data_points, num_points,
                                data_verts,  num_polys);
  printf("sphere_div has written %s \n", out_name);
  }
  return 0;
}

static double d3len(vert_typ a, vert_typ b)
{
  double x1, y1, z1, x2, y2, z2;
  x1 = a[0];  y1 = a[1];  z1 = a[2];
  x2 = b[0];  y2 = b[1];  z2 = b[2];
  return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
}

static void tetrahedron_compute()
{
  double phia, the120, the, r;
  double phiaa  = -19.471220333; /* generator */
  int i;
  
  /* based on unit sphere x*x+y*y+z*z=1           */
  /* x=r*cos(the)*cos(phi)           the 0 to 2Pi */
  /* y=r*sin(the)*cos(phi)    phi -Pi/2 to Pi/2   */
  /* z=r*sin(phi)          radians=Pi*degrees/180 */
  r = 1.0;
  phia = Pi*phiaa/180.0; /* 1 set of three points */
  vertices[0][0] = 0.0;
  vertices[0][1] = 0.0;
  vertices[0][2] = r;
  printf("i=%d, x=%6.3f, y=%6.3f, z=%6.3f \n",
          0, vertices[0][0], vertices[0][1], vertices[0][2]);
  the120 = Pi*120.0/180;
  the = 0.0;
  for(i=1; i<4; i++)
  {
    vertices[i][0]=r*cos(the)*cos(phia);
    vertices[i][1]=r*sin(the)*cos(phia);
    vertices[i][2]=r*sin(phia);
    printf("i=%d, x=%6.3f, y=%6.3f, z=%6.3f \n",
            i, vertices[i][0], vertices[i][1], vertices[i][2]);
    the = the+the120;
  }

  for(i=0; i<4; i++)
  {
    data_points[num_points].x = vertices[i][0];
    data_points[num_points].y = vertices[i][1];
    data_points[num_points].z = vertices[i][2];
    num_points++;
  }

  /* map vertices to 4 faces */
  polygon(0,1,2);
  polygon(0,2,3);
  polygon(0,3,1);
  polygon(1,2,3);
}

static void octahedron_compute()
{
  double phia, the90, the, r;
  int i;
  r = 1.0;
  phia = 0.0; /* 1 set of four points */
  the90 = Pi*90.0/180;
  vertices[0][0]=0.0;
  vertices[0][1]=0.0;
  vertices[0][2]=r;
  printf("i=%d, x=%6.3f, y=%6.3f, z=%6.3f \n",
          0, vertices[0][0], vertices[0][1], vertices[0][2]);
  vertices[5][0]=0.0;
  vertices[5][1]=0.0;
  vertices[5][2]=-r;

  the = 0.0;
  for(i=1; i<5; i++)
  {
    vertices[i][0]=r*cos(the)*cos(phia);
    vertices[i][1]=r*sin(the)*cos(phia);
    vertices[i][2]=r*sin(phia);
    printf("i=%d, x=%6.3f, y=%6.3f, z=%6.3f \n",
            i, vertices[i][0], vertices[i][1], vertices[i][2]);
    the = the+the90;
  }
  printf("i=%d, x=%6.3f, y=%6.3f, z=%6.3f \n",
          5, vertices[5][0], vertices[5][1], vertices[5][2]);

  printf("num_points=%d, num_polys=%d, num_edges=%d \n",
         num_points, num_polys, num_edges);
  num_polys = 0;
  printf("num_points=%d, num_polys=%d, num_edges=%d \n",
         num_points, num_polys, num_edges);
  
  for(i=0; i<6; i++)
  {
    data_points[num_points].x = vertices[i][0];
    data_points[num_points].y = vertices[i][1];
    data_points[num_points].z = vertices[i][2];
    num_points++;
  }

  /* map vertices to 8 faces */
  polygon(0,1,2);
  polygon(0,2,3);
  polygon(0,3,4);
  polygon(0,4,1);
  polygon(5,2,1);
  polygon(5,3,2);
  polygon(5,4,3);
  polygon(5,1,4);
}

static void polygon(int a, int b, int c)
{
  int i;

  /* insert next polygon */
  i = 4*num_polys;
  data_verts[i]   = 3;
  data_verts[i+1] = a+1; /* .dat file starts with 1, not zero */
  data_verts[i+2] = b+1;
  data_verts[i+3] = c+1;
  link_verts[num_polys][0] = insert_edge(a,b);
  link_verts[num_polys][1] = insert_edge(b,c);
  link_verts[num_polys][2] = insert_edge(c,a);
  num_polys++;
}

static int insert_edge(int a, int b)
{
  int i, first, second;

  if(a<b) { first=a; second=b; }
  else    { first=b; second=a; }
  for(i=0; i<num_edges; i++)
  {
    if(edge_list[i][0]==first && edge_list[i][1]==second) return i;
  }
  edge_list[num_edges][0] = first; /* raw, not +1 */
  edge_list[num_edges][1] = second;
  edge_list[num_edges][2] = 0; /* for mid point vertex */
  num_edges++;
  return num_edges-1;
}

static void split_poly(int i) /* make four faces out of one */
{
  int a, b, c, dl, d, el, e, fl, f; /* original a,b,c  mids  d,e,f */
  /* new d,e,f   a,d,f   b,e,d  c,f,e   */
  a  = data_verts[4*i+1]-1;
  b  = data_verts[4*i+2]-1;
  c  = data_verts[4*i+3]-1;
  dl = link_verts[i][0];
  d  = edge_list[dl][2];
  el = link_verts[i][1];
  e  = edge_list[el][2];
  fl = link_verts[i][2];
  f  = edge_list[fl][2];
  printf("a=%d, b=%d, c=%d, d=%d, e=%d, f=%d, dl=%d, el=%d, fl=%d \n",
         a, b, c, d, e, f, dl, el, fl);
  data_verts[4*i+1] = d+1;
  data_verts[4*i+2] = e+1;
  data_verts[4*i+3] = f+1;
  data_verts[4*num_polys]   = 3;
  data_verts[4*num_polys+1] = a+1;
  data_verts[4*num_polys+2] = d+1;
  data_verts[4*num_polys+3] = f+1;
  num_polys++;
  data_verts[4*num_polys]   = 3;
  data_verts[4*num_polys+1] = b+1;
  data_verts[4*num_polys+2] = e+1;
  data_verts[4*num_polys+3] = d+1;
  num_polys++;
  data_verts[4*num_polys]   = 3;
  data_verts[4*num_polys+1] = c+1;
  data_verts[4*num_polys+2] = f+1;
  data_verts[4*num_polys+3] = e+1;
  num_polys++;
}

static void split(int splits)
{
  int i, j, k, s;
  vert_typ mid;
  vert_typ origin = {0.0, 0.0, 0.0};
  int points[100];
  double len, len_one;
  double avgx, avgy, avgz;
  double lenn = 0.1;
  int npts, pts, pt, newk;
  double nx, ny, nz;
  int debug = 3;

  /* tetrahedron_compute(); */
  octahedron_compute();
  printf("num_points=%d, num_polys=%d, num_edges=%d \n",
         num_points, num_polys, num_edges);
  for(s=0; s<splits; s++)
  {
    printf("splits=%d, s=%d \n", splits, s);
    /* initial .dat structure built, now split */
    /* build link_verts and edge list */
    num_edges = 0;
    for(i=0; i<num_polys; i++)
    {
      link_verts[i][0] = insert_edge(data_verts[4*i+1]-1, data_verts[4*i+2]-1);
      link_verts[i][1] = insert_edge(data_verts[4*i+2]-1, data_verts[4*i+3]-1);
      link_verts[i][2] = insert_edge(data_verts[4*i+3]-1,data_verts[4*i+1]-1);
    }
    for(i=0; i<num_edges; i++)
    {
      printf("edge %d, %d  %d \n", i, edge_list[i][0], edge_list[i][1]);
    }

    for(i=0; i<num_edges; i++) /* split edges */
    {
      j = edge_list[i][0];
      k = edge_list[i][1];
      mid[0]=(data_points[j].x+data_points[k].x)/2.0; /* x */
      mid[1]=(data_points[j].y+data_points[k].y)/2.0; /* y */
      mid[2]=(data_points[j].z+data_points[k].z)/2.0; /* z */
      len = d3len(mid, origin);
      mid[0] = mid[0]/len;
      mid[1] = mid[1]/len;
      mid[2] = mid[2]/len;
      printf("midx=%6.3f, y=%6.3f, z=%6.3f, len=%6.3f \n",
              mid[0], mid[1], mid[2], len);
      edge_list[i][2] = num_points;
      data_points[num_points].x = mid[0];
      data_points[num_points].y = mid[1];
      data_points[num_points].z = mid[2];
      num_points++;
    }

    printf("num_points=%d, num_polys=%d, num_edges=%d \n",
           num_points, num_polys, num_edges);
    for(i=0; i<num_edges; i++)
    {
      printf("edge %d,  v1=%d, v2=%d, mid=%d \n",
             i, edge_list[i][0], edge_list[i][1], edge_list[i][2]);
    }
    for(i=0; i<num_polys; i++)
    {
      printf("link %d, d=%d,  e=%d, f=%d \n",
             i, link_verts[i][0], link_verts[i][1], link_verts[i][2]);
    }


    j = num_polys; /* only loop through old */
    for(i=0; i<j; i++) /* split polys */
    {
      split_poly(i);
    }
    printf("num_points=%d, num_polys=%d, num_edges=%d \n",
           num_points, num_polys, num_edges);
    for(i=0; i<num_polys; i++)
    {
      printf("polygon(%d,%d,%d); \n",
              data_verts[4*i+1]-1, data_verts[4*i+2]-1, data_verts[4*i+3]-1);
    }
    printf("num_points=%d, num_polys=%d, num_edges=%d \n",
           num_points, num_polys, num_edges);
  }
  num_edges = num_edges*4;
  printf("num_points=%d, num_polys=%d, num_edges=%d \n",
         num_points, num_polys, num_edges);

  // spike generator
  if(spike>0)
  {
    printf("making Medieval flail, spiked ball\n");
    // copy data points, more added use new_point, new_poly counts
    // each poly becomes 3 triangles, new point normal from center

    k = 0; // do ccw
    for(i=0; i<num_polys; i++)
    {
      pts = data_verts[k++];
      for(j=0; j<pts; j++)
      {
	 pt = data_verts[k++];
	 points[j] = pt;
      }
      // k is past last point
      for(j=0; j<pts; j++)
      {
	data_verts[k-j-1] = points[j];
      }
    }
    // copy existing points, add new based on normal
    for(i=0; i<num_points; i++)
    {
      new_points[i].x = data_points[i].x;
      new_points[i].y = data_points[i].y;
      new_points[i].z = data_points[i].z;
    }
    new_point = num_points; // first index to place new
    k = 0;
    newk = 0; 
    for(i=0; i<num_polys; i++)
    {
      pts = data_verts[k++]; // should be 3
      new_verts[newk++] = pts;
      new_verts[newk++] = new_point+1; // 1 subtracted on use
      new_verts[newk++] = data_verts[k+2];
      new_verts[newk++] = data_verts[k+0];
      new_verts[newk++] = pts;
      new_verts[newk++] = new_point+1;
      new_verts[newk++] = data_verts[k+1];
      new_verts[newk++] = data_verts[k+2];
      new_verts[newk++] = pts;
      new_verts[newk++] = new_point+1;
      new_verts[newk++] = data_verts[k+0];
      new_verts[newk++] = data_verts[k+1];

      avgx = 0.0;
      avgy = 0.0;
      avgz = 0.0;
      compute_normal(k, &nx, &ny, &nz);
      for(j=0; j<pts; j++)
      {
        pt = data_verts[k++];
        avgx += data_points[pt-1].x;
        avgy += data_points[pt-1].y;
        avgz += data_points[pt-1].z;
        if(debug>2)
	{
          printf("averaging j=%d, k=%d, pt=%d\n", j, k, pt);
          printf("x=%f, y=%f, z=%f\n", data_points[pt-1].x,
                 data_points[pt-1].y, data_points[pt-1].z); 
	}
      }
      avgx /= (double)pts;
      avgy /= (double)pts;
      avgz /= (double)pts;;
      if(debug>1) printf("avgx=%f, avgy=%f, avgz=%f\n", avgx, avgy, avgz);
      new_points[new_point].x = avgx+lenn*nx;
      new_points[new_point].y = avgy+lenn*ny;
      new_points[new_point].z = avgz+lenn*nz;
      new_point++;
      new_polys = new_polys+3; // added 3 triangles and one vertex
      printf("new_point=%d, new_polys=%d, newk=%d, k=%d\n",
             new_point, new_polys, newk, k);
      printf("new x=%f, y=%f, z=%f\n", new_points[new_point-1].x,
	     new_points[new_point-1].y, new_points[new_point-1].z);
    } // end num_polys
  } // end spike
} // end split

static void compute_normal(int ivert, 
                           double *nxout, double *nyout, double *nzout)
{
  int j, k, pt;
  double ax, ay, az, bx, by, bz;
  double nx, ny, nz, s;
  double v[3][3];

  k = ivert;
  for(j=0; j<3; j++) /* use first three points, assume planar */
  {
    pt = data_verts[k++];
    v[j][0] = data_points[pt-1].x; // 2,1,0 if backwards
    v[j][1] = data_points[pt-1].y;
    v[j][2] = data_points[pt-1].z;
  }
  /* 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 = sqrt(nx*nx+ny*ny+nz*nz);
  nx = nx / s;
  ny = ny / s;
  nz = nz / s;
  *nxout = nx;
  *nyout = ny;
  *nzout = nz;
} /* end compute normal */