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

/*
	Create a casini asteroid for PovRay
	This is a two stage process, first apply some bulges
	Then apply crater impacts 
*/

#define N 300
#define M 200
TRIFACE *mesh = NULL;
int n = 0;

/* Impacts */
int impacts = 30;
double impact_rmin = 0.05, impact_rmax = 0.3;
double impact_dmin = 0.005, impact_dmax = 0.02;

/* Bulges */
int bulges = 30;
double bulge_rmin = 0.3,bulge_rmax = 0.7;
double bulge_dmin = -0.075,bulge_dmax = 0.075;

XYZ Eval(double,double,double,int,double);
XYZ XRotate(XYZ,double);
void GenRing(FILE *,double,double,double,double,double);

int main(int argc,char **argv)
{
	int i,j,k;
	int thesign;
	double a=0.95,b=1,tstart,tstop,t1,t2,theta1,theta2;
	double d,delta,len,radius,scale;
	XYZ p,q[4],norm;
	FILE *fptr;

	if (argc < 3) {
		fprintf(stderr,"Usage: %s bulges impacts\n",argv[0]);
		exit(-1);
	}
	bulges  = atoi(argv[1]);
	impacts = atoi(argv[2]);
	srandom(12345);

	/* Create the mesh */
   tstart = 0;
   if (a <= b)
      tstop = PI;
   else
      tstop = 0.5 * asin(b*b/(a*a));
	for (j=0;j<M;j++) {       /* Rotation about x */
      theta1 = TWOPI * j / (double)M;
      theta2 = TWOPI * ((j+1)%M) / (double)M;

      for (i=0;i<N;i++) {    /* Points on x-y plane */
         thesign = 1;
         if (a <= b) {
            t1 = tstart + (tstop - tstart) * i / (double)N;
            t2 = tstart + (tstop - tstart) * (i+1) / (double)N;
         } else {
            if (i < N/2) {
               t1 = tstart + 2 * (tstop - tstart) * i / (double)N;
               t2 = tstart + 2 * (tstop - tstart) * (i+1) / (double)N;
            } else {
               t1 = tstart + 2 * (tstop - tstart) * (i-N/2) / (double)N;
               t2 = tstart + 2 * (tstop - tstart) * (i+1-N/2) / (double)N;
               thesign = -1;
            }
         }

         /* Calculate the 4 points making up a facet */
         q[0] = Eval(t1,a,b,thesign,theta1);
         q[1] = Eval(t2,a,b,thesign,theta1);
         q[2] = Eval(t2,a,b,thesign,theta2); 
         q[3] = Eval(t1,a,b,thesign,theta2);
   
			if (i != N-1) {
				mesh = realloc(mesh,(n+1)*sizeof(TRIFACE));
				mesh[n].p[0] = q[0];
         	mesh[n].p[1] = q[2];
         	mesh[n].p[2] = q[1];
				n++;
			}
			if (i != 0) {
				mesh = realloc(mesh,(n+1)*sizeof(TRIFACE));
         	mesh[n].p[0] = q[0];
         	mesh[n].p[1] = q[3];
         	mesh[n].p[2] = q[2];
				n++;
			}
      }
   }
	fprintf(stderr,"Created %d triangles\n",n);

	/* Apply the bulges */
	for (k=0;k<bulges;k++) {
		fprintf(stderr,".");
		i = random() % n;
		j = random() % 3;
		p = mesh[i].p[j];
      scale  = bulge_dmin + drand48() * (bulge_dmax - bulge_dmin);
      radius = bulge_rmin + drand48() * (bulge_rmax - bulge_rmin);
      for (i=0;i<n;i++) {
         norm = CalcNormal(mesh[i].p[0],mesh[i].p[1],mesh[i].p[2]);
         for (j=0;j<3;j++) {
            if ((d = VectorLength(p,mesh[i].p[j])) >= radius)
               continue;
            len = scale * 0.5 * (1 + cos(PI * d / radius));
            mesh[i].p[j].x += norm.x * len;
            mesh[i].p[j].y += norm.y * len;
            mesh[i].p[j].z += norm.z * len;
         }
      }
	}
	fprintf(stderr,"\n");

   /* Apply the craters */
   for (k=0;k<impacts;k++) {
      fprintf(stderr,"*");
      i = random() % n;
      j = random() % 3;
      p = mesh[i].p[j];
		norm = CalcNormal(mesh[i].p[0],mesh[i].p[1],mesh[i].p[2]); 
      scale  = impact_dmin + drand48() * (impact_dmax - impact_dmin);
      radius = impact_rmin + drand48() * (impact_rmax - impact_rmin);
      for (i=0;i<n;i++) {
			/* norm = CalcNormal(mesh[i].p[0],mesh[i].p[1],mesh[i].p[2]); */
         for (j=0;j<3;j++) {
            d = VectorLength(p,mesh[i].p[j]);
				len = scale * cos(PI * d / radius) * exp(-(d*d)/(radius*radius));
            mesh[i].p[j].x -= norm.x * len;
            mesh[i].p[j].y -= norm.y * len;
            mesh[i].p[j].z -= norm.z * len;
         }
      }
   }
   fprintf(stderr,"\n");

	/* Write to PovRay */
	for (i=0;i<n;i++) {
      printf("triangle {\n");
		for (j=0;j<3;j++) {
     		printf("\t<%lf,%lf,%lf>", 
				mesh[i].p[j].x,mesh[i].p[j].y,mesh[i].p[j].z);
			if (j < 2)
				printf(",");
			printf("\n");
		}
      printf("\ttexture { asteroidmat }\n");
		printf("\trotate <0,clock*360,0>\n");
      printf("}\n");
	}

	/* Create rings 
	GenRing( 0.0,1.2,1.0,0.0,0.0);
	GenRing( 0.3,1.0,0.0,1.0,0.0);
	GenRing(-0.3,1.0,0.0,1.0,0.0);
	*/
	fptr = fopen("ring.inc","w");
	GenRing(fptr,0.0,1.2,0.0,1.0,0.0);
	fclose(fptr);
}

void GenRing(FILE *fptr,double x,double r,double red,double green,double blue)
{
	int i;
	XYZ q[2];

   for (i=0;i<360;i++) {
      q[0].x = x;
      q[0].y = r * cos(i*DTOR);
      q[0].z = r * sin(i*DTOR);
      q[1].x = x;
      q[1].y = r * cos((i+1)*DTOR);
      q[1].z = r * sin((i+1)*DTOR);
      fprintf(fptr,"cylinder {\n");
      fprintf(fptr,"\t<%g,%g,%g>, <%g,%g,%g>, 0.01\n",
            q[0].x,q[0].y,q[0].z,q[1].x,q[1].y,q[1].z);
      fprintf(fptr,"\ttexture { \n");
      fprintf(fptr,"\t\tpigment {\n");
      fprintf(fptr,"\t\t\trgb <%g,%g,%g>\n",red,green,blue);
      fprintf(fptr,"\t\t}\n");
      fprintf(fptr,"\t\tfinish { ringfinish }\n");
      fprintf(fptr,"\t}\n");
		fprintf(fptr,"\trotate <0,clock*360,0>\n");
      fprintf(fptr,"}\n");
      fprintf(fptr,"sphere {\n");
      fprintf(fptr,"\t<%g,%g,%g>, 0.01\n",
            q[0].x,q[0].y,q[0].z);
      fprintf(fptr,"\ttexture { \n");
		fprintf(fptr,"\t\tpigment {\n");
		fprintf(fptr,"\t\t\trgb <%g,%g,%g>\n",red,green,blue);
		fprintf(fptr,"\t\t}\n");
		fprintf(fptr,"\t\tfinish { ringfinish }\n");
		fprintf(fptr,"\t}\n");
		fprintf(fptr,"\trotate <0,clock*360,0>\n");
      fprintf(fptr,"}\n");
   }
}

/*
   Rotate a point about the x axis
*/
XYZ XRotate(XYZ p,double theta)
{
   XYZ q;

   q.x =  p.x;
   q.y =  p.y * cos(theta) + p.z * sin(theta);
   q.z = -p.y * sin(theta) + p.z * cos(theta);

   return(q);
}

/*
   Evaluate the Cassinian Oval at t given a and b
   t is assumed to be between the correct ranges
   For a <= b, 0 < t < TWOPI
   For a >  b, -to < t < to where to = 0.5 * asin(b^2/a^2)
*/
XYZ Eval(double t,double a,double b,int sign,double theta)
{
   XYZ p;
   double c1,c2;

   c1 = b*b*b*b - a*a*a*a*sin(2*t)*sin(2*t);
   if (c1 <= 0) /* Shouldn't happen with correct usage */
      c2 = a * a * cos(2*t);
   else
      c2 = a * a * cos(2*t) + sign * sqrt(c1);
   p.x = cos(t) * sqrt(c2);
   p.y = sin(t) * sqrt(c2);
   p.z = 0.0;

   /* Translate single loop case to origin */
   if (a > b)
      p.x -= a;

	p = XRotate(p,theta);
   return(p);
}