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

extern const char *__progname;

typedef struct pt PT;
typedef struct tri TRI;

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

struct tri {
  int px[3];
  } ;

static double cos36;
static double sin36;
static int inters;
static int npt;
static PT *pts;
static int nextpt;
static int ntri;
static TRI *tris;
static int nexttri;
static int *mesh;
static int lines[12][12][2];

static double sqr(double x)
{
 return(x*x);
}

static double dist(double d)
{
 double x1;
 double z1;
 double x2;
 double y2;
 double z2;

 x1 = sqrt(1-sqr(d));
 z1 = d;
 x2 = x1 * cos36;
 y2 = x1 * sin36;
 z2 = -d;
 return((sqr(x1)+sqr(1-z1)) - (sqr(x2-x1)+sqr(y2)+sqr(z2-z1)));
}

static double solve(double (*fn)(double), double l, double h)
{
 double m;
 double f;

 if (((*fn)(l) < 0) && ((*fn)(h) > 0))
  {
  }
 else if (((*fn)(l) > 0) && ((*fn)(h) < 0))
  { m = l;
    l = h;
    h = m;
  }
 else
  { fprintf(stderr,"invalid solve\n");
    exit(1);
  }
 while (1)
  { m = (h + l) / 2;
    if ((m == l) || (m == h)) return(m);
    f = (*fn)(m);
    if (f == 0) return(m);
    if (f < 0) l = m; else h = m;
  }
}

static PT epszero(PT p)
{
 if (fabs(p.x) < 1e-10) p.x = 0;
 if (fabs(p.y) < 1e-10) p.y = 0;
 if (fabs(p.z) < 1e-10) p.z = 0;
 return(p);
}

static void setpt(int inx, PT p)
{
 p = epszero(p);
 pts[inx] = p;
}

static double interp1(double a, double f, double b)
{
 return((a*(1-f))+(b*f));
}

static PT interp3(PT a, double f, PT b)
{
 return(MAKEPT(interp1(a.x,f,b.x),interp1(a.y,f,b.y),interp1(a.z,f,b.z)));
}

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 PT unit(PT v)
{
 return(smul3(1/norm(v),v));
}

static void makeline(int a, int b)
{
 int i;

 if (inters < 1) return;
 if (lines[a][b][0] < 0)
  { lines[a][b][0] = nextpt;
    lines[a][b][1] = 1;
    lines[b][a][0] = nextpt + inters - 1;
    lines[b][a][1] = -1;
    for (i=0;i<inters;i++)
     { pts[nextpt++] = epszero(unit(interp3(pts[a],(i+1)/(double)(inters+1),pts[b])));
     }
  }
}

static void triangle(int i1, int i2, int i3)
{
 int i;
 int j;
 int k;
 int x;

 makeline(i1,i2);
 makeline(i2,i3);
 makeline(i3,i1);
 mesh[x=0] = i1;
 j = lines[i1][i2][0];
 i = lines[i1][i2][1];
 for (k=0;k<inters;k++,j+=i) mesh[x+=inters+3] = j;
 mesh[x+=inters+3] = i2;
 j = lines[i2][i3][0];
 i = lines[i2][i3][1];
 for (k=0;k<inters;k++,j+=i) mesh[--x] = j;
 mesh[--x] = i3;
 j = lines[i3][i1][0];
 i = lines[i3][i1][1];
 for (k=0;k<inters;k++,j+=i) mesh[x-=inters+2] = j;
 if (inters >= 2)
  { for (i=0;i<inters-1;i++)
     { for (j=0;j<=i;j++)
	{ mesh[((i+2)*(inters+2))+j+1] = nextpt;
	  pts[nextpt++] = epszero(unit(interp3(interp3(pts[i1],
						       (i+2)/(double)(inters+1),
						       pts[i3]),
					       (j+1)/(double)(i+2),
					       interp3(pts[i1],
						       (i+2)/(double)(inters+1),
						       pts[i2]))));
	}
     }
  }
 for (i=0;i<=inters;i++)
  { for (j=0;j<=i;j++)
     { tris[nexttri++] = (TRI){ .px = { mesh[(i*(inters+2))+j],
					mesh[((i+1)*(inters+2))+j+1],
					mesh[((i+1)*(inters+2))+j] } };
     }
    for (j=0;j<i;j++)
     { tris[nexttri++] = (TRI){ .px = { mesh[(i*(inters+2))+j],
					mesh[(i*(inters+2))+j+1],
					mesh[((i+1)*(inters+2))+j+1] } };
     }
  }
}

static void dump(void)
{
 int i;

 printf("ModelV1\n");
 printf("1 %d %d %d\n",npt,npt,ntri);
 printf("White 0 0 0 .7 .7 .7 10 .5\n");
 for (i=0;i<npt;i++) printf("%g %g %g\n",pts[i].x,pts[i].y,pts[i].z);
 for (i=0;i<npt;i++) printf("%d s %g %g %g\n",i,pts[i].x,pts[i].y,pts[i].z);
 for (i=0;i<ntri;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 d;
 double r;
 int i;
 int j;

 if (ac != 2)
  { fprintf(stderr,"Usage: %s intermediate-point-count\n",__progname);
    exit(1);
  }
 inters = atoi(av[1]);
 npt = 12 + (20 * inters) + (10 * inters * inters);
 pts = malloc(npt*sizeof(PT));
 ntri = 20 + (40 * inters) + (20 * inters * inters);
 tris = malloc(ntri*sizeof(TRI));
 mesh = malloc((inters+2)*(inters+2)*sizeof(int));
 cos36 = cos(M_PI/5);
 sin36 = sin(M_PI/5);
 /* This is grossly inelegant; there certainly ought to be a closed-form
    solution for d, and indeed I seem to recall working it out once.  But
    since I seem to be unable to do it again.... */
 d = solve(&dist,0,1);
 r = sqrt(1-sqr(d));
 setpt(0,MAKEPT(0,0,1));
 setpt(11,MAKEPT(0,0,-1));
 for (i=0;i<5;i++)
  { setpt(i+1,MAKEPT(r*cos(i*2*M_PI/5),r*sin(i*2*M_PI/5),d));
    setpt(i+6,MAKEPT(r*cos((i+i+1)*M_PI/5),r*sin((i+i+1)*M_PI/5),-d));
  }
 for (i=0;i<12;i++) for (j=0;j<12;j++) lines[i][j][0] = -1;
 nextpt = 12;
 for (i=0;i<4;i++)
  { triangle(0,i+2,i+1);
    triangle(11,i+6,i+7);
    triangle(i+1,i+2,i+6);
    triangle(i+2,i+7,i+6);
  }
 triangle(0,1,5);
 triangle(11,10,6);
 triangle(5,1,10);
 triangle(1,6,10);
 if ((nextpt != npt) || (nexttri != ntri)) abort();
 dump();
 return(0);
}