#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "paulslib.h"

/*
   Create a borromean ring in STL format
*/

void MakeSphere(double,double,double,double);

FILE *f1,*f2;
int NSEGMENTS = 100;
int SPHEREN = 16;

int main(int argc,char **argv)
{
   int i;
   double x,y,z;
   double mu,r1,r2;

   if (argc > 1)
      NSEGMENTS = atoi(argv[1]);
   if (argc > 2)
      SPHEREN = atoi(argv[2]);

   f1 = fopen("borromean.geom","w");
   f2 = fopen("borromean.stl","w");

   r1 = sqrt(3.0) / 3.0;
   r2 = 0.2;

   for (i=0;i<NSEGMENTS;i++) {
      mu = i * TWOPI / (double)NSEGMENTS;

      x = cos(mu);
      y = sin(mu) + r1;
      z = -cos(3*mu) / 3;
      MakeSphere(x,y,z,r2);

      x = cos(mu) + 0.5; 
      y = sin(mu) - r1/2;
      z = -cos(3*mu) / 3;
      MakeSphere(x,y,z,r2);

      x = cos(mu) - 0.5;
      y = sin(mu) - r1/2;
      z = -cos(3*mu) / 3;
      MakeSphere(x,y,z,r2);
   }

   fclose(f1);
   fclose(f2);
}

void MakeSphere(double cx,double cy,double cz,double r)
{
   int i,j;
   double t1,t2,t3,t4;
   XYZ e[4],normal;

   for (j=0;j<SPHEREN/2;j++) {
      t1 = -PID2 + j * PI / (SPHEREN/2);
      t2 = -PID2 + (j + 1) * PI / (SPHEREN/2);

      for (i=0;i<SPHEREN;i++) {
         t3 = i * TWOPI / SPHEREN;
         t4 = ((i+1) % SPHEREN) * TWOPI / SPHEREN;

         /* The sphere at the origin */
         e[0].x = cos(t1) * cos(t3);
         e[0].y = sin(t1);
         e[0].z = cos(t1) * sin(t3);
         e[1].x = cos(t2) * cos(t3);
         e[1].y = sin(t2);
         e[1].z = cos(t2) * sin(t3);
         e[2].x = cos(t2) * cos(t4);
         e[2].y = sin(t2);
         e[2].z = cos(t2) * sin(t4);
         e[3].x = cos(t1) * cos(t4);
         e[3].y = sin(t1);
         e[3].z = cos(t1) * sin(t4);

         /* geom format for face 1 */
         fprintf(f1,"f3n ");
         fprintf(f1,"%g %g %g ",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f1,"%g %g %g ",cx+r*e[2].x,cy+r*e[2].y,cz+r*e[2].z); 
         fprintf(f1,"%g %g %g ",cx+r*e[1].x,cy+r*e[1].y,cz+r*e[1].z); 
         fprintf(f1,"%g %g %g ",e[0].x,e[0].y,e[0].z);
         fprintf(f1,"%g %g %g ",e[2].x,e[2].y,e[2].z);
         fprintf(f1,"%g %g %g ",e[1].x,e[1].y,e[1].z);
         fprintf(f1,"1 0 0\n");

         /* STL format for face 1 */
         normal = CalcNormal(e[0],e[2],e[1]);
         fprintf(f2,"facet normal %g %g %g\n",normal.x,normal.y,normal.z);
         fprintf(f2,"outer loop\n");
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[2].x,cy+r*e[2].y,cz+r*e[2].z);
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[1].x,cy+r*e[1].y,cz+r*e[1].z);
         fprintf(f2,"endloop\n");
         fprintf(f2,"endfacet\n");

         /* geom format for face 2 */
         fprintf(f1,"f3n ");
         fprintf(f1,"%g %g %g ",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f1,"%g %g %g ",cx+r*e[3].x,cy+r*e[3].y,cz+r*e[3].z);
         fprintf(f1,"%g %g %g ",cx+r*e[2].x,cy+r*e[2].y,cz+r*e[2].z);
         fprintf(f1,"%g %g %g ",e[0].x,e[0].y,e[0].z);
         fprintf(f1,"%g %g %g ",e[3].x,e[3].y,e[3].z);
         fprintf(f1,"%g %g %g ",e[2].x,e[2].y,e[2].z);
         fprintf(f1,"1 0 0\n");

         /* STL format for face 2 */
         normal = CalcNormal(e[0],e[3],e[2]);
         fprintf(f2,"facet normal %g %g %g\n",normal.x,normal.y,normal.z);
         fprintf(f2,"outer loop\n");
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[3].x,cy+r*e[3].y,cz+r*e[3].z);
         fprintf(f2,"   vertex %g %g %g\n",cx+r*e[2].x,cy+r*e[2].y,cz+r*e[2].z);
         fprintf(f2,"endloop\n");
         fprintf(f2,"endfacet\n");

      }
   }
}