/* make_cyl_spiral_635.c  write out a cyl_spiral_635.dat for run6 */

#include <math.h>
#include <stdio.h> 

#define LINES  128     /* storage max for nx       */
#define POINTS 128     /* storage max for ny       */
#define TRIS 16384     /* storage max for points   */

int main(int argc, char *argv[])
{
  struct point { float x, y, z; int index; };
  struct point raw[LINES][POINTS];
  struct tri_poly { int p1, p2, p3; };
  struct tri_poly tri_raw[TRIS];
  FILE *fp;
  int i, j;
  int points, polys;
  double x, y, z;
  int nx = 63;     /* last == first */
  int ny = 15;
  double pi = 3.141592653589793238462643383279502884197;

  /* p306 curves and surfaces */
  double a = 0.1; /* inner radius */
  double b = 1.0; /* max z */
  double c = 0.5; /* outer radius */
  double n = 3;   /* cycles */
  double u = 0.0; /* parametric in u and v */
  double du = 2.0*pi/(double)ny;
  double v = 0.0;
  double dv = 2.0*pi/(double)(nx-1);
  
  printf("make_cyl_spiral_635 running \n");
  
  points = 1; /* thats the standard */

  for(i=0; i<nx; i++)
  {
    u = 0.0;
    for(j=0; j<ny; j++)
    {
      raw[i][j].x = a*cos(n*v)*(1.0+cos(u))+c*cos(n*v);
      raw[i][j].y = a*sin(n*v)*(1.0+cos(u))+c*sin(n*v);
      raw[i][j].z = b*v/(2.0*pi) + a*sin(u);
      raw[i][j].index = points;
      points++;
      u = u + du;
    }
    v = v + dv;
  }
  points--; /* fix up count */
  printf("raw built \n");

  polys = 0;
  for(i=0; i<(nx-1); i++)
  {
    for(j=0; j<ny; j++)
    {
      tri_raw[polys].p1 = raw[i][j].index;
      tri_raw[polys].p2 = raw[i][(j+1)%ny].index;
      tri_raw[polys].p3 = raw[(i+1)%nx][(j+1)%ny].index;
      polys++;
      tri_raw[polys].p1 = raw[i][j].index;
      tri_raw[polys].p2 = raw[(i+1)%nx][(j+1)%ny].index;
      tri_raw[polys].p3 = raw[(i+1)%nx][j].index;
      polys++;
    }
  }
  printf("tri_raw built \n");
  printf("Outputing points= %d polys %d \n", points, polys);
  fp = fopen("cyl_spiral_635.dat","w");
  if(fp==NULL)
  {
    printf("Can not open %s \n", "cyl_spiral_635.dat");
    return 1;
  }
  fprintf(fp, "%d %d \n", points, polys);
  for(i=0; i<nx; i++)   /* loop around surface */
  {
    for(j=0; j<ny; j++)  /* walk around cross section generating skin  */
    {
      fprintf(fp,"%f %f %f \n", raw[i][j].x, raw[i][j].y, raw[i][j].z);
    }
  }
  for(i=0; i<polys; i++)   /* now output 3 point polygons */
  {
    fprintf(fp, "3 %d %d %d \n", tri_raw[i].p1, tri_raw[i].p2, tri_raw[i].p3); 
  }
  fclose(fp);
  printf("wrote  cyl_spiral_635.dat  for  run6 \n");
  return 0;
} /* end of make_cyl_spiral_635.c */