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

#define ARCN 15
#define SCALE .25

typedef enum {
	  PK_CORNER = 1,
	  PK_SURFACE
	  } PTKIND;

typedef struct xyz XYZ;
typedef struct tri TRI;
typedef struct pt PT;
typedef struct loc LOC;

struct xyz {
  double x;
  double y;
  double z;
  } ;

struct loc {
  LOC *link;
  XYZ p;
  int n;
  } ;

struct pt {
  PT *link;
  LOC *l;
  PTKIND kind;
  XYZ normal;
  int n;
  } ;

struct tri {
  TRI *link;
  PT *a;
  PT *b;
  PT *c;
  int n;
  } ;

#define XHP ((XYZ){1,0,0})
#define YHP ((XYZ){0,1,0})
#define ZHP ((XYZ){0,0,1})
#define XHN ((XYZ){-1,0,0})
#define YHN ((XYZ){0,-1,0})
#define ZHN ((XYZ){0,0,-1})

static LOC *locs;
static int nlocs;
static PT *pts;
static int npts;
static TRI *tris;
static int ntris;

static double xyz_length(XYZ) __attribute__((__const__));
static double xyz_length(XYZ v)
{
 return(sqrt((v.x*v.x)+(v.y*v.y)+(v.z*v.z)));
}

static XYZ xyz_smul(XYZ, double) __attribute__((__const__));
static XYZ xyz_smul(XYZ v, double s)
{
 return((XYZ){v.x*s,v.y*s,v.z*s});
}

static XYZ xyz_normalize(XYZ) __attribute__((__const__));
static XYZ xyz_normalize(XYZ v)
{
 return(xyz_smul(v,1/xyz_length(v)));
}

static XYZ xyz_add(XYZ, XYZ) __attribute__((__const__));
static XYZ xyz_add(XYZ a, XYZ b)
{
 return( (XYZ) { .x = a.x + b.x,
		 .y = a.y + b.y,
		 .z = a.z + b.z } );
}

static XYZ xyz_sub(XYZ, XYZ) __attribute__((__const__));
static XYZ xyz_sub(XYZ a, XYZ b)
{
 return( (XYZ) { .x = a.x - b.x,
		 .y = a.y - b.y,
		 .z = a.z - b.z } );
}

#if 0
static XYZ xyz_cross(XYZ, XYZ) __attribute__((__const__));
static XYZ xyz_cross(XYZ a, XYZ b)
{
 return( (XYZ) { .x = (a.y * b.z) - (a.z * b.y),
		 .y = (a.z * b.x) - (a.x * b.z),
		 .z = (a.x * b.y) - (a.y * b.x) } );
}
#endif

static LOC *loc_at(XYZ p)
{
 LOC *l;

 l = malloc(sizeof(LOC));
 l->p = p;
 l->link = locs;
 locs = l;
 nlocs ++;
 return(l);
}

static PT *pt_s(LOC *l, XYZ n)
{
 PT *v;

 v = malloc(sizeof(PT));
 v->l = l;
 v->kind = PK_SURFACE;
 v->normal = n;
 v->link = pts;
 pts = v;
 npts ++;
 return(v);
}

static PT *pt_c(LOC *l)
{
 PT *v;

 v = malloc(sizeof(PT));
 v->l = l;
 v->kind = PK_CORNER;
 v->link = pts;
 pts = v;
 npts ++;
 return(v);
}

static void add_tri(PT *a, PT *b, PT *c)
{
 TRI *t;

 (void)((volatile PT *)a)->l;
 (void)((volatile PT *)b)->l;
 (void)((volatile PT *)c)->l;
 t = malloc(sizeof(TRI));
 t->a = a;
 t->b = b;
 t->c = c;
 t->link = tris;
 tris = t;
 ntris ++;
}

static void add_prism_row(PT **pv1, PT **pv2, int nseg)
{
 for (;nseg>0;nseg--,pv1++,pv2++)
  { add_tri(pv1[0],pv1[1],pv2[0]);
    add_tri(pv1[1],pv2[1],pv2[0]);
  }
}

static void add_fan(PT *hub, PT **v, int nseg, int inc)
{
 for (;nseg>0;nseg--,v+=inc) add_tri(hub,v[0],v[inc]);
}

static void setup_mesh(void)
{
 double qc[ARCN+1][2];
 int i;
 int j;
 PT *row1[ARCN+3];
 PT *row2[ARCN+3];
 PT *corner1;
 PT *corner2;
 PT *corner3;
 PT *corner4;
 PT *midline[ARCN+1];
 PT *arc1[ARCN+1][ARCN];
 PT *arc1end[ARCN];
 PT *arc2[ARCN+1][ARCN];
 PT *arc2end[ARCN];
 PT *row3[ARCN+3];
 PT *row4[ARCN+3];
 PT *sharp[12];
 PT *row1c[ARCN+3];
 PT *row4c[ARCN+3];

 locs = 0;
 nlocs = 0;
 pts = 0;
 npts = 0;
 tris = 0;
 ntris = 0;
 for (i=0;i<=ARCN;i++)
  { qc[i][0] = cos(i*M_PI/(2*ARCN));
    qc[i][1] = sin(i*M_PI/(2*ARCN));
  }
 row1[0] = pt_s(loc_at((XYZ){1,0,1}),(XYZ){0,1,0});
 row1[ARCN+2] = pt_s(loc_at((XYZ){0,-1,1}),(XYZ){-1,0,0});
 for (i=0;i<=ARCN;i++)
  { row1[i+1] = pt_s(
	loc_at((XYZ){SCALE*(1-qc[i][1]),SCALE*(qc[i][0]-1),1}),
	(XYZ){-qc[i][1],qc[i][0],0} );
  }
 for (i=0;i<ARCN+3;i++)
  { row2[i] = pt_s(
	loc_at(xyz_sub(row1[i]->l->p,(XYZ){0,0,1})),
	row1[i]->normal );
  }
 for (i=0;i<ARCN+3;i++)
  { row3[i] = pt_s(
	loc_at((XYZ){-row2[i]->l->p.x,row2[i]->l->p.z,row2[i]->l->p.y}),
	(XYZ){-row2[i]->normal.x,row2[i]->normal.z,row2[i]->normal.y} );
  }
 for (i=0;i<ARCN+3;i++)
  { row4[i] = pt_s(
	loc_at(xyz_add(row3[i]->l->p,(XYZ){0,1,0})),
	row3[i]->normal );
  }
 midline[0] = pt_s(row3[ARCN+1]->l,(XYZ){0,1,0});
 for (i=1;i<ARCN;i++)
  { midline[i] = pt_s(
	loc_at((XYZ){0,row2[i+1]->l->p.y,-row2[i+1]->l->p.x}),
	(XYZ){0,qc[i][0],qc[i][1]} );
  }
 midline[ARCN] = pt_s(row2[ARCN+1]->l,(XYZ){0,0,1});
 for (i=0;i<ARCN;i++) arc1[0][i] = row2[i+1];
 for (i=0;i<ARCN;i++) arc1[ARCN][i] = midline[i];
 for (j=1;j<ARCN;j++)
  { for (i=0;i<ARCN;i++)
     { arc1[j][i] = pt_s(
		loc_at((XYZ){ arc1[0][i]->l->p.x * qc[j][0],
			      arc1[0][i]->l->p.y,
			      - arc1[0][i]->l->p.x * qc[j][1] }),
		(XYZ){ - qc[i][1] * qc[j][0],
		       qc[i][0],
		       qc[i][1] * qc[j][1] } );
     }
  }
 for (j=0;j<ARCN;j++)
  { arc1end[j] = pt_s(
	midline[ARCN]->l,
	xyz_normalize(xyz_add( (XYZ){ - qc[ARCN][1] * qc[j][0],
				      qc[ARCN][0],
				      qc[ARCN][1] * qc[j][1] },
			       (XYZ){ - qc[ARCN][1] * qc[j+1][0],
				      qc[ARCN][0],
				      qc[ARCN][1] * qc[j+1][1] } )) );
  }
 for (i=0;i<ARCN;i++) arc2[0][i] = row3[i+1];
 for (i=0;i<ARCN;i++) arc2[ARCN][i] = midline[ARCN-i];
 for (j=1;j<ARCN;j++)
  { for (i=0;i<ARCN;i++)
     { arc2[j][i] = pt_s(
		loc_at((XYZ){ arc2[0][i]->l->p.x * qc[j][0],
			      arc2[0][i]->l->p.x * qc[j][1],
			      arc2[0][i]->l->p.z }),
		(XYZ){ qc[i][1] * qc[j][0],
		       qc[i][1] * qc[j][1],
		       qc[i][0] } );
     }
  }
 for (j=0;j<ARCN;j++)
  { arc2end[j] = pt_s(
	midline[0]->l,
	xyz_normalize(xyz_add( (XYZ){ qc[ARCN][1] * qc[j][0],
				      qc[ARCN][1] * qc[j][1],
				      qc[ARCN][0] },
			       (XYZ){ qc[ARCN][1] * qc[j+1][0],
				      qc[ARCN][1] * qc[j+1][1],
				      qc[ARCN][0] } )) );
  }
 corner1 = pt_s(loc_at((XYZ){1,0,-1}),(XYZ){0,1,0});
 corner2 = pt_s(row3[ARCN+2]->l,(XYZ){0,1,0});
 corner3 = pt_s(loc_at((XYZ){-1,-01,0}),(XYZ){0,0,1});
 corner4 = pt_s(row2[ARCN+2]->l,(XYZ){0,0,1});
 sharp[0] = pt_c(row1[0]->l);
 sharp[1] = pt_c(row1[ARCN+2]->l);
 sharp[2] = pt_c(loc_at((XYZ){1,-1,1}));
 sharp[3] = pt_c(corner1->l);
 sharp[4] = pt_c(loc_at((XYZ){1,-1,-1}));
 sharp[5] = pt_c(corner3->l);
 sharp[6] = pt_c(loc_at((XYZ){-1,-1,-1}));
 sharp[7] = pt_c(corner2->l);
 sharp[8] = pt_c(row4[ARCN+2]->l);
 sharp[9] = pt_c(row4[0]->l);
 sharp[10] = pt_c(loc_at((XYZ){-1,1,-1}));
 sharp[11] = pt_c(corner4->l);
 add_prism_row(&row1[0],&row2[0],ARCN+2);
 add_prism_row(&row4[0],&row3[0],ARCN+2);
 add_tri(corner1,row2[0],row2[1]);
 add_tri(corner1,midline[0],corner2);
 add_tri(corner3,row3[0],row3[1]);
 add_tri(corner3,midline[ARCN],corner4);
 for (j=0;j<ARCN;j++)
  { add_prism_row(&arc1[j][0],&arc1[j+1][0],ARCN-1);
    add_tri(corner1,arc1[j][0],arc1[j+1][0]);
    add_tri(arc1end[j],arc1[j+1][ARCN-1],arc1[j][ARCN-1]);
    add_prism_row(&arc2[j][0],&arc2[j+1][0],ARCN-1);
    add_tri(corner3,arc2[j][0],arc2[j+1][0]);
    add_tri(arc2end[j],arc2[j+1][ARCN-1],arc2[j][ARCN-1]);
  }
 row1c[0] = sharp[0];
 for (i=1;i<ARCN+2;i++) row1c[i] = pt_c(row1[i]->l);
 row1c[ARCN+2] = sharp[1];
 row4c[0] = sharp[9];
 for (i=1;i<ARCN+2;i++) row4c[i] = pt_c(row4[i]->l);
 row4c[ARCN+2] = sharp[8];
 add_fan(sharp[2],&row1c[ARCN+2],ARCN+2,-1);
 add_fan(sharp[10],&row4c[ARCN+2],ARCN+2,-1);
 add_tri(sharp[0],sharp[3],sharp[4]);
 add_tri(sharp[0],sharp[4],sharp[2]);
 add_tri(sharp[4],sharp[6],sharp[5]);
 add_tri(sharp[4],sharp[5],sharp[11]);
 add_tri(sharp[4],sharp[11],sharp[1]);
 add_tri(sharp[4],sharp[1],sharp[2]);
 add_tri(sharp[6],sharp[10],sharp[9]);
 add_tri(sharp[6],sharp[9],sharp[5]);
 add_tri(sharp[6],sharp[4],sharp[3]);
 add_tri(sharp[6],sharp[3],sharp[7]);
 add_tri(sharp[6],sharp[7],sharp[8]);
 add_tri(sharp[6],sharp[8],sharp[10]);
}

static void dump_output(void)
{
 LOC *l;
 PT *p;
 TRI *t;
 int i;

 for (i=0,l=locs;l;l=l->link,i++) l->n = i;
 if (i != nlocs) abort();
 for (i=0,p=pts;p;p=p->link,i++) p->n = i;
 if (i != npts) abort();
 for (i=0,t=tris;t;t=t->link,i++) t->n = i;
 if (i != ntris) abort();
 printf("@\nModelV1\n");
 printf("1 %d %d %d\n",nlocs,npts,ntris);
 printf("White 0 0 0 .7 .7 .7 10 .5\n");
 for (l=locs;l;l=l->link) printf("%g %g %g\n",l->p.x,l->p.y,l->p.z);
 for (p=pts;p;p=p->link)
  { switch (p->kind)
     { case PK_CORNER:
	  printf("%d c\n",p->l->n);
	  break;
       case PK_SURFACE:
	  printf("%d s %g %g %g\n",p->l->n,p->normal.x,p->normal.y,p->normal.z);
	  break;
       default:
	  abort();
	  break;
     }
  }
 for (t=tris;t;t=t->link) printf("0 %d %d %d\n",t->a->n,t->b->n,t->c->n);
}

int main(void);
int main(void)
{
 setup_mesh();
 dump_output();
 return(0);
}