#define NRING 50

#define ROOT_HALF .707106781186547524400844362104849039284835937688474036588
#define ROOT_3_4 .866025403784438646763723170752936183471402626905190314

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

extern const char *__progname;

#include "1d.h"
#include "3d.h"
#include "allocv.h"

typedef struct csrv CSRV;
typedef struct tri TRI;
typedef struct point POINT;

struct csrv {
  PT3 cs;
  PT3 norm;
  } ;

struct point {
  int pt;
  PT3 normal;
  } ;

struct tri {
  int material;
#define MAT_WHITE 0
#define MAT_LIGHT 1
#define MAT_DARK  2
  int px[3];
  } ;

static int ring0[NRING+1];
static int ring0_offset;
static int ring1[NRING+1];
static int ring2[NRING+1];
static int *rcur;
static int *rprev;
static int ring_offset;

static ALLOCV *pts;
static ALLOCV *points;
static ALLOCV *tris;

static double length;
static double slant;
static double limpness;
static double dent_offset;
static double dent_depth_a;
static double dent_depth_b;
static double dent_width;
static double swell;
static double swellpower;
static double tipslant;
static double tiplen;
static double holeadj1;
static double holeadj2;

static double end;
static PT3 basecentre;
static PT3 shaftdir;
static PT3 shaftdn;

static int cur_material;

#define DEGTORAD(a) (((a)*M_PI)/180)
#define RADTODEG(a) (((a)*180)/M_PI)

static void init_arrays(void)
{
 pts = allocv_init(sizeof(PT3));
 points = allocv_init(sizeof(POINT));
 tris = allocv_init(sizeof(TRI));
}

static int newpt(PT3 loc)
{
 int x;

 x = allocv_new(pts);
 *(PT3 *)allocv_el(pts,x) = loc;
 return(x);
}

static int newpoint(int loc, PT3 norm)
{
 int x;

 if (norm.x != norm.x) abort();
 x = allocv_new(points);
 if (! allocv_okx(pts,loc)) abort();
 *(POINT *)allocv_el(points,x) = (POINT){ .pt=loc, .normal=unit3(norm) };
 return(x);
}

static int newtri(int a, int b, int c)
{
 int x;

 x = allocv_new(tris);
 if ( !allocv_okx(points,a) ||
      !allocv_okx(points,b) ||
      !allocv_okx(points,c) ) abort();
 *(TRI *)allocv_el(tris,x) = (TRI){ .material=cur_material,
				    .px={ [0]=a, [1]=b, [2]=c } };
 return(x);
}

/*
static void newtri2(int a, int b, int c)
{
 newtri(a,b,c);
 newtri(a,c,b);
}
*/

static double swellbase(double a)
{
 return(.1 + (tipslant * .5 * (1 - cos(a*2*M_PI))));
}

static double swellbase_max(void)
{
 return(.1 + tipslant);
}

static double swellbase_d(double a)
{
 double d;

 d = M_PI * tipslant * sin(a*M_PI);
 if (a > .5) d = -d;
 if (a == .5) d = 0;
 return(d);
}

static double dent_r(double f, double t)
{
 return(1 - (f*( exp(-sqr((t/dent_width)+dent_offset)) +
		 exp(-sqr((t/dent_width)-dent_offset)) )));
}

static double dent_drdt(double f, double t)
{
 return(-f*( ( exp(-sqr((t/dent_width)+dent_offset)) *
	       -(2/dent_width)*((t/dent_width)+dent_offset) ) +
	     ( exp(-sqr((t/dent_width)-dent_offset)) *
	       -(2/dent_width)*((t/dent_width)-dent_offset) ) ));
}

static double swell_amount(double t)
{
 return(swell * (1 - pow((1-(2*fabs(t))),swellpower)));
}

static CSRV shaft_cs(double t, double l)
{
 double a;
 double r;
 double drdt;
 double c;
 double s;
 double f;

 t -= .5;
 a = t * M_PI * 2;
 c = cos(a);
 s = sin(a);
 f = interp1(dent_depth_a,l,dent_depth_b);
 r = dent_r(f,t);
 drdt = dent_drdt(f,t);
 return((CSRV){ .cs = MAKEPT3(r*c,r*s,0),
		.norm = unit3(MAKEPT3((drdt*s)+(r*M_PI*2*c),-((drdt*c)+(r*M_PI*2*-s)),0)) });
}

static CSRV swell_cs(double t0, int which)
{
 double t;
 double a;
 double r;
 double c;
 double s;
 double sw;
 double z;
 PT3 norm;

 static PT3 norm_0(void)
  { return(unit3(add3(smul3(r*swellbase_d(t0)/(M_PI*2),
			    unit3(MAKEPT3(r*M_PI*2*-s,r*M_PI*2*c,0))),
		      MAKEPT3(0,0,-1))));
  }

 static PT3 norm_2(void)
  { return(unit3(MAKEPT3(r*M_PI*2*c,-(r*M_PI*2*-s),0)));
  }

 t = t0 - .5;
 a = t * M_PI * 2;
 c = cos(a);
 s = sin(a);
 r = 1;
 sw = swell_amount(t);
 z = swellbase(t0);
 switch (which)
  { case 0:
       r += sw / 3;
       norm = norm_0();
       break;
    case 1:
       r += sw * ((1/3.0) + (ROOT_HALF*2/3.0));
       z += sw * (1-ROOT_HALF)*2/3.0;
       norm = unit3(add3(norm_0(),norm_2()));
       break;
    case 2:
       r += sw;
       z += sw * 2 / 3.0;
       norm = norm_2();
       break;
  }
 return((CSRV){ .cs = MAKEPT3(r*c,r*s,z),
		.norm = norm });
}

static CSRV end_cs(double t0, double frac)
{
 double t;
 double a;
 double c;
 double s;
 double r0;
 double sw;
 double base;
 double r;
 PT3 norm;

 t = t0 - .5;
 a = t * M_PI * 2;
 c = cos(a);
 s = sin(a);
 sw = interp1(swell_amount(t),frac*4,swell);
 if (sw > swell) sw = swell;
 r0 = 1 + sw;
 base = swellbase(t0) + (sw * 2 / 3.0);
 r = r0 * sqrt(1-frac);
 norm = unit3(MAKEPT3(r*M_PI*2*c,-(r*M_PI*2*-s),0));
 norm.z = r0 / (2 * sqrt(1-frac));
 return((CSRV){ .cs = MAKEPT3(r*c,r*s,interp1(base,frac,end)),
		.norm = unit3(norm) });
}

static CSRV hole_cs(double t0, int which)
{
 double t;
 double a;
 double c;
 double s;

 t = t0 - .5;
 a = t * M_PI * 2;
 c = cos(a);
 s = sin(a);
 switch (which)
  { case 0:
       return((CSRV){ .cs = MAKEPT3(c*.2,s*.03,end+interp1(holeadj1,(1+c)/2,holeadj2)),
		      .norm = MAKEPT3(0,0,1) });
       break;
    case 1:
       return((CSRV){ .cs = MAKEPT3(c*.15,s*.02,end+interp1(holeadj1,(1+c)/2,holeadj2)-.02),
		      .norm = unit3(MAKEPT3(-c,-s,.75)) });
       break;
  }
 abort();
}

static void initring(PT3 cen, double zdist, PT3 axis, PT3 side)
{
 int i;
 CSRV cs;
 PT3 rx;
 PT3 ry;
 PT3 r;

 axis = unit3(axis);
 cen = add3(cen,smul3(zdist,axis));
 rx = unit3(sub3(side,smul3(dot3(side,axis),axis)));
 ry = unit3(cross(axis,rx));
 for (i=NRING-1;i>=0;i--)
  { cs = shaft_cs(i/(double)NRING,0);
    r = add3(smul3(cs.cs.x,rx),smul3(cs.cs.y,ry));
    ring1[i] = newpoint(newpt(add3(cen,r)),add3(add3(smul3(cs.norm.x,rx),smul3(cs.norm.y,ry)),smul3(cs.norm.z,axis)));
  }
 ring1[NRING] = ring1[0];
 rcur = &ring1[0];
 rprev = &ring2[0];
 ring_offset = 0;
}

static void swaprings(void)
{
 int *t;

 t = rcur;
 rcur = rprev;
 rprev = t;
}

static void ring(PT3 cen, double zdist, PT3 axis, PT3 side, CSRV (*csf)(double))
{
 int i;
 CSRV cs;
 PT3 rx;
 PT3 ry;
 PT3 r;

 ring_offset = (ring_offset & 1) ^ 1;
 axis = unit3(axis);
 cen = add3(cen,smul3(zdist,axis));
 rx = unit3(sub3(side,smul3(dot3(side,axis),axis)));
 ry = unit3(cross(axis,rx));
 swaprings();
 for (i=NRING-1;i>=0;i--)
  { cs = (*csf)((i+i+ring_offset)/(2.0*NRING));
    r = add3(add3(smul3(cs.cs.x,rx),smul3(cs.cs.y,ry)),smul3(cs.cs.z,axis));
    rcur[i] = newpoint( newpt(add3(cen,r)),
			add3( add3( smul3(cs.norm.x,rx),
				    smul3(cs.norm.y,ry) ),
			      smul3(cs.norm.z,axis) ) );
  }
 rcur[NRING] = rcur[0];
 for (i=NRING-1;i>=0;i--)
  { newtri(rcur[i],rcur[i+1],rprev[i+ring_offset]);
    newtri(rcur[i+(ring_offset^1)],rprev[i+1],rprev[i]);
  }
}

static void endcap(PT3 cen, double zdist, PT3 axis)
{
 int i;
 int cx;

 cx = newpoint(newpt(add3(cen,smul3(zdist,axis))),axis);
 for (i=NRING;i>0;i--) newtri(cx,rcur[i],rcur[i-1]);
}

static void dump(void)
{
 int i;
 int n;

 printf("ModelV1\n");
 printf("3 %d %d %d\n",allocv_n(pts),allocv_n(points),allocv_n(tris));
 printf("White 0 0 0 .7 .7 .7 10 .5\n");
 printf("Light 0 0 0 .4 .4 .4 10 .1\n");
 printf("Dark 0 0 0 .2 .2 .2 0 0\n");
 n = allocv_n(pts);
 for (i=0;i<n;i++)
  { PT3 *e;
    e = allocv_el(pts,i);
    printf("%g %g %g\n",e->x,e->y,e->z);
  }
 n = allocv_n(points);
 for (i=0;i<n;i++)
  { POINT *e;
    e = allocv_el(points,i);
    printf("%d s %g %g %g\n",e->pt,e->normal.x,e->normal.y,e->normal.z);
  }
 n = allocv_n(tris);
 for (i=0;i<n;i++)
  { TRI *e;
    e = allocv_el(tris,i);
    printf("%d %d %d %d\n",e->material,e->px[0],e->px[1],e->px[2]);
  }
}

#if 0
static void arcpoints(PT3 cen, PT3 zaxis, PT3 xaxis, double radius,
 double angmin, double angmax, int npoints, PT3 *points)
{
 PT3 yaxis;
 int i;
 double a;

 if (npoints < 2) abort();
 npoints --;
 zaxis = unit3(zaxis);
 xaxis = unit3(sub3(xaxis,smul3(dot3(xaxis,zaxis),zaxis)));
 yaxis = unit3(cross(zaxis,xaxis));
 for (i=npoints;i>=0;i--)
  { a = interp1(angmin,i/(double)npoints,angmax);
    points[i] = add3( cen,
		      add3( smul3(cos(a)*radius,xaxis),
			    smul3(sin(a)*radius,yaxis) ) );
  }
}
#endif

int main(int, char **);
int main(int ac, char **av)
{
 int i;
 PT3 cen;
 PT3 dir;
 PT3 dn;
 PT3 rt;

 if (ac != 14)
  { fprintf(stderr,"Usage: %s\
 length\
 slant\
 limpness\
 dent-offset\
 dent-depth1\
 dent-depth2\
 dent-width\
 swell\
 swellpower\
 tipslant\
 tiplen\
 holeadj1\
 holeadj2\
\n",__progname);
    exit(1);
  }
 i = 1;
 length = atof(av[i++]);
 slant = atof(av[i++]);
 limpness = atof(av[i++]);
 dent_offset = atof(av[i++]);
 dent_depth_a = atof(av[i++]);
 dent_depth_b = atof(av[i++]);
 dent_width = atof(av[i++]);
 swell = atof(av[i++]);
 swellpower = atof(av[i++]);
 tipslant = atof(av[i++]);
 tiplen = atof(av[i++]);
 holeadj1 = atof(av[i++]);
 holeadj2 = atof(av[i++]);
 end = swellbase_max() + swell + tiplen;
 shaftdir = MAKEPT3(sin(DEGTORAD(slant)),0,cos(DEGTORAD(slant)));
 shaftdn = MAKEPT3(-shaftdir.z,0,shaftdir.x);
 basecentre = smul3(.5,shaftdir);
 cur_material = MAT_WHITE;
 init_arrays();
 initring(basecentre,0,shaftdir,shaftdn);
 for (i=NRING;i>=0;i--) ring0[i] = rcur[i];
 ring0_offset = ring_offset;
 cen = basecentre;
 dir = shaftdir;
 dn = shaftdn;
 rt = MAKEPT3(0,1,0);
#ifdef STRAIGHT_SHAFT
 ring( cen, 0, dir, dn,
       ({ static CSRV foo(double a)
	   { return(shaft_cs(a,0));
	   }
	  &foo;
	}) );
 cen = add3(cen,smul3(length,dir));
#else
 for (i=1;i<=32;i++)
  { double cura;
    double la;
    cen = add3(cen,smul3(length/32,dir));
    cura = asin(norm3(cross(dir,MAKEPT3(-1,0,0))));
    la = limpness * length / 32;
    dir = unit3(rotate3(dir,rt,(cura<la)?cura:la));
    dn = unit3(sub3(dn,smul3(dot3(dn,dir),dir)));
    ring( cen, 0, dir, dn,
	  ({ static CSRV foo(double a)
	      { return(shaft_cs(a,i/32.0));
	      }
	     &foo;
	   }) );
  }
#endif
 ring( cen, 0, dir, dn,
       ({ static CSRV foo(double a)
	   { return(shaft_cs(a,1));
	   }
	  &foo;
	}) );
 ring( cen, 0, dir, dn,
       ({ static CSRV foo(double a)
	   { CSRV v;
	     v = shaft_cs(a,1);
	     v.cs.z += swellbase(a);
	     return(v);
	   }
	  &foo;
	}) );
 for (i=0;i<NRING;i++) rcur[i] = newpoint(((POINT *)allocv_el(points,rcur[i]))->pt,MAKEPT3(-1,0,0));
 rcur[NRING] = rcur[0];
 for (i=0;i<3;i++)
  { ring( cen, 0, dir, dn,
	  ({ static CSRV foo(double a)
	      { return(swell_cs(a,i));
	      }
	     &foo;
	   }) );
  }
 for (i=1;i<8;i++)
  { ring( cen, 0, dir, dn,
	  ({ static CSRV foo(double a)
	      { return(end_cs(a,i/8.0));
	      }
	     &foo;
	   }) );
  }
 ring( cen, 0, dir, dn,
       ({ static CSRV foo(double a)
	   { return(hole_cs(a,0));
	   }
	  &foo;
	}) );
 cur_material = MAT_DARK;
 endcap(cen,end-tiplen,dir);
/*
 cur_material = MAT_WHITE;
 ring_offset = ring0_offset;
 for (i=NRING;i>=0;i--) rcur[i] = ring0[i];
*/
 dump();
 return(0);
}