#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 double angle;
static double cangle;
static double sangle;
static double tangle2;
static int circpts;
static int qcircpts;
static int bendpts;

#define GETPT(n) ((PT3 *)allocv_el(pts,(n)))
#define GETPOINT(n) ((POINT *)allocv_el(points,(n)))
#define GETTRI(n) ((TRI *)allocv_el(tris,(n)))

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;

 if ((a == b) || (b == c) || (c == a)) abort();
 if (GETPOINT(a)->pt == GETPOINT(b)->pt) abort();
 if (GETPOINT(b)->pt == GETPOINT(c)->pt) abort();
 if (GETPOINT(c)->pt == GETPOINT(a)->pt) abort();
 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 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 void gen(void)
{
 int i;
 int i2q;
 int j;
 double a;
 double ca;
 double sa;
 PT3 n0;
 PT3 n1;
 PT3 offset;
 PT3 n;
 int do_v;
 int cpt[2];
 int cpoint[2];
 int endpt[2][circpts];
 int endpoint[2][circpts];
 int inpt[2][circpts];
 int inpoint[2][circpts];
 int bendpt[(2*qcircpts)-1][bendpts-1];
 int bendpoint[(2*qcircpts)-1][bendpts+1];
 double cosv[circpts];
 double sinv[circpts];
 int veept[(2*qcircpts)-1];
 int veepoint[2][(2*qcircpts)-1];

 do_v = (angle > M_PI/2);
 for (i=circpts-1;i>=0;i--)
  { a = (i * M_PI * 2) / circpts;
    cosv[i] = sin(a); /* cos(a-90°) */
    sinv[i] = - cos(a); /* sin(a-90°) */
  }
 n0 = MAKEPT3(-1,0,0);
 n1 = MAKEPT3(cangle,sangle,0);
 cpt[0] = newpt(MAKEPT3(-10,0,0));
 cpoint[0] = newpoint(cpt[0],n0);
 cpt[1] = newpt(MAKEPT3(10*cangle,10*sangle,0));
 cpoint[1] = newpoint(cpt[1],n1);
 for (i=circpts-1;i>=0;i--)
  { endpt[0][i] = newpt(MAKEPT3(-10,cosv[i],sinv[i]));
    endpoint[0][i] = newpoint(endpt[0][i],n0);
    endpt[1][i] = newpt(MAKEPT3( (10*cangle) - (sangle*cosv[i]),
				 (10*sangle) + (cangle*cosv[i]),
				 sinv[i] ));
    endpoint[1][i] = newpoint(endpt[1][i],n1);
  }
 cur_material = MAT_WHITE;
 newtri(endpoint[0][circpts-1],endpoint[0][0],cpoint[0]);
 for (i=circpts-2;i>=0;i--) newtri(endpoint[0][i],endpoint[0][i+1],cpoint[0]);
 newtri(endpoint[1][0],endpoint[1][circpts-1],cpoint[1]);
 for (i=circpts-2;i>=0;i--) newtri(endpoint[1][i+1],endpoint[1][i],cpoint[1]);
 for (i=(2*qcircpts)-1;i>=0;i--)
  { inpt[0][i] = newpt(MAKEPT3(-cosv[i]*tangle2,cosv[i],sinv[i]));
    inpt[1][i] = inpt[0][i];
  }
 j = 2 * qcircpts;
 offset = smul3(-10,n0);
 for (i=circpts-1;i>=j;i--)
  { inpt[0][i] = newpt(add3(offset,*GETPT(endpt[0][i])));
  }
 offset = smul3(-10,n1);
 for (i=circpts-1;i>=j;i--)
  { inpt[1][i] = newpt(add3(offset,*GETPT(endpt[1][i])));
  }
 if (do_v)
  { for (i=(2*qcircpts)-1;i>0;i--) veept[i-1] = inpt[0][i];
    offset = smul3(.1,n0);
    for (i=(2*qcircpts)-1;i>0;i--)
     { inpt[0][i] = newpt(add3(*GETPT(inpt[0][i]),offset));
     }
    offset = smul3(.1,n1);
    for (i=(2*qcircpts)-1;i>0;i--)
     { inpt[1][i] = newpt(add3(*GETPT(inpt[1][i]),offset));
     }
  }
 for (i=circpts-1;i>=0;i--)
  { n = MAKEPT3(0,cosv[i],sinv[i]);
    endpoint[0][i] = newpoint(endpt[0][i],n);
    inpoint[0][i] = newpoint(inpt[0][i],n);
    if (do_v && (i < 2*qcircpts) && (i > 0))
     { veepoint[0][i-1] = newpoint(veept[i-1],n);
     }
  }
 for (i=circpts-1;i>=0;i--)
  { n = MAKEPT3(-sangle*cosv[i],cangle*cosv[i],sinv[i]);
    endpoint[1][i] = newpoint(endpt[1][i],n);
    inpoint[1][i] = newpoint(inpt[1][i],n);
    if (do_v && (i < 2*qcircpts) && (i > 0))
     { veepoint[1][i-1] = newpoint(veept[i-1],n);
     }
  }
 for (i=(2*qcircpts)-2;i>=0;i--)
  { i2q = i + 1 + (2 * qcircpts);
    for (j=bendpts-2;j>=0;j--)
     { a = angle * (j+1) / bendpts;
       ca = cos(a);
       sa = sin(a);
       n = MAKEPT3(-sa*cosv[i2q],ca*cosv[i2q],sinv[i2q]);
       bendpoint[i][j+1] = newpoint(bendpt[i][j]=newpt(n),n);
     }
    bendpoint[i][0] = inpoint[0][i2q];
    bendpoint[i][bendpts] = inpoint[1][i2q];
  }
 newtri(endpoint[0][0],endpoint[0][circpts-1],inpoint[0][0]);
 newtri(endpoint[0][circpts-1],inpoint[0][circpts-1],inpoint[0][0]);
 for (i=circpts-2;i>=0;i--)
  { newtri(endpoint[0][i+1],endpoint[0][i],inpoint[0][i+1]);
    newtri(endpoint[0][i],inpoint[0][i],inpoint[0][i+1]);
  }
 newtri(endpoint[1][0],inpoint[1][0],endpoint[1][circpts-1]);
 newtri(endpoint[1][circpts-1],inpoint[1][0],inpoint[1][circpts-1]);
 for (i=circpts-2;i>=0;i--)
  { newtri(endpoint[1][i+1],inpoint[1][i+1],endpoint[1][i]);
    newtri(endpoint[1][i],inpoint[1][i+1],inpoint[1][i]);
  }
 i = (2 * qcircpts) - 2;
 for (j=bendpts-1;j>=0;j--)
  { newtri(inpoint[0][2*qcircpts],bendpoint[0][j+1],bendpoint[0][j]);
    newtri(inpoint[0][0],bendpoint[i][j],bendpoint[i][j+1]);
  }
 for (i--;i>=0;i--)
  { for (j=bendpts-1;j>=0;j--)
     { newtri(bendpoint[i][j+1],bendpoint[i+1][j+1],bendpoint[i][j]);
       newtri(bendpoint[i][j],bendpoint[i+1][j+1],bendpoint[i+1][j]);
     }
  }
 if (do_v)
  { cur_material = MAT_DARK;
    i = 2 * qcircpts;
    newtri(inpoint[0][0],veepoint[0][0],inpoint[0][1]);
    newtri(inpoint[1][0],inpoint[1][1],veepoint[1][0]);
    newtri(inpoint[0][i],inpoint[0][i-1],veepoint[0][i-2]);
    newtri(inpoint[1][i],veepoint[1][i-2],inpoint[1][i-1]);
    for (i-=2;i>0;i--)
     { newtri(veepoint[0][i-1],inpoint[0][i+1],inpoint[0][i]);
       newtri(veepoint[0][i-1],veepoint[0][i],inpoint[0][i+1]);
       newtri(veepoint[1][i-1],inpoint[1][i],inpoint[1][i+1]);
       newtri(veepoint[1][i-1],inpoint[1][i+1],veepoint[1][i]);
     }
  }
}

int main(int, char **);
int main(int ac, char **av)
{
 if (ac != 4)
  { fprintf(stderr,"Usage: %s angle quartercircpts bendpts\n",__progname);
    exit(1);
  }
 angle = DEGTORAD(atof(av[1]));
 cangle = cos(angle);
 sangle = sin(angle);
 tangle2 = tan(angle/2);
 qcircpts = atoi(av[2]);
 circpts = qcircpts * 4;
 bendpts = atoi(av[3]);
 init_arrays();
 gen();
 dump();
 return(0);
}