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

#define N 20

#include "3darith.h"

typedef struct pt PT;

struct pt {
  XYZ p;
  int x;
  int nl;
  int edge;
#define MAXLINKS 6
  PT *links[MAXLINKS];
  } ;

#define NPT (((N+1)*(N+2))/2)

static PT p[NPT];

static void find_p(void)
{
 double arcc[N+1];
 double arcs[N+1];
 PT *v[N+2];
 int i;
 int j;
 int ptx;
 PT *pp;

 void link_between(PT *a, PT *b)
  { if (!a || !b) return;
    if (a->nl >= MAXLINKS) abort();
    a->links[a->nl++] = b;
    if (b->nl >= MAXLINKS) abort();
    b->links[b->nl++] = a;
  }

 static void draw_projection(
	double cx, double cy, double scale,
	XYZ cp, XYZ eyedir, XYZ up )
  { int x;
    int l;
    PT *pp;
    XYZ rt;
    void project(XYZ pt)
     { pt = xyzsub(pt,cp);
       printf("%g %g",cx+(scale*dot(pt,rt)),cy+(scale*dot(pt,up)));
     }
    printf("gsave\n");
    eyedir = xyzunit(eyedir);
    up = xyzunit(subtract_component(up,eyedir));
    rt = xyzcross(up,eyedir);
    for (x=NPT-1;x>=0;x--)
     { pp = &p[x];
       for (l=pp->nl-1;l>=0;l--)
	{ project(pp->p);
	  printf(" moveto ");
	  project(pp->links[l]->p);
	  printf(" lineto\n");
	}
     }
    printf(".1 setlinewidth stroke grestore\n");
  }

 void dump(void)
  { draw_projection( 306, 300, 275,
		     (XYZ){1/3.,1/3.,1/3.},
		     (XYZ){1,0,0},
		     (XYZ){0,0,1} );
    draw_projection( 306, 600, 275,
		     (XYZ){1/3.,1/3.,1/3.},
		     (XYZ){0,1,0},
		     (XYZ){0,0,1} );
    printf("showpage\n");
  }

 int adjust(void)
  { PT *ap;
    int x;
    XYZ s;
    int i;
    double d;
    double maxd;
    int rv;
    PT *lp;
    rv = 0;
    for (x=NPT-1;x>=0;x--)
     { ap = &p[x];
       if (ap->edge) continue;
       s = (XYZ){0,0,0};
       maxd = 0;
       for (i=ap->nl-1;i>=0;i--)
	{ lp = ap->links[i];
	  s = xyzadd(s,lp->p);
	  d = xyzlength(xyzsub(lp->p,ap->p));
	  if (d > maxd) maxd = d;
	}
       s = xyzunit(xyzscale(s,1./ap->nl));
       d = xyzlength(xyzsub(ap->p,s));
       if (d > maxd/100) rv = 1;
       ap->p = s;
     }
    return(rv);
  }

 for (i=N;i>=0;i--)
  { arcc[i] = cos(i*M_PI/(N*2));
    arcs[i] = sin(i*M_PI/(N*2));
  }
 for (i=NPT-1;i>=0;i--) p[i].x = i;
 ptx = 0;
 for (i=N;i>=0;i--)
  { pp = &p[ptx++];
    pp->p.x = arcc[i];
    pp->p.y = 0;
    pp->p.z = arcs[i];
    pp->nl = 0;
    pp->edge = 1;
    v[i] = pp;
    if (i < N) link_between(pp,v[i+1]);
  }
 v[N+1] = 0;
 for (j=N-1;j>=0;j--)
  { for (i=j;i>=0;i--)
     { pp = &p[ptx++];
       if (i == 0)
	{ pp->p.x = arcc[N-j];
	  pp->p.y = arcs[N-j];
	  pp->p.z = 0;
	  pp->edge = 1;
	}
       else if (i == j)
	{ pp->p.x = 0;
	  pp->p.y = arcs[N-j];
	  pp->p.z = arcc[N-j];
	  pp->edge = 1;
	}
       else
	{ pp->p = xyzunit((XYZ){1,1,1});
	  pp->edge = 0;
	}
       p->nl = 0;
       link_between(pp,v[i+2]);
       link_between(pp,v[i+1]);
       link_between(pp,v[i]);
       v[i+1] = pp;
     }
    bcopy(&v[1],&v[0],(j+1)*sizeof(PT *));
    v[j+1] = 0;
  }
 if (ptx != NPT) abort();
 do dump(); while (adjust());
}

int main(void);
int main(void)
{
 find_p();
 return(0);
}