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

#define TWOPI 6.283185307179586476925287
#define RADIUS 0.02

FILE *f1,*f2;
int NSEGMENTS = 400;
int SPHEREN = 16;
double r2 = 0.15;
int nfacet = 0;

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

int main(int argc,char **argv)
{
   int i;
   double x,y,z;
   double xlast,ylast,zlast;
   double xmin=1e32,xmax=-1e32,ymin=1e32,ymax=-1e32,zmin=1e32,zmax=-1e32;
   double dd,ddmax=0,ddmin=1e32;
   double mu;
   int nlongitude,nmeridian;

   /* Command line arguments */
   if (argc < 3) {
      fprintf(stderr,"usage: %s nmeridian nlongitude\n",argv[0]);
      exit(-1);
   }
   nmeridian  = atoi(argv[1]);
   nlongitude = atoi(argv[2]);

   /* Open two files, one for the geom file, one for the stl file */
   f1 = fopen("knot3.geom","w");
   f2 = fopen("knot3.stl","w");
   fprintf(f2,"solid\n");

   for (i=0;i<NSEGMENTS;i++) {

      mu = i * TWOPI * nmeridian / (double)NSEGMENTS;
      x = cos(mu) * (1 + cos(nlongitude*mu/(double)nmeridian) / 2.0); 
      y = sin(mu) * (1 + cos(nlongitude*mu/(double)nmeridian) / 2.0);
      z = sin(nlongitude*mu/(double)nmeridian) / 2.0;

      MakeSphere(x,y,z,r2);

      if (i > 0) {
         dd = (x-xlast)*(x-xlast) + (y-ylast)*(y-ylast) + (z-zlast)*(z-zlast);
         ddmin = (dd < ddmin ? dd : ddmin);
         ddmax = (dd > ddmax ? dd : ddmax);
      }
      xlast = x;
      ylast = y;
      zlast = z;

      xmin = (x < xmin ? x : xmin);
      xmax = (x > xmax ? x : xmax);
      ymin = (y < ymin ? y : ymin);
      ymax = (y > ymax ? y : ymax);
      zmin = (z < zmin ? z : zmin);
      zmax = (z > zmax ? z : zmax);
   }
   fprintf(f2,"endsolid\n");

   fprintf(stderr,"Sphere separation range = %g -> %g\n",ddmin,ddmax);
   fprintf(stderr,"Sphere radius = %g\n",r2);
   fprintf(stderr,"Number of spheres = %d\n",NSEGMENTS);
   fprintf(stderr,"Total facets = %d\n",nfacet); 
   fprintf(stderr,"X range = %g -> %g\n",xmin,xmax);
   fprintf(stderr,"Y range = %g -> %g\n",ymin,ymax);
   fprintf(stderr,"Z range = %g -> %g\n",zmin,zmax);

   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,"\tvertex %g %g %g\n",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f2,"\tvertex %g %g %g\n",cx+r*e[2].x,cy+r*e[2].y,cz+r*e[2].z);
         fprintf(f2,"\tvertex %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,"\tvertex %g %g %g\n",cx+r*e[0].x,cy+r*e[0].y,cz+r*e[0].z);
         fprintf(f2,"\tvertex %g %g %g\n",cx+r*e[3].x,cy+r*e[3].y,cz+r*e[3].z);
         fprintf(f2,"\tvertex %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");

         nfacet += 2;
      }
   }
}