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

#define MAXNUM_POINTS		20000		 /* maximum number of points */
#define MAXNUM_TRIANGLES	20000	      /* maximum number of triangles */

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

long rand_i_seed;				       /* random number seed */


struct Point{			      
 double x;
 double y;
 double z;					      
} Point;

struct Triangle{		
  long p1;	
  long p2;
  long p3;
};

char *progname;


/* Display an error message on stderr, prefixed with the program name, and
 * suffixed with strerror()'s value; then sleep for 2 seconds before returning.
 */
/*
void e_msg(char *fmt,...) {
 va_list ap;
 va_start(ap,fmt);
 fprintf(stderr,"%s: ",progname);
 vfprintf(stderr,fmt,ap);
 fprintf(stderr," (error: %s)\n",strerror(errno));
 va_end(ap);
}
*/

/* Similar to e_msg() above, but no suffix
 */
/*
void i_msg(char *fmt,...) {
 va_list ap;
 va_start(ap,fmt);
 fprintf(stderr,"%s: ",progname);
 vfprintf(stderr,fmt,ap);
 fprintf(stderr,"\n");
 va_end(ap);
}
*/


/* Uniform Random Number Generator
 *
 * From Numerical Recipies
 */
double rand_uniform(void) {
 int j;
 long k;
 static long iy=0;
 static long iv[NTAB];
 double temp;
 if (rand_i_seed<=0||!iy) {
  if ((-rand_i_seed)<1) rand_i_seed=1; else rand_i_seed=-rand_i_seed;
  for(j=NTAB+7;j>=0;j--) {
   k=(rand_i_seed)/IQ;
   rand_i_seed=IA*(rand_i_seed-k*IQ)-IR*k;
   if (rand_i_seed<0) rand_i_seed+=IM;
   if (j<NTAB) iv[j]=rand_i_seed;
  }
  iy=iv[0];
 }
 k=(rand_i_seed)/IQ;
 rand_i_seed=IA*(rand_i_seed-k*IQ)-IR*k;
 if (rand_i_seed<0) rand_i_seed+=IM;
 j=iy/NDIV;
 iy=iv[j];
 iv[j]=rand_i_seed;
 if ((temp=AM*iy)>RNMX) return(RNMX); else return(temp);
}


/* Gaussian Random Number Generator
 *
 * From Numerical Recipies
 */
double rand_gaussian(void) {
 static int iset=0;
 static double gset;
 double fac,rsq,v1,v2;
 if (iset==0) {
  do {
   v1=2.0*rand_uniform()-1.0;
   v2=2.0*rand_uniform()-1.0;
   rsq=v1*v1+v2*v2;
  } while (rsq>=1.0||rsq==0.0);
  fac=sqrt(-2.0*log(rsq)/rsq);
  gset=v1*fac;
  iset=1;
  return(v2*fac);
 } else {
  iset=0;
  return(gset);
 }
}



void vec_split(long p1,long p2,long t1,long t2,long *npoints,long *ntriangles,
               struct Point p[],struct Triangle t[],float ptrb,float ptrb_pwr,
	       float sx,float sy,float sz) {
 float dx,dy,dz,l,r1,r2,r,theta;
 /* calcute unperturbed midpoint - and perturber scale */
 /*
 r1=sqrt(p[p1].x*p[p1].x+p[p1].y*p[p1].y+p[p1].z*p[p1].z);
 r2=sqrt(p[p2].x*p[p2].x+p[p2].y*p[p2].y+p[p2].z*p[p2].z);
 dx=p[p1].x-p[p2].x; dy=p[p1].y-p[p2].y; dz=p[p1].z-p[p2].z;
 l=(dx*dx+dy*dy+dz*dz)*ptrb;
 dx=0.5*(p[p1].x+p[p2].x); dy=0.5*(p[p1].y+p[p2].y); dz=0.5*(p[p1].z+p[p2].z);
 r=sqrt(dx*dx+dy*dy+dz*dz);
 r=(r1+r2)/(2*r);
 dx*=r; dy*=r; dz*=r;
 */
 /*
 dx=p[p1].x-p[p2].x; dy=p[p1].y-p[p2].y; dz=p[p1].z-p[p2].z;
 l=sqrt(dx*dx+dy*dy+dz*dz)*ptrb;
 r1=sqrt(p[p1].x*p[p1].x+p[p1].y*p[p1].y+p[p1].z*p[p1].z);
 r2=sqrt(p[p2].x*p[p2].x+p[p2].y*p[p2].y+p[p2].z*p[p2].z);
 theta=((p[p1].x*p[p2].x+p[p1].y*p[p2].y+p[p1].z*p[p2].z)/(r1*r2));
 theta=acos(theta);
 theta=1/(2*cos(theta*0.5));
 dx=(p[p1].x+p[p2].x)*theta; dy=(p[p1].y+p[p2].y)*theta; dz=(p[p1].z+p[p2].z)*theta;
 */
 dx=p[p1].x-p[p2].x; dy=p[p1].y-p[p2].y; dz=p[p1].z-p[p2].z;
 l=sqrt(dx*dx+dy*dy+dz*dz);
 r1=sqrt(p[p1].x*p[p1].x/(sx*sx)+p[p1].y*p[p1].y/(sy*sy)+p[p1].z*p[p1].z/(sz*sz));
 r2=sqrt(p[p2].x*p[p2].x/(sx*sx)+p[p2].y*p[p2].y/(sy*sy)+p[p2].z*p[p2].z/(sz*sz));
 theta=((p[p1].x*p[p2].x/(sx*sx)+p[p1].y*p[p2].y/(sy*sy)+p[p1].z*p[p2].z/(sz*sz))/(r1*r2));
 theta=acos(theta);
 theta=1/(2*cos(theta*0.5));
 l=pow(ptrb*l,ptrb_pwr);
 if(l>1) l=1;
 dx=(p[p1].x+p[p2].x)*theta; dy=(p[p1].y+p[p2].y)*theta; dz=(p[p1].z+p[p2].z)*theta;
 r=1+rand_gaussian()*l;
 while(r<=0){
    r=1+rand_gaussian()*l;
 }
 dx=dx*r; dy=dy*r;  dz=dz*r;
 *npoints=*npoints+1;
 p[*npoints].x=dx; p[*npoints].y=dy; p[*npoints].z=dz;
 *ntriangles=*ntriangles+1;
 t[*ntriangles].p1=t[t1].p1;
 t[*ntriangles].p2=t[t1].p2;
 t[*ntriangles].p3=t[t1].p3;
 if (t[t1].p1==p1) t[t1].p1=*npoints;
 else if (t[t1].p2==p1) t[t1].p2=*npoints;
 else if (t[t1].p3==p1) t[t1].p3=*npoints;
 /*else i_msg("vec_split: error [t1,p1]");*/
 if (t[*ntriangles].p1==p2) t[*ntriangles].p1=*npoints;
 else if (t[*ntriangles].p2==p2) t[*ntriangles].p2=*npoints;
 else if (t[*ntriangles].p3==p2) t[*ntriangles].p3=*npoints;
 /*else i_msg("vec_split: error [ntriangles,p2]");*/
 *ntriangles=*ntriangles+1;
 t[*ntriangles].p1=t[t2].p1;
 t[*ntriangles].p2=t[t2].p2;
 t[*ntriangles].p3=t[t2].p3;
 if (t[t2].p1==p1) t[t2].p1=*npoints;
 else if (t[t2].p2==p1) t[t2].p2=*npoints;
 else if (t[t2].p3==p1) t[t2].p3=*npoints;
 /*else i_msg("vec_split: error [t2,p1]");*/
 if (t[*ntriangles].p1==p2) t[*ntriangles].p1=*npoints;
 else if (t[*ntriangles].p2==p2) t[*ntriangles].p2=*npoints;
 else if (t[*ntriangles].p3==p2) t[*ntriangles].p3=*npoints;
 /*else i_msg("vec_split: error [t2,p2]");*/
}


void vec_select(long npoints,long ntriangles,struct Point p[],
                struct Triangle t[],long *p1,long *p2,long *t1,long *t2) {
 double max_len,dx,dy,dz,len;
 long i,new;
 max_len=0;
 for (i=0; i<=ntriangles; i++) {
  new=0;						 /* measure vectors: */
  dx=p[t[i].p1].x-p[t[i].p2].x;
  dy=p[t[i].p1].y-p[t[i].p2].y;
  dz=p[t[i].p1].z-p[t[i].p2].z;
  len=dx*dx+dy*dy+dz*dz;
  if (len>max_len) {
   *p1=t[i].p1; *p2=t[i].p2; *t1=i; new=1; max_len=len;
  }
  dx=p[t[i].p2].x-p[t[i].p3].x;
  dy=p[t[i].p2].y-p[t[i].p3].y;
  dz=p[t[i].p2].z-p[t[i].p3].z;
  len=dx*dx+dy*dy+dz*dz;
  if (len>max_len) {
   *p1=t[i].p2; *p2=t[i].p3; *t1=i; new=1; max_len=len;
  }
  dx=p[t[i].p3].x-p[t[i].p1].x;
  dy=p[t[i].p3].y-p[t[i].p1].y;
  dz=p[t[i].p3].z-p[t[i].p1].z;
  len=dx*dx+dy*dy+dz*dz;
  if (len>max_len) {
   *p1=t[i].p3; *p2=t[i].p1; *t1=i; new=1; max_len=len;
  }
  if (new==0) {			 /* if not shorter, test triangle for match: */
   if((t[i].p1==*p1)||(t[i].p2==*p1)||(t[i].p3==*p1)) {
    if((t[i].p1==*p2)||(t[i].p2==*p2)||(t[i].p3==*p2))
     *t2=i;
   }
  }
 }
}


void calc_normals(long npoints,long ntriangles,struct Point p[],
                  struct Point v[],struct Triangle t[]) {
 long i,*n;
 double x1,y1,z1,x2,y2,z2,vx,vy,vz,l;
 n=(long *)malloc(sizeof(long)*npoints);
 if (n==NULL) {/*e_msg("calc_normals: allocation failed");*/ return;}
 for(i=0; i<=npoints; i++) {v[i].x=0; v[i].y=0; v[i].z=0; n[i]=0;}
 for(i=0; i<=ntriangles; i++) {
  x1=p[t[i].p1].x-p[t[i].p2].x;
  y1=p[t[i].p1].y-p[t[i].p2].y;
  z1=p[t[i].p1].z-p[t[i].p2].z;
  x2=p[t[i].p2].x-p[t[i].p3].x;
  y2=p[t[i].p2].y-p[t[i].p3].y;
  z2=p[t[i].p2].z-p[t[i].p3].z;
  vx=y1*z2-y2*z1; vy=x2*z1-x1*z2; vz=x1*y2-x2*y1;
  l=sqrt(vx*vx+vy*vy+vz*vz);
  vx/=l; vy/=l; vz/=l;
  v[t[i].p1].x=v[t[i].p1].x*n[t[i].p1]+vx;
  v[t[i].p1].y=v[t[i].p1].y*n[t[i].p1]+vy;
  v[t[i].p1].z=v[t[i].p1].z*n[t[i].p1]+vz;
  n[t[i].p1]=n[t[i].p1]+1;
  v[t[i].p2].x=v[t[i].p2].x*n[t[i].p2]+vx;
  v[t[i].p2].y=v[t[i].p2].y*n[t[i].p2]+vy;
  v[t[i].p2].z=v[t[i].p2].z*n[t[i].p2]+vz;
  n[t[i].p2]=n[t[i].p2]+1;
  v[t[i].p3].x=v[t[i].p3].x*n[t[i].p3]+vx;
  v[t[i].p3].y=v[t[i].p3].y*n[t[i].p3]+vy;
  v[t[i].p3].z=v[t[i].p3].z*n[t[i].p3]+vz;
  n[t[i].p3]=n[t[i].p3]+1;
 }
 free(n);
}

/* Set up initial array of points and triangles and seed the random number
 * generator, then split vectors "subdivide" times, adding noise, and finally
 * calculate the surface normals and output the mesh in POV object format.
 */
void main(int argc,char *argv[]) {
 struct Point points[MAXNUM_POINTS];
 struct Point pvecs[MAXNUM_POINTS];
 struct Triangle triangles[MAXNUM_TRIANGLES];
 int subdivisions;
 float xsize,ysize,zsize,ptrb,ptrb_pwr;
 long npoints,ntriangles;
 long p1,p2,t1,t2,pt;
 int i;
 progname=argv[0];
 subdivisions = 1994;
 rand_i_seed = getpid();
 xsize = 1;
 ysize = 1;
 zsize = 1;
 ptrb = 0.05;
 ptrb_pwr = 0;

 if (argc==1){
   fprintf(stderr,"Usage: %s subdivisions seed x y z prtb pwr\n",argv[0]);
   fprintf(stderr,"Using defaults\n");
 } 
 if(argc>1){
   subdivisions=atoi(argv[1]);
 }
 if(argc>2){
   rand_i_seed=atol(argv[2]);
 }
 if(argc>5){
   xsize=atof(argv[3]);
   ysize=atof(argv[4]);
   zsize=atof(argv[5]);
 }
 if(argc>6){
   ptrb=atof(argv[6]);
 }
 if(argc>7){
   ptrb_pwr = atof(argv[7]);
 }
 points[0].x=xsize;  points[0].y=0; points[0].z=0;
 points[1].x=-xsize; points[1].y=0; points[1].z=0;
 points[2].x=0; points[2].y=ysize;  points[2].z=0;
 points[3].x=0; points[3].y=-ysize; points[3].z=0;
 points[4].x=0; points[4].y=0; points[4].z=zsize;
 points[5].x=0; points[5].y=0; points[5].z=-zsize;
 npoints=5;
 triangles[0].p1=0; triangles[0].p2=4; triangles[0].p3=2;
 triangles[1].p1=0; triangles[1].p2=3; triangles[1].p3=4;
 triangles[2].p1=0; triangles[2].p2=5; triangles[2].p3=3;
 triangles[3].p1=0; triangles[3].p2=2; triangles[3].p3=5;
 triangles[4].p1=1; triangles[4].p2=2; triangles[4].p3=4;
 triangles[5].p1=1; triangles[5].p2=4; triangles[5].p3=3;
 triangles[6].p1=1; triangles[6].p2=3; triangles[6].p3=5;
 triangles[7].p1=1; triangles[7].p2=5; triangles[7].p3=2;
 ntriangles=7;
 for(i=0; i<subdivisions; i++) {
  vec_select(npoints,ntriangles,points,triangles,&p1,&p2,&t1,&t2);
  vec_split(p1,p2,t1,t2,&npoints,&ntriangles,points,triangles,ptrb,ptrb_pwr,xsize,ysize,zsize);
  if ((i%100)==0) fprintf(stderr,"i=%d [points: %ld]\n",i,npoints);
 }
 calc_normals(npoints,ntriangles,points,pvecs,triangles);
 printf("#declare ASTEROID = mesh {\n");
 for(i=0; i<=ntriangles; i++) {
  printf(" smooth_triangle {\n");
  pt=triangles[i].p1;
  printf("  <%f,%f,%f>,",points[pt].x,points[pt].y,points[pt].z);
  printf("<%f,%f,%f>,\n",pvecs[pt].x,pvecs[pt].y,pvecs[pt].z);
  pt=triangles[i].p2;
  printf("  <%f,%f,%f>,",points[pt].x,points[pt].y,points[pt].z);
  printf("<%f,%f,%f>,\n",pvecs[pt].x,pvecs[pt].y,pvecs[pt].z);
  pt=triangles[i].p3;
  printf("  <%f,%f,%f>,",points[pt].x,points[pt].y,points[pt].z);
  printf("<%f,%f,%f>\n",pvecs[pt].x,pvecs[pt].y,pvecs[pt].z);
  printf(" }\n");
 }
 printf(" texture {pigment {rgb <1,1,1>}}\n}\n");
}