#define NRING 25

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

extern const char *__progname;

typedef struct pt PT;
typedef struct tri TRI;
typedef struct point POINT;

struct pt {
  double x;
  double y;
  double z;
  } ;
#define MAKEPT(xv,yv,zv) ((PT){.x=(xv),.y=(yv),.z=(zv)})

struct point {
  int pt;
  PT normal;
  } ;

struct tri {
  int px[3];
  } ;

static int ring1[NRING+1];
static int ring2[NRING+1];
static int *rcur;
static int *rprev;

static PT pts[35*NRING];
static int npts;
static POINT points[40*NRING];
static int npoints;
static TRI tris[80*NRING];
static int ntris;

static PT smul3(double f, PT p)
{
 return(MAKEPT(p.x*f,p.y*f,p.z*f));
}

static double norm(PT p)
{
 return(sqrt((p.x*p.x)+(p.y*p.y)+(p.z*p.z)));
}

static double dot(PT a, PT b)
{
 return((a.x*b.x)+(a.y*b.y)+(a.z*b.z));
}

static PT cross(PT a, PT b)
{
 return(MAKEPT((a.y*b.z)-(a.z*b.y),(a.z*b.x)-(a.x*b.z),(a.x*b.y)-(a.y*b.x)));
}

static PT sub3(PT a, PT b)
{
 return(MAKEPT(a.x-b.x,a.y-b.y,a.z-b.z));
}

static PT add3(PT a, PT b)
{
 return(MAKEPT(a.x+b.x,a.y+b.y,a.z+b.z));
}

static PT unit(PT v)
{
 return(smul3(1/norm(v),v));
}

static int newpt(PT loc)
{
 int x;

 x = npts ++;
 if (x >= sizeof(pts)/sizeof(pts[0])) abort();
 pts[x] = loc;
 return(x);
}

static int newpoint(int loc, PT norm)
{
 int x;

 x = npoints ++;
 if (x >= sizeof(points)/sizeof(points[0])) abort();
 if ((loc < 0) || (loc >= npts)) abort();
 points[x] = (POINT){ .pt=loc, .normal=unit(norm) };
 return(x);
}

static int newtri(int a, int b, int c)
{
 int x;

 x = ntris ++;
 if (x >= sizeof(tris)/sizeof(tris[0])) abort();
 if ( (a < 0) || (a >= npoints) ||
      (b < 0) || (b >= npoints) ||
      (c < 0) || (c >= npoints) ) abort();
 tris[x] = (TRI){ .px={ [0]=a, [1]=b, [2]=c } };
 return(x);
}

static void initring(PT cen, PT axis, PT side, double rad)
{
 int i;
 int cx;
 PT r[NRING];
 int pv[NRING];
 PT norm;
 double c;
 double s;
 PT rx;
 PT ry;

 norm = smul3(-1,axis);
 axis = unit(axis);
 rx = unit(sub3(side,smul3(dot(side,axis),axis)));
 ry = unit(cross(axis,rx));
 cx = newpoint(newpt(cen),norm);
 for (i=NRING-1;i>=0;i--)
  { c = cos((i*2*M_PI)/NRING);
    s = sin((i*2*M_PI)/NRING);
    r[i] = add3(smul3(c*rad,rx),smul3(s*rad,ry));
    pv[i] = newpt(add3(cen,r[i]));
    ring1[i] = newpoint(pv[i],norm);
  }
 ring1[NRING] = ring1[0];
 for (i=NRING;i>0;i--) newtri(cx,ring1[i-1],ring1[i]);
 for (i=NRING-1;i>=0;i--) ring1[i] = newpoint(pv[i],r[i]);
 ring1[NRING] = ring1[0];
 rcur = &ring1[0];
 rprev = &ring2[0];
}

static void swaprings(void)
{
 int *t;

 t = rcur;
 rcur = rprev;
 rprev = t;
}

static void ring(PT cen, PT axis, PT side, double rad, int offset)
{
 int i;
 double c;
 double s;
 PT rx;
 PT ry;
 PT r;

 if ((offset < 0) || (offset > 1)) abort();
 axis = unit(axis);
 rx = unit(sub3(side,smul3(dot(side,axis),axis)));
 ry = unit(cross(axis,rx));
 swaprings();
 for (i=NRING-1;i>=0;i--)
  { c = cos(((i+i+offset)*M_PI)/NRING);
    s = sin(((i+i+offset)*M_PI)/NRING);
    r = add3(smul3(c*rad,rx),smul3(s*rad,ry));
    rcur[i] = newpoint(newpt(add3(cen,r)),r);
  }
 rcur[NRING] = rcur[0];
 for (i=NRING-1;i>=0;i--)
  { newtri(rcur[i],rcur[i+1],rprev[i+offset]);
    newtri(rcur[i+(offset^1)],rprev[i+1],rprev[i]);
  }
}

static void endcap(PT cen, PT axis)
{
 int i;
 int cx;

 for (i=NRING-1;i>=0;i--) rcur[i] = newpoint(points[rcur[i]].pt,axis);
 rcur[NRING] = rcur[0];
 cx = newpoint(newpt(cen),axis);
 for (i=NRING;i>0;i--) newtri(cx,rcur[i],rcur[i-1]);
}

static void dump(void)
{
 int i;

 printf("ModelV1\n");
 printf("1 %d %d %d\n",npts,npoints,ntris);
 printf("White 0 0 0 .7 .7 .7 10 .5\n");
 for (i=0;i<npts;i++) printf("%g %g %g\n",pts[i].x,pts[i].y,pts[i].z);
 for (i=0;i<npoints;i++) printf("%d s %g %g %g\n", points[i].pt,
				points[i].normal.x,
				points[i].normal.y,
				points[i].normal.z);
 for (i=0;i<ntris;i++) printf("0 %d %d %d\n",tris[i].px[0],tris[i].px[1],tris[i].px[2]);
}

int main(int, char **);
int main(int ac, char **av)
{
 double a;
 double r;
 double sa;
 double ca;
 int i;

 if (ac != 3)
  { fprintf(stderr,"Usage: %s angle-in-degrees tube-radius\n",__progname);
    exit(1);
  }
 a = atof(av[1]) * M_PI / 180;
 r = atof(av[2]);
 npts = 0;
 npoints = 0;
 ntris = 0;
 initring(MAKEPT(-r*2.5,-1,0),MAKEPT(1,0,0),MAKEPT(0,-1,0),r);
 for (i=-9;i<0;i++) ring(MAKEPT(i*r/4.0,-1,0),MAKEPT(1,0,0),MAKEPT(0,-1,0),r,(-i)&1);
 for (i=0;i<=10;i++)
  { sa = sin((a*i)/10);
    ca = cos((a*i)/10);
    ring(MAKEPT(sa,-ca,0),MAKEPT(ca,sa,0),MAKEPT(sa,-ca,0),r,i&1);
  }
 for (i=1;i<=10;i++)
  { ring(MAKEPT(sa+(i*ca*r/4.0),-ca+(i*sa*r/4.0),0),MAKEPT(ca,sa,0),MAKEPT(sa,-ca,0),r,i&1);
  }
 endcap(MAKEPT(sa+(ca*r*2.5),-ca+(sa*r*2.5),0),MAKEPT(ca,sa,0));
 dump();
 return(0);
}