// If the user forgot to specify, use defaults
#ifndef (Smooth)
#local Smooth=true;
#end

#ifndef (Seed)
#local Seed=0;
#end

#local S=seed(Seed); // just in case the user didn't declare it properly

#ifndef (Splits)
#local Splits=1;
#end

#ifndef (Size)
#local Size=10;
#end

#ifndef (Rough)
#local Rough=.2;
#end

#ifndef (Omega)
#local Omega=.5;
#end

#local Tmax=20;
#local Emax=30;
#local Vmax=12;

#local I=0;
#while (I<Splits)
  #local Vmax=Vmax+Emax;
  #local Emax=Emax*2+Tmax*3;
  #local Tmax=Tmax*4;
#local I=I+1;
#end

//VC now contains the number of vertices in the final mesh
//EC now contains the number of edges in the final mesh
//TC now contains the number of triangles in the final mesh

#local VL=array[Vmax] // Location of each vertex
#if (Smooth!=false)
#local VN=array[Vmax] // Normal for each vertex
#end
#local E1=array[Emax] // lower-numbered endpoint for edge
#local E2=array[Emax] // higher-numbered endpoint for edge
#local T1=array[Tmax] // Edge 1 for triangle
#local T2=array[Tmax] // Edge 2 for triangle
#local T3=array[Tmax] // Edge 3 for triangle

#local VC=12;
#local EC=30;
#local TC=20;

//now determine the original vertices

#local qA=sqrt((3+sqrt(5))/(5+sqrt(5)));
#local qB=sqrt(2/(5+sqrt(5)));
#local VL[ 0]=< qA, qB,0>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 1]=<-qA, qB,0>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 2]=< qA,-qB,0>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 3]=<-qA,-qB,0>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 4]=<0, qA, qB>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 5]=<0,-qA, qB>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 6]=<0, qA,-qB>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 7]=<0,-qA,-qB>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 8]=< qB,0, qA>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[ 9]=< qB,0,-qA>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[10]=<-qB,0, qA>*Size*(1+Rough*(rand(S)-rand(S)));
#local VL[11]=<-qB,0,-qA>*Size*(1+Rough*(rand(S)-rand(S)));

#local I=12;
#while (I<Vmax)
  #local VL[I]=<0,0,0>;
#local I=I+1;
#end

//and the original edges
#local E1[ 0]= 0; #local E2[ 0]= 2;
#local E1[ 1]= 1; #local E2[ 1]= 3;
#local E1[ 2]= 4; #local E2[ 2]= 6;
#local E1[ 3]= 5; #local E2[ 3]= 7;
#local E1[ 4]= 8; #local E2[ 4]=10;
#local E1[ 5]= 9; #local E2[ 5]=11;
#local E1[ 6]= 2; #local E2[ 6]= 5;
#local E1[ 7]= 2; #local E2[ 7]= 7;
#local E1[ 8]= 2; #local E2[ 8]= 8;
#local E1[ 9]= 2; #local E2[ 9]= 9;
#local E1[10]= 1; #local E2[10]= 4;
#local E1[11]= 1; #local E2[11]= 6;
#local E1[12]= 1; #local E2[12]=10;
#local E1[13]= 1; #local E2[13]=11;
#local E1[14]= 0; #local E2[14]= 4;
#local E1[15]= 0; #local E2[15]= 8;
#local E1[16]= 0; #local E2[16]= 9;
#local E1[17]= 0; #local E2[17]= 6;
#local E1[18]= 3; #local E2[18]= 5;
#local E1[19]= 3; #local E2[19]= 7;
#local E1[20]= 3; #local E2[20]=10;
#local E1[21]= 3; #local E2[21]=11;
#local E1[22]= 4; #local E2[22]= 8;
#local E1[23]= 4; #local E2[23]=10;
#local E1[24]= 5; #local E2[24]= 8;
#local E1[25]= 5; #local E2[25]=10;
#local E1[26]= 6; #local E2[26]= 9;
#local E1[27]= 6; #local E2[27]=11;
#local E1[28]= 7; #local E2[28]= 9;
#local E1[29]= 7; #local E2[29]=11;

#local I=30;
#while (I<Emax)
  #local E1[I]=0;
  #local E2[I]=0;
#local I=I+1;
#end

#local T1[ 0]= 0; #local T2[ 0]= 8; #local T3[ 0]=15;
#local T1[ 1]= 0; #local T2[ 1]=16; #local T3[ 1]= 9;
#local T1[ 2]= 9; #local T2[ 2]=28; #local T3[ 2]= 7;
#local T1[ 3]= 3; #local T2[ 3]= 6; #local T3[ 3]= 7;
#local T1[ 4]= 6; #local T2[ 4]=24; #local T3[ 4]= 8;
#local T1[ 5]= 3; #local T2[ 5]=19; #local T3[ 5]=18;
#local T1[ 6]= 4; #local T2[ 6]=24; #local T3[ 6]=25;
#local T1[ 7]=14; #local T2[ 7]=15; #local T3[ 7]=22;
#local T1[ 8]=16; #local T2[ 8]=17; #local T3[ 8]=26;
#local T1[ 9]= 5; #local T2[ 9]=29; #local T3[ 9]=28;
#local T1[10]=14; #local T2[10]= 2; #local T3[10]=17;
#local T1[11]= 5; #local T2[11]=26; #local T3[11]=27;
#local T1[12]=19; #local T2[12]=29; #local T3[12]=21;
#local T1[13]=18; #local T2[13]=20; #local T3[13]=25;
#local T1[14]= 4; #local T2[14]=23; #local T3[14]=22;

#local T1[15]=10; #local T2[15]=23; #local T3[15]=12;
#local T1[16]=11; #local T2[16]= 2; #local T3[16]=10;
#local T1[17]=13; #local T2[17]=27; #local T3[17]=11;
#local T1[18]= 1; #local T2[18]=21; #local T3[18]=13;
#local T1[19]=12; #local T2[19]=20; #local T3[19]= 1;

#local I=30;
#while (I<Tmax)
  #local T1[I]=0;
  #local T2[I]=0;
  #local T3[I]=0;
#local I=I+1;
#end

//The icosahedron is now as complete as we need it to be
//for the fractal divisions, which commence now.

//init the tR value for the level about to happen
#local tR=Rough*Size;

#local Depth=0;
#while (Depth<Splits)

#local tR=tR*Omega;

//logic to determine new center points for each edge
//and then split the edges

#local I=0;
#while (I<EC)
#local P1=VL[E1[I]];
#local P2=VL[E2[I]];

  #local VL[I+VC]=vnormalize(P1+P2)*
    ((vlength(P1)+vlength(P2))/2+(rand(S)-rand(S))*tR);

  #local E2[I+EC]=I+VC;
  #local E1[I+EC]=E2[I];
  #local E2[I]=I+VC;

#local I=I+1;
#end

//now go through each triangle, making new edges and daughter
//triangles

#local I=0;
#while (I<TC)

#local Edge1=T1[I];
#local Edge2=T2[I];
#local Edge3=T3[I];
#local Edge4=T1[I]+EC;
#local Edge5=T2[I]+EC;
#local Edge6=T3[I]+EC;
#local Edge7=EC+EC+I;
#local Edge8=EC+EC+TC+I;
#local Edge9=EC+EC+TC+TC+I;


//logic to make sure that edges can be assigned in the
//proper order
#if (E1[Edge1]!=E1[Edge6])
#local TEMP=Edge1; #local Edge1=Edge4; #local Edge4=TEMP;
#end

#if (E1[Edge1]!=E1[Edge6])
#local TEMP=Edge3; #local Edge3=Edge6; #local Edge6=TEMP;
#end

#if (E1[Edge1]!=E1[Edge6])
#local TEMP=Edge1; #local Edge1=Edge4; #local Edge4=TEMP;
#end

#if (E1[Edge2]!=E1[Edge4])
#local TEMP=Edge2; #local Edge2=Edge5; #local Edge5=TEMP;
#end

//create the interior edges
#local E1[Edge7]=E2[Edge2];
#local E2[Edge7]=E2[Edge3];
#local E1[Edge8]=E2[Edge3];
#local E2[Edge8]=E2[Edge1];
#local E1[Edge9]=E2[Edge1];
#local E2[Edge9]=E2[Edge2];

//create the four daughter triangles
#local T1[I]=Edge1;
#local T2[I]=Edge8;
#local T3[I]=Edge6;
#local T1[I+TC]=Edge2;
#local T2[I+TC]=Edge9;
#local T3[I+TC]=Edge4;
#local T1[I+TC+TC]=Edge3;
#local T2[I+TC+TC]=Edge7;
#local T3[I+TC+TC]=Edge5;
#local T1[I+TC+TC+TC]=Edge7;
#local T2[I+TC+TC+TC]=Edge8;
#local T3[I+TC+TC+TC]=Edge9;

#local I=I+1;
#end

//calculate the new number of vertices, edges, and triangles
#local VC=VC+EC;
#local EC=EC+EC+TC+TC+TC;
#local TC=TC+TC+TC+TC;

#local Depth=Depth+1;
#end

#if (Smooth!=false)

//Now the triangle normals get calculated

// clear out the normals for each vertex
#local I= 0; #while (I<Vmax) #local VN[I]=<0,0,0>; #local I=I+1; #end

//Now we calculate the surface normals for each triangle

#local I=0;
#while (I<TC)

#local V1=E1[T1[I]];
#local V2=E2[T1[I]];
#local V3=E1[T2[I]];
#local V4=E2[T2[I]];

//a bit of shuffling logic to ensure that V1 and V3
//are the same vertex

#if (V1!=V3)
  #local TEMP=V3; #local V3=V4; #local V4=TEMP;
#end

#if (V1!=V3)
  #local TEMP=V1; #local V1=V2; #local V2=TEMP;
#end

#if (V1!=V3)
  #local TEMP=V3; #local V3=V4; #local V4=TEMP;
#end

#local Norm=vcross(VL[V2]-VL[V1],VL[V4]-VL[V1]);
#local VN[V1]=VN[V1]+Norm;
#local VN[V2]=VN[V2]+Norm;
#local VN[V4]=VN[V4]+Norm;

#local I=I+1;
#end

#end // end of normal calculation
// Now that we have the normals for smoothing, we can
//actually declare the triangle mesh

mesh {
#local I=0;
#while (I<TC)

#local V1=E1[T1[I]];
#local V2=E2[T1[I]];
#local V3=E1[T2[I]];
#local V4=E2[T2[I]];

//a bit of shuffling logic to ensure that V1 and V3
//are the same vertex

#if (V1!=V3)
  #local TEMP=V3; #local V3=V4; #local V4=TEMP;
#end

#if (V1!=V3)
  #local TEMP=V1; #local V1=V2; #local V2=TEMP;
#end

#if (V1!=V3)
  #local TEMP=V3; #local V3=V4; #local V4=TEMP;
#end

#if (Smooth=false)
  triangle { VL[V1],VL[V2],VL[V4] }
#else
  smooth_triangle { VL[V1],VN[V1],VL[V2],VN[V2],VL[V4],VN[V4] }
#end
#local I=I+1;
#end
}
