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

extern const char *__progname;

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

typedef struct tri TRI;
typedef struct point POINT;

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 ALLOCV *pts;
static ALLOCV *points;
static ALLOCV *tris;

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=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 gen_cylinder(
	PT3 axis1,
	PT3 axis2,
	PT3 (*ptfxn)(double, double),
	PT3 rdir1,
	double rad1,
	double rad2,
	int sides,
	int rings )
{
 PT3 x;
 PT3 y;
 PT3 z;
 PT3 zhat;
 int r;
 int rx;
 int c1;
 int c2;
 PT3 n1;
 PT3 n2;
 int end1[sides];
 int end2[sides];
 int vp3[rings+1][sides];
 int vpt[rings+1][sides];

 static PT3 place_pt(int ring, double angle)
  { PT3 p;
    p = (*ptfxn)(ring/(double)rings,angle);
    return( add3( add3( axis1,
			smul3(p.z,z) ),
		  add3( smul3(p.x,x),
			smul3(p.y,y) ) ) );
  }

 static int t(int a, int b, int c)
  { PT3 n;
    POINT *pa;
    POINT *pb;
    POINT *pc;
    PT3 *p3a;
    PT3 *p3b;
    PT3 *p3c;
    pa = allocv_el(points,a);
    pb = allocv_el(points,b);
    pc = allocv_el(points,c);
    p3a = allocv_el(pts,pa->pt);
    p3b = allocv_el(pts,pb->pt);
    p3c = allocv_el(pts,pc->pt);
    n = unit3(cross(sub3(*p3c,*p3a),sub3(*p3b,*p3a)));
    if ((n.x != n.x) || (n.y != n.y) || (n.z != n.z)) abort();
    pa->normal = add3(pa->normal,n);
    pb->normal = add3(pb->normal,n);
    pc->normal = add3(pc->normal,n);
    return(newtri(a,b,c));
  }

 z = sub3(axis2,axis1);
 zhat = unit3(z);
 x = smul3(rad1,unit3(sub3(rdir1,smul3(dot3(rdir1,zhat),zhat))));
 y = smul3(rad2,unit3(cross(z,x)));
 for (r=rings;r>=0;r--)
  { for (rx=sides-1;rx>=0;rx--)
     { vp3[r][rx] = newpt(place_pt(r,((rx+(r*.5))*M_PI*2)/sides));
       vpt[r][rx] = newpoint(vp3[r][rx],MAKEPT3(0,0,0));
     }
  }
 n1 = unit3(sub3(axis1,axis2));
 n2 = smul3(-1,n1);
 c1 = newpoint(newpt(axis1),n1);
 c2 = newpoint(newpt(axis2),n2);
 for (rx=sides-1;rx>=0;rx--)
  { end1[rx] = newpoint(vp3[0][rx],n1);
    end2[rx] = newpoint(vp3[rings][rx],n2);
  }
 t(end1[sides-1],end1[0],c1);
 t(end2[0],end2[sides-1],c2);
 for (rx=sides-2;rx>=0;rx--)
  { t(end1[rx],end1[rx+1],c1);
    t(end2[rx+1],end2[rx],c2);
  }
 for (r=rings-1;r>=0;r--)
  { cur_material = MAT_WHITE;
    t(vpt[r][sides-1],vpt[r+1][sides-1],vpt[r][0]);
    cur_material = MAT_WHITE;
    t(vpt[r+1][sides-1],vpt[r+1][0],vpt[r][0]);
    for (rx=sides-1;rx>=1;rx--)
     { cur_material = MAT_WHITE;
       t(vpt[r][rx-1],vpt[r+1][rx-1],vpt[r][rx]);
       cur_material = MAT_WHITE;
       t(vpt[r+1][rx-1],vpt[r+1][rx],vpt[r][rx]);
     }
  }
}

static void normalize_normals(void)
{
 int i;

 for (i=allocv_n(points)-1;i>=0;i--)
  { POINT *e;
    e = allocv_el(points,i);
    e->normal = unit3(e->normal);
  }
}

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

static PT3 cyl_torsoblend(double along, double angle, double offset)
{
 PT3 p1;
 PT3 p2;

 p1 = MAKEPT3(cos(angle),sin(angle),0);
 if (cos(angle) > 0) p1.x += offset; else p1.x -= offset;
 p2 = MAKEPT3(cos(angle)*(offset+1),sin(angle),1);
 return(interp3(p1,along,p2));
}

#define LEG_SIDES 17
#define LEG_RINGS 8
#define LEG_LENGTH 5
#define LEG_OFFSET 1.5

#define CROTCH_ZSCALE 2
#define CROTCH_SIDES 17
#define CROTCH_RINGS 8
#define CROTCH_BARREL 1

#define TORSO_BLEND 3
#define TORSO_SIDES 17
#define TORSO_RINGS 6

int main(void);
int main(void)
{
 init_arrays();
 gen_cylinder( MAKEPT3(-LEG_OFFSET,0,-LEG_LENGTH/2.0),
	       MAKEPT3(-LEG_OFFSET,0,LEG_LENGTH/2.0),
	       ({ static PT3 foo(double along, double angle)
		   { return(MAKEPT3(cos(angle),sin(angle),along));
		   }
		  &foo;
		}),
	       MAKEPT3(0,1,0), 1, 1,
	       LEG_SIDES, LEG_RINGS );
 gen_cylinder( MAKEPT3(LEG_OFFSET,0,-LEG_LENGTH/2.0),
	       MAKEPT3(LEG_OFFSET,0,LEG_LENGTH/2.0),
	       ({ static PT3 foo(double along, double angle)
		   { return(MAKEPT3(cos(angle),sin(angle),along));
		   }
		  &foo;
		}),
	       MAKEPT3(0,1,0), 1, 1,
	       LEG_SIDES, LEG_RINGS );
 gen_cylinder( MAKEPT3(-LEG_OFFSET,0,LEG_LENGTH/2.0),
	       MAKEPT3(LEG_OFFSET,0,LEG_LENGTH/2.0),
	       ({ static PT3 foo(double along, double angle)
		   { double bf;
		     bf = (along - .5) * (along - .5) * CROTCH_BARREL;
		     return(MAKEPT3(cos(angle)*(1-bf),sin(angle),along));
		   }
		  &foo;
		}),
	       MAKEPT3(0,0,1), 2, 1,
	       CROTCH_SIDES, CROTCH_RINGS );
 gen_cylinder( MAKEPT3(0,0,LEG_LENGTH/2.0),
	       MAKEPT3(0,0,(LEG_LENGTH/2.0)+TORSO_BLEND),
	       ({ static PT3 foo(double along, double angle)
		   { return(cyl_torsoblend(along,angle,LEG_OFFSET));
		   }
		  &foo;
		}),
	       MAKEPT3(1,0,0), 1, 1,
	       TORSO_SIDES, TORSO_RINGS );
 normalize_normals();
 dump();
 return(0);
}