/* make_helix.c  write out a helix_635.dat for run4 */
/* actually, a troidal helix, the toroid is made by make_torrus_635.c */

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

#define LINES  259      /* storage max for nx  */
#define POINTS 259      /* storage max for ny */
#define TRIS 68000      /* storage max for points   */


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 debug;             /* 0 for no debug print in xterm */
  float phi, theta;
  float phi1, theta1;
  int i, j;
  int points, polys;
  int status;
  float x0, y0, z0;
  float x1, y1, z1;
  float x2, y2, z2;
  float x3, y3, z3;
  float x2d, y2d, z2d, r2d, r2n;
  float vx2, vy2, vz2;
  float rx2, ry2, rz2;
  double pi = 3.141592653589793238462643383279502884197;
  float r1 = 8.0;  /* major radius of torrus */
  float r2 = 4.0;  /* minor radius of torrus */
  float r3 = 1.0;  /* minor radius of helix */
  float F  = 8.0;  /* wrapping factor of r2 around r1 */

  /*   for big helix  */
  float phi_step = pi/128.0;
  int nx = 257;

  float theta_step = pi/8.0;
  int ny = 17;

  /* for smooth shaded helix, smaller steps */
  phi_step = pi/64.0;
  nx = 129;
  theta_step = pi/8.0;
  ny = 17;


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

  phi1 = 0.0;
  for(i=0; i<nx; i++)   /* loop around toride at radius r2, around that at r2 */
  {                     /* this makes x1+x2, the center of the generated figure */
    phi = phi1;
    phi1 = phi1 + phi_step;
    x1 =r1*sin(phi);
    y1 =r1*cos(phi);
    z1 = 0.0;

    x2 = x1 + r2*sin(phi)*cos(phi*F); /* F is number of helix loops */
    y2 = y1 + r2*cos(phi)*cos(phi*F);
    z2 = z1 + r2*         sin(phi*F);

    /* the derivitive of x1+x2 to get velocity vector direction */
    x2d =  r1*cos(phi) + r2*cos(phi)*cos(phi*F) - F*r2*sin(phi)*sin(phi*F);
    y2d = -r1*sin(phi) - r2*sin(phi)*cos(phi*F) - F*r2*cos(phi)*sin(phi*F);
    z2d = F*r2*cos(phi*F);
    r2d = sqrt(x2d*x2d+y2d*y2d+z2d*z2d); /* normalize */
    x2d = x2d/r2d;
    y2d = y2d/r2d;
    z2d = z2d/r2d;

    /* r2,x2,y2,z2 only, thus subtract out r1,x1,y1,z1 */
    r2n = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
    rx2 = (x2-x1)/r2n;
    ry2 = (y2-y1)/r2n;
    rz2 = (z2-z1)/r2n; /* now have r2 vector, normal to velocity */

    vx2 = ry2*z2d - rz2*y2d; /* cross product */
    vy2 = rz2*x2d - rx2*z2d;
    vz2 = rx2*y2d - ry2*x2d; /* this and r2 vector define plane of cross section */

    theta1 = 0.0;
    for(j=0; j<ny; j++)  /* walk around cross section generating skin */
    {
      theta = theta1;
      theta1 = theta1 + theta_step;
      x3 = r3*(rx2*cos(theta) + vx2*sin(theta));
      y3 = r3*(ry2*cos(theta) + vy2*sin(theta));
      z3 = r3*(rz2*cos(theta) + vz2*sin(theta));

      raw[i][j].x =x2+x3;  /* actually sum of x1+x2+x3, the final point on surface */
      raw[i][j].y =y2+y3;
      raw[i][j].z =z2+z3;
      raw[i][j].index = points;
      points++;
    }
  }
  points--;
  printf("raw built \n");

  polys = 0;             /* now build set of points defining surface */
  for(i=0; i<(nx-1); i++)       
  {
    /* reuse last point of i and j as first point */
    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("helix_635.dat","w");
  if(fp==NULL)
  {
    printf("Can not open %s \n", "helix_635.dat");
    return 1;
  }
  fprintf(fp, "%d %d \n", points, polys);
  for(i=0; i<nx; i++)   /* loop around toride at radius r2, around that at r2 */
  {                     /* this makes x1+x2, the center of the generated figure */

    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);

  return 0;
} /* end of make_helix_635.c */