/*
 * Args, with example values:
 *	50	nring (count, triangles per ring)
 *	20	shaftrings (count, rings in main shaft)
 *	6	diprings (count, rings in dip)
 *	4	swellrings (count, rings in swell)
 *	8	tiprings (count, rings in tip)
 *	3	dimplerings (count, rings in dimple)
 *	1	radius (distance)
 *	10	length (distance)
 *	20	upslant (degrees, 0=horizontal, 90=straight up)
 *	.1	stiffness (1=straight, 0=limp)
 *	30	headslant (degrees, 0=square across)
 *	5	dip (just below head, percent of radius)
 *	.05	diplen1 (distance, swell edge to min radius)
 *	.15	diplen2 (distance, normal shaft to min radius)
 *	8	swell (head edge, percent of radius)
 *	.05	swellround (radius of curvature of swell edge)
 *	.3	tiplen (distance, notch max to tip)
 *	.5	tipround (radius of curvature of tip)
 *	10	notchwidth (circumferential degrees, centre to edge)
 *	.3	dimpleoffset (distance, from centre towards notch)
 *	.4	dimplelen (distance)
 *	.1	dimplewidth (distance)
 *	.1	dimpledepth (distance)
 *
 * Origin is centre of base.  upslant=0 extends towards +X.  Gravity
 *  pulls towards -Z.
 */

#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;
static int ring0_offset;
static int *ring1;
static int *ring2;
static int *rcur;
static int *rprev;
static int ring_offset;

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

static int nring;
static int shaftrings;
static int diprings;
static int swellrings;
static int tiprings;
static int dimplerings;
static double radius;
static double length;
static double upslant;
static double stiffness;
static double headslant;
static double dip;
static double diplen1;
static double diplen2;
static double swell;
static double swellround;
static double tiplen;
static double tipround;
static double notchwidth;
static double dimpleoffset;
static double dimplelen;
static double dimplewidth;
static double dimpledepth;

static int newpt(PT3 loc)
{
 int x;

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

#if 0

static int cur_material;

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

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

#endif

static void unreasonable(const char *)
	__attribute__((__noreturn__));
static void unreasonable(const char *msg)
{
 fprintf(stderr,"%s: unreasonable arguments: (%s)\n",__progname,msg);
 exit(1);
}

static void check_args(void)
{
 if (nring < 4) unreasonable("nring must be >= 4");
 if (shaftrings < 2) unreasonable("shaftrings must be >= 2");
 if (diprings < 2) unreasonable("diprings must be >= 2");
 if (swellrings < 2) unreasonable("swellrings must be >= 2");
 if (tiprings < 2) unreasonable("tiprings must be >= 2");
 if (dimplerings < 1) unreasonable("dimplerings must be >= 1");
 if (radius < 1e-3) unreasonable("radius must be >= .001");
 if (length < 1e-3) unreasonable("length must be >= .001");
 if ((upslant < -89.9) || (upslant > 89.9)) unreasonable("upslant must be [-89.9..89.9]");
 if ((stiffness < 0) || (stiffness > 1)) unreasonable("stiffness must be [0..1]");
 if ((headslant < 0) || (headslant > 80)) unreasonable("headslant must be [0..80]");
 if ((dip < -100) || (dip > 100)) unreasonable("dip must be [-100..100]");
 if (diplen1 < 0) unreasonable("diplen1 must be >= 0");
 if (diplen2 < 0) unreasonable("diplen2 must be >= 0");
}

static void init_arrays(void)
{
 ring0 = malloc((nring+1)*sizeof(int));
 ring1 = malloc((nring+1)*sizeof(int));
 ring2 = malloc((nring+1)*sizeof(int));
 pts = allocv_init(sizeof(PT3));
 points = allocv_init(sizeof(POINT));
 tris = allocv_init(sizeof(TRI));
}

static void genring(PT3 c, PT3 x, PT3 y, int offset)
{
 int i;
 double a;
 double s;
 double c;
 PT3 r;

 for (i=0;i<nring;i++)
  { a = ((i + (offset?.5:0)) * 2 * M_PI) / nring;
    s = sin(a);
    c = cos(a);
    r = add3(smul3(c,x),smul3(s,y));
    rcur[i] = newpoint(add3(c,r),r);
  }
}

static void gen_shaft(void)
{
 PT3 o;
 PT3 x;
 PT3 y;
 PT3 z;
 int i;

 ring0_offset = 0;
 ring_offset = 0;
 o = MAKEPT3(0,0,0);
 z = MAKEPT3(cos(DEGTORAD(upslant)),0,sin(DEGTORAD(upslant)));
 y = MAKEPT3(0,1,0);
 x = cross(y,z);
 rcur = ring0;
 rprev = ring1;
 genring(o,x,y,0);
 for (i=shaftrings;i>0;i--)
  { 
stiffness = 0 -> straight, radius of curvature = infinity
stiffness = 1 -> limp, radius of curvature = shaft radius
radius of curvature = shaft radius / stiffness
  }
}

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

 if (ac != 24)
  { fprintf(stderr,"\
Arguments:\n\
 nring (count, triangles per ring)\n\
 shaftrings (count, rings in main shaft)\n\
 diprings (count, rings in dip)\n\
 swellrings (count, rings in swell)\n\
 tiprings (count, rings in tip)\n\
 dimplerings (count, rings in dimple)\n\
 radius (distance)\n\
 length (distance)\n\
 upslant (degrees above horizontal)\n\
 stiffness (0=straight, 1=limp)\n\
 headslant (degrees, 0=straight across)\n\
 dip (below head, percent of radius)\n\
 diplen1 (distance, swell edge to min radius)\n\
 diplen2 (distance, min radius to shaft)\n\
 swell (head edge, percent of radius)\n\
 swellround (distance, radius of curvature of swell edge)\n\
 tiplen (distance, notch to tip)\n\
 tipround (radius of curvature of tip)\n\
 notchwidth (circumferential degrees, middle to edge)\n\
 dimpleoffset (distance, centreline towards notch)\n\
 dimplelen (distance, length of dimple)\n\
 dimplewidth (distance, width of dimple)\n\
 dimpledepth (distance, depth of dimple)\n\
\n");
    exit(1);
  }
 i = 1;
 nring = atoi(av[i++]);
 shaftrings = atoi(av[i++]);
 diprings = atoi(av[i++]);
 swellrings = atoi(av[i++]);
 tiprings = atoi(av[i++]);
 dimplerings = atoi(av[i++]);
 radius = atof(av[i++]);
 length = atof(av[i++]);
 upslant = atof(av[i++]);
 stiffness = atof(av[i++]);
 headslant = atof(av[i++]);
 dip = atof(av[i++]);
 diplen1 = atof(av[i++]);
 diplen2 = atof(av[i++]);
 swell = atof(av[i++]);
 swellround = atof(av[i++]);
 tiplen = atof(av[i++]);
 tipround = atof(av[i++]);
 notchwidth = atof(av[i++]);
 dimpleoffset = atof(av[i++]);
 dimplelen = atof(av[i++]);
 dimplewidth = atof(av[i++]);
 dimpledepth = atof(av[i++]);
 check_args();
 init_arrays();
 gen_shaft();
#if 0
 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();
#endif
 return(0);
}