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

#define R0		0.03
#define R1		1.0
#if 0
/* Big shell parameters */
#define PTSPERCIRC	16
#define CIRCPERREV	20
#define REVS		8
#else 
/* Small shell parameters */
#define R0		1.0
#define R1		1.0
#define PTSPERCIRC	12
#define CIRCPERREV	16
#define REVS		1.0
#endif


#define FRACTIONDY	1.5	/* number of radii to drop per revolution */
#define FRACTDECAY	1.0

void main()
{
  int ptn, circn, revn;
  int n, npt, npt2, ncirc, ncirc2;
  float r, tuber;
  float theta, phi;
  float dthetadcirc = (2.0 * 3.1415926536) / CIRCPERREV;
  float dphidpt = (2.0 * 3.1415926536) / PTSPERCIRC;
  float y;
  float ptx, pty, ptz;
  int mode;
  float drdeps = (R1 - R0) / ((float) CIRCPERREV * REVS);
  float drdcirc, ddrdcirc;

  drdcirc = exp(log(R1/R0)/ (CIRCPERREV * REVS));
  ddrdcirc = exp(log(FRACTDECAY) /  (CIRCPERREV * REVS));
  
#define POINT 	0
#define NORMAL 	1
#define DIFFUSE 2
#define SPECULAR 3
#define INDEX 	4

  printf("#Inventor V2.0 ascii\n"
	 "Separator {\n");


  for (mode = POINT; mode <= INDEX; mode++) {

    r =R0;
    theta = 0;
    phi = 0;
    y = 0;
    n=0; 
    

    /* Print initial info to tell the .iv file which component we
       are specifying.  The program makes 5 passes: position, normal, 
       color, specular, or the vertex numbers to connect for each
       triangle.
       */
    switch (mode) {
    case POINT:
      printf("Coordinate3 {\n"
	     "point [\n");
      break;
    case NORMAL:
      printf("NormalBinding {\n"
	     "value	PER_VERTEX_INDEXED\n"
	     "}\n"
	     "Normal {\n"
	     "vector [\n");
      break;
    case DIFFUSE:
      printf("MaterialBinding {\n"
	     "value       PER_VERTEX_INDEXED\n"
	     "}\n"
	     "Material {\n"
	     "diffuseColor [\n");
      break;
    case SPECULAR:
      printf("specularColor [\n");
      break;
    case INDEX:
      printf("IndexedFaceSet {\n"
	     "coordIndex [\n");
      break;
    }

    /* Now go through and, for each vertex, print out the relevant
       info.
       */
    for (revn = 0; revn < REVS; revn++) {
      for (circn = 0; circn < CIRCPERREV; circn++) {
	/* Increment, unless we're finishing the last...*/
	if (revn < REVS-1 || circn < CIRCPERREV-2) {
	  theta += dthetadcirc;
	  if (revn > 0) {
	    y -= (FRACTIONDY * r) / CIRCPERREV;
	  }
	  r *= drdcirc;
	  drdcirc *= ddrdcirc;
	  tuber=r;
	} else {
	  /* last ring. make r 0 */
	  tuber=0;
	}
	for (ptn = 0; ptn < PTSPERCIRC; ptn++) {
	  switch (mode) {
	  case POINT:
	    ptx = cos(theta) * (r + tuber*cos(phi));
	    pty = y + tuber*sin(phi);
	    ptz = sin(theta) * (r + tuber*cos(phi));
	    printf("%f %f %f,\n", ptx, pty, ptz);
	    break;
	  case NORMAL:
	    ptx = cos(theta)*cos(phi);
	    pty = sin(phi);
	    ptz = sin(theta)*cos(phi);
	    printf("%f %f %f,\n", ptx, pty, ptz);
	    break;
	  case DIFFUSE:
	    ptx = .2+.7*cos(-1.33*theta+phi)*cos(-1.33*theta+phi);
	    pty =    .9*cos(-1.33*theta+phi)*cos(-1.33*theta+phi);
	    ptz =    .9*cos(-1.33*theta+phi)*cos(-1.33*theta+phi);
	    printf("%f %f %f,\n", ptx, pty, ptz);
	    break;
	  case SPECULAR:
	    ptx = .05 + .3*(1.0-cos(-1.33*theta+phi)*cos(-1.33*theta+phi));
	    pty = .05 + .3*(1.0-cos(-1.33*theta+phi)*cos(-1.33*theta+phi));
	    ptz = .05 + .3*(1.0-cos(-1.33*theta+phi)*cos(-1.33*theta+phi));
	    printf("%f %f %f,\n", ptx, pty, ptz);

	    break;
	  case INDEX:
	    npt = n%PTSPERCIRC;
	    npt2 = (n+PTSPERCIRC-1)% PTSPERCIRC;
	    ncirc = revn*PTSPERCIRC*CIRCPERREV + PTSPERCIRC * circn;
	    ncirc2 = revn*PTSPERCIRC*CIRCPERREV + PTSPERCIRC * (circn+1);
	    if (n < PTSPERCIRC * CIRCPERREV * REVS - CIRCPERREV+2) {
	      printf("\t%d, %d, %d, -1,   \t%d, %d, %d, -1,\n",
		     npt2+ncirc, 
		     npt+ncirc, 
		     npt + ncirc2,
		     npt2+ncirc2, 
		     npt2+ncirc, 
		     npt + ncirc2);	
	    }
	    n++;
	    break;
	  }
	    
	  phi += dphidpt;

	}
      }
    }
    switch (mode) {
    case POINT:
    case NORMAL:
    case SPECULAR:
    case INDEX:
      printf("]}\n");
      break;
    case DIFFUSE:
      printf("]\n");
      break;
    }
  }

  printf("}\n");
  
}