/* make_spiral_635.c  write out a spiral_635.dat for run6 */
/* spiral trough  X and Y -8 to 8,  Z about -0.5 to 8.8   */

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

#define LINES  258     /* storage max for nx       */
#define POINTS 258     /* storage max for ny       */
#define TRIS 66000     /* storage max for points   */

static double spiral(double x, double y); /* z = */

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;

  float x0 = -8.0;
  float dx = 0.5;
  int nx = 33;

  float y0 = -8.0;
  float dy = 0.5;
  float ny = 33;

  printf("make_spiral_635 running \n");
  
  points = 1; /* thats the standard */

  x = x0 - dx;
  for(i=0; i<nx; i++)
  {
    x = x + dx;
    y = y0 - dy;
    for(j=0; j<ny; j++)
    {
      y = y + dy;
      z = (float)spiral((double)x, (double)y);
      raw[i][j].x = x;
      raw[i][j].y = y;
      raw[i][j].z = z;
      raw[i][j].index = points;
      points++;
    }
  }

  points--; /* fix up count */
  printf("raw built \n");

  polys = 0;
  for(i=0; i<(nx-1); i++)
  {
    for(j=0; j<(ny-1); j++)
    {
      tri_raw[polys].p1 = raw[i][j].index;
      tri_raw[polys].p2 = raw[i][j+1].index;
      tri_raw[polys].p3 = raw[i+1][j+1].index;
      polys++;
      tri_raw[polys].p1 = raw[i][j].index;
      tri_raw[polys].p2 = raw[i+1][j+1].index;
      tri_raw[polys].p3 = raw[i+1][j].index;
      polys++;
    }
  }
  printf("tri_raw built \n");
  printf("Outputing points= %d polys %d \n", points, polys);
  fp = fopen("spiral_635.dat","w");
  if(fp==NULL)
  {
    printf("Can not open %s \n", "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  spiral_635.dat  for  run6 \n");
  return 0;
} /* end of make_spiral_635.c */

static double spiral(double x, double y)
{
  double halfpi = 1.57079;
  double r, theta, tz;

  r = sqrt(x*x+y*y);
  if( r < 0.01) theta = 0.0;
  else          theta = atan2(y,x);
  tz = theta/halfpi;
  if(tz < 0.0) tz = tz + 4.0;
  return 0.75*r + sin((r+tz)*halfpi);
} /* end of spiral */