/*
 * Args, with example values:
 *	50	nring (count, triangles per ring)
 *	20	shaftrings (count, rings in main shaft)
 *	3	diprings (count, rings in dip half)
 *	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	limpness (0=straight, 1=limp)
 *	30	headslant (degrees, 0=square across)
 *	5	dip (just below head, percent of radius)
 *	.15	diplen1 (distance, normal shaft to min radius)
 *	.05	diplen2 (distance, swell edge 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 "2d.h"
#include "3d.h"
#include "allocv.h"

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

typedef struct ptn PTN;
typedef struct tri TRI;
typedef struct point POINT;

struct ptn {
  PT3 pt;
  PT3 n;
  } ;

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 PTN **ptmtx;
static int *ring0;

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

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

static int cur_material;

static double sindeg(double a)
{
 return(sin(DEGTORAD(a)));
}

static double cosdeg(double a)
{
 return(cos(DEGTORAD(a)));
}

static double tandeg(double a)
{
 return(tan(DEGTORAD(a)));
}

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 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 ((limpness < 0) || (limpness > 1)) unreasonable("limpness 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)
{
 int tr;
 int i;

 pts = allocv_init(sizeof(PT3));
 points = allocv_init(sizeof(POINT));
 tris = allocv_init(sizeof(TRI));
 tr = shaftrings + (2 * diprings) + swellrings + tiprings + dimplerings + 1;
 ptmtx = malloc(tr*sizeof(*ptmtx));
 ptmtx[0] = malloc((tr*nring)*sizeof(**ptmtx));
 for (i=1;i<tr;i++) ptmtx[i] = ptmtx[i-1] + nring;
 ring0 = malloc(nring*sizeof(*ring0));
}

static PTN profile_pt(PT2 (*profile)(int, double), int ringno, int odd, int circum)
{
 PT2 pp;
 double ca;
 double cc;
 double cs;
 double sla;
 PT3 rc;
 PT3 x;
 PT3 y;
 PT3 z;
 PT3 r;
 double maxsl;

 ca = (((circum * 2) + (odd?1:0)) * M_PI) / nring;
 pp = (*profile)(ringno,ca);
 cs = sin(ca);
 cc = cos(ca);
 if (limpness < 1e-3)
  { z = MAKEPT3(cosdeg(upslant),0,sindeg(upslant));
    rc = smul3(pp.x,z);
  }
 else
  { double radcurv;
    radcurv = radius / limpness;
    maxsl = DEGTORAD(upslant+90) * radcurv;
    if (pp.x < maxsl)
     { sla = pp.x / radcurv;
       rc = MAKEPT3( (sindeg(upslant)-sin(DEGTORAD(upslant)-sla)) * radcurv,
		     0,
		     (cos(DEGTORAD(upslant)-sla)-cosdeg(upslant)) * radcurv );
       z = MAKEPT3(cos(DEGTORAD(upslant)-sla),0,sin(DEGTORAD(upslant)-sla));
     }
    else
     { rc = MAKEPT3( (sindeg(upslant)+1) * radcurv,
		     0,
		     (-cosdeg(upslant) * radcurv) - (pp.x - maxsl) );
       z = MAKEPT3(0,0,-1);
     }
  }
 y = MAKEPT3(0,-1,0);
 x = cross(y,z);
 r = add3(smul3(cc,x),smul3(cs,y));
 return((PTN){.pt=add3(rc,smul3(pp.y,r)),.n=r});
}

static PT2 profile_shaft(int ringno, double ca)
{
 return( (PT2){ .x = (length + (radius * (1 - cos(ca)) * tandeg(headslant))) * (ringno / (double)shaftrings),
		.y = 1 } );
}

static double notchfactor(double ca)
{
 double d;

 d = (ca - M_PI) / DEGTORAD(notchwidth);
 return( (fabs(d) < 1) ? (1-cos(d*M_PI))/2 : 1 );
}

static PT2 profile_dip1(int ringno, double ca)
{
 PT2 s;
 double nf;
 double rf;

 s = profile_shaft(shaftrings,ca);
 nf = notchfactor(ca);
 rf = ringno / (double)diprings;
 s.x += diplen1 * rf;
 s.y = 1 - (((1 - cos(rf*M_PI)) / 2) * nf * dip * .01);
 return(s);
}

static PT2 profile_dip2(int ringno, double ca)
{
 PT2 s;
 double nf;
 double rf;

 s = profile_shaft(shaftrings,ca);
 nf = notchfactor(ca);
 rf = (ringno * M_PI) / (2.0 * diprings);
 s.x += diplen1 + (diplen2 * sin(rf));
 s.y = 1 - (cos(rf) * nf * dip * .01);
 return(s);
}

static PTN shaft_pt(int ringno, int odd, int circum)
{
 return(profile_pt(&profile_shaft,ringno,odd,circum));
}

static PTN dip1_pt(int ringno, int odd, int circum)
{
 return(profile_pt(&profile_dip1,ringno,odd,circum));
}

static PTN dip2_pt(int ringno, int odd, int circum)
{
 return(profile_pt(&profile_dip2,ringno,odd,circum));
}

static void gen_base(void)
{
 PT3 bn;
 int op;
 int *rp;
 int i;

 bn = MAKEPT3(-cosdeg(upslant),0,-sindeg(upslant));
 op = newpoint(newpt(MAKEPT3(0,0,0)),bn);
 rp = malloc((nring+1)*sizeof(*rp));
 for (i=nring-1;i>=0;i--) rp[i] = newpoint(ring0[i],bn);
 rp[nring] = rp[0];
 for (i=nring;i>0;i--) newtri(op,rp[i-1],rp[i]);
 free(rp);
}

static void gen_rings(int ringcount)
{
 int *r1;
 int *r2;
 int *rc;
 int *rp;
 int *rt;
 int r;
 int c;

 r1 = malloc((nring+1)*sizeof(*r1));
 r2 = malloc((nring+1)*sizeof(*r2));
 rc = r1;
 rp = r2;
 for (c=nring-1;c>=0;c--)
  { rc[c] = newpoint(ring0[c],ptmtx[r][c].n);
  }
 rc[nring] = rc[0];
 for (r=1;r<ringcount;r++)
  { rt = rc;
    rc = rp;
    rp = rt;
    for (c=nring-1;c>=0;c--)
     { rc[c] = newpoint(newpt(ptmtx[r][c].pt),ptmtx[r][c].n);
     }
    rc[nring] = rc[0];
    if (r & 1)
     { for (c=nring-1;c>=0;c--)
	{ newtri(rp[c],rc[c],rp[c+1]);
	  newtri(rp[c+1],rc[c],rc[c+1]);
	}
     }
    else
     { for (c=nring-1;c>=0;c--)
	{ newtri(rp[c],rc[c],rc[c+1]);
	  newtri(rp[c],rc[c+1],rp[c+1]);
	}
     }
  }
}

static void gen(void)
{
 int r;
 int i;
 int c;

 r = 0;
 for (i=0;i<=shaftrings;i++)
  { for (c=nring-1;c>=0;c--)
     { ptmtx[r][c] = shaft_pt(i,r&1,c);
     }
    r ++;
  }
 for (i=1;i<=diprings;i++)
  { for (c=nring-1;c>=0;c--)
     { ptmtx[r][c] = dip1_pt(i,r&1,c);
     }
    r ++;
  }
 for (i=1;i<=diprings;i++)
  { for (c=nring-1;c>=0;c--)
     { ptmtx[r][c] = dip2_pt(i,r&1,c);
     }
    r ++;
  }
 for (c=nring-1;c>=0;c--) ring0[c] = newpt(ptmtx[0][c].pt);
 gen_base();
 gen_rings(r);
}

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

int main(int, char **);
int main(int ac, char **av)
{
 int i;

 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 half)\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\
 limpness (0=straight, 1=limp)\n\
 headslant (degrees, 0=straight across)\n\
 dip (below head, percent of radius)\n\
 diplen1 (distance, shaft to min radius)\n\
 diplen2 (distance, min radius to swell edge)\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++]);
 limpness = 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();
 dump();
 return(0);
}