#macro ScalarRange(Seed)
  (rand(Seed)-rand(Seed))
#end

#macro VectorRange(Seed)
  (<rand(Seed),rand(Seed),rand(Seed)>-<rand(Seed),rand(Seed),rand(Seed)>)
#end

#macro RoundedMid(Array,K,Ia,Ja,Ib,Jb,Ir,Jr)
  #local Array[K][Ir][Jr]=vnormalize(Array[K][Ia][Ja]+Array[K][Ib][Jb])*((vlength(Array[K][Ia][Ja])+vlength(Array[K][Ib][Jb]))/2+ScalarRange(S)*tR);
#end

#macro RoughMid(Array,K,Ia,Ja,Ib,Jb,Ir,Jr)
  #local Array[K][Ir][Jr]=(Array[K][Ia][Ja]+Array[K][Ib][Jb]))/2+VectorRange(S)*tR);
#end

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

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

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

#ifndef (Size)
#local Size=10;
#end
/*
#ifndef (Start)
#local Start=1;
#end
*/
#ifndef (Rough)
#local Rough=.2;
#end

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

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

#local SC=1; // a commonly used value

#local SM=1; // just as commonly used

#local I=0; #while (I<Splits) #local SM=SM+SM; #local I=I+1; #end

#local P=array[5][SM+1][SM+SM+1]

#if (Smooth) #local N=array[5][SM+1][SM+SM+1] #end

// initialize the fractals.  There are five fractal arrays,
// which fit together to cover the surface of the sphere

#local qA=(sqrt(5)+1)/4;
#local qB=1/2/sqrt(qA*qA+1/4);
#local qA=qA/sqrt(qA*qA+1/4);

#local P[0][0][0]=< qA, qB,0>*Size*(ScalarRange(S)*Rough+1); //0
#local P[1][0][0]=P[0][0][0];
#local P[2][0][0]=P[0][0][0];
#local P[3][0][0]=P[0][0][0];
#local P[4][0][0]=P[0][0][0];

#local P[2][1][1]=<-qA, qB,0>*Size*(ScalarRange(S)*Rough+1); //1
#local P[3][0][2]=P[2][1][1];

#local P[4][1][0]=< qA,-qB,0>*Size*(ScalarRange(S)*Rough+1); //2
#local P[0][0][1]=P[4][1][0];

#local P[0][1][2]=<-qA,-qB,0>*Size*(ScalarRange(S)*Rough+1); //3
#local P[1][1][2]=P[0][1][2];
#local P[2][1][2]=P[0][1][2];
#local P[3][1][2]=P[0][1][2];
#local P[4][1][2]=P[0][1][2];

#local P[1][1][0]=<0, qA, qB>*Size*(ScalarRange(S)*Rough+1); //4
#local P[2][0][1]=P[1][1][0];

#local P[0][1][1]=<0,-qA, qB>*Size*(ScalarRange(S)*Rough+1); //5
#local P[1][0][2]=P[0][1][1];

#local P[2][1][0]=<0, qA,-qB>*Size*(ScalarRange(S)*Rough+1); //6
#local P[3][0][1]=P[2][1][0];

#local P[4][1][1]=<0,-qA,-qB>*Size*(ScalarRange(S)*Rough+1); //7
#local P[0][0][2]=P[4][1][1];

#local P[0][1][0]=< qB,0, qA>*Size*(ScalarRange(S)*Rough+1); //8
#local P[1][0][1]=P[0][1][0];

#local P[3][1][0]=< qB,0,-qA>*Size*(ScalarRange(S)*Rough+1); //9
#local P[4][0][1]=P[3][1][0];

#local P[1][1][1]=<-qB,0, qA>*Size*(ScalarRange(S)*Rough+1); //10
#local P[2][0][2]=P[1][1][1];

#local P[3][1][1]=<-qB,0,-qA>*Size*(ScalarRange(S)*Rough+1); //11
#local P[4][0][2]=P[3][1][1];

//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 make room for the centerpoints of each edge

#local SC=SC+SC;

#local I=SC;   #while (I>-1)
#local J=SC*2; #while (J>-1)
#local K=0; #while (K<5)

  #local P[K][I][J]=P[K][I/2][J/2];

#local K=K+1; #end
#local J=J-2; #end
#local I=I-2; #end

// logic to calculate each new mid point

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

#local K=0;
#while (K<5)

#local J=0;
#while (J<SC*2)

#if (Rounding=false)
  RoughMid(P,K,I,J, I+2,J, I+1,J)
  RoughMid(P,K,I,J, I,J+2, I,J+1)
  RoughMid(P,K,I+2,J, I,J+2, I+1,J+1)
#else
  RoundedMid(P,K,I,J, I+2,J, I+1,J)
  RoundedMid(P,K,I,J, I,J+2, I,J+1)
  RoundedMid(P,K,I+2,J, I,J+2, I+1,J+1)
#end

#local J=J+2;
#end

#if (Rounding=false)
  RoughMid(P,K,I,J, I+2,J, I+1,J)
#else
  RoundedMid(P,K,I,J, I+2,J, I+1,J)
#end

#local K=K+1;
#end

#local I=I+2;
#end

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

#local P[0][I][0]    =P[1][0][I];
#local P[0][SC][I]   =P[1][0][SC+I];
#local P[0][SC][SC+I]=P[1][I][SC+SC];

#local P[1][I][0]    =P[2][0][I];
#local P[1][SC][I]   =P[2][0][SC+I];
#local P[1][SC][SC+I]=P[2][I][SC+SC];

#local P[2][I][0]    =P[3][0][I];
#local P[2][SC][I]   =P[3][0][SC+I];
#local P[2][SC][SC+I]=P[3][I][SC+SC];

#local P[3][I][0]    =P[4][0][I];
#local P[3][SC][I]   =P[4][0][SC+I];
#local P[3][SC][SC+I]=P[4][I][SC+SC];

#local P[4][I][0]    =P[0][0][I];
#local P[4][SC][I]   =P[0][0][SC+I];
#local P[4][SC][SC+I]=P[0][I][SC+SC];

#local I=I+1; #end

#local Depth=Depth+1;
#end

// Conversion done to here

#if (Smooth!=false)

//Now the triangle normals get calculated

// clear out the normals for each vertex
#local I=0; #while (I<SC+1)
#local J=0; #while (J<SC+SC+1)
#local K=0; #while (K<5)

#local N[K][I][J]=<0,0,0>;

#local K=K+1; #end
#local J=J+1; #end
#local I=I+1; #end

//Now we calculate the surface normals for each triangle

#local I=0; #while (I<SC)
#local J=0; #while (J<SC+SC)
#local K=0; #while (K<5)

#local V1=P[K][I+1][J]-P[K][I][J];
#local V2=P[K][I][J+1]-P[K][I][J];
#local N[K][I][J]=N[K][I][J]+vcross(V1,V2);
#local N[K][I+1][J]=N[K][I+1][J]+vcross(V1,V2);
#local N[K][I][J+1]=N[K][I][J+1]+vcross(V1,V2);

#local V1=P[K][I][J+1]-P[K][I+1][J+1];
#local V2=P[K][I+1][J]-P[K][I+1][J+1];
#local N[K][I+1][J+1]=N[K][I+1][J+1]+vcross(V1,V2);
#local N[K][I+1][J]=N[K][I+1][J]+vcross(V1,V2);
#local N[K][I][J+1]=N[K][I][J+1]+vcross(V1,V2);

#local K=K+1; #end
#local J=J+1; #end
#local I=I+1; #end

#local I=1; #while (I<SC)

#local N[0][I][0]=N[0][I][0]+N[1][0][I];
#local N[1][0][I]=N[0][I][0];

#local N[0][SC][I]=  N[0][SC][I]+N[1][0][SC+I];
#local N[1][0][SC+I]=N[0][SC][I];

#local N[0][SC][SC+I]=N[0][SC][SC+I]+N[1][I][SC+SC];
#local N[1][I][SC+SC]=N[0][SC][SC+I];

#local N[1][I][0]=N[1][I][0]+N[2][0][I];
#local N[2][0][I]=N[1][I][0];

#local N[1][SC][I]=  N[1][SC][I]+N[2][0][SC+I];
#local N[2][0][SC+I]=N[1][SC][I];

#local N[1][SC][SC+I]=N[1][SC][SC+I]+N[2][I][SC+SC];
#local N[2][I][SC+SC]=N[1][SC][SC+I];

#local N[2][I][0]=N[2][I][0]+N[3][0][I];
#local N[3][0][I]=N[2][I][0];

#local N[2][SC][I]=  N[2][SC][I]+N[3][0][SC+I];
#local N[3][0][SC+I]=N[2][SC][I];

#local N[2][SC][SC+I]=N[2][SC][SC+I]+N[3][I][SC+SC];
#local N[3][I][SC+SC]=N[2][SC][SC+I];

#local N[3][I][0]=N[3][I][0]+N[4][0][I];
#local N[4][0][I]=N[3][I][0];

#local N[3][SC][I]=  N[3][SC][I]+N[4][0][SC+I];
#local N[4][0][SC+I]=N[3][SC][I];

#local N[3][SC][SC+I]=N[3][SC][SC+I]+N[4][I][SC+SC];
#local N[4][I][SC+SC]=N[3][SC][SC+I];

#local N[4][I][0]=N[4][I][0]+N[0][0][I];
#local N[0][0][I]=N[4][I][0];

#local N[4][SC][I]=  N[4][SC][I]+N[0][0][SC+I];
#local N[0][0][SC+I]=N[4][SC][I];

#local N[4][SC][SC+I]=N[4][SC][SC+I]+N[0][I][SC+SC];
#local N[0][I][SC+SC]=N[4][SC][SC+I];

#local I=I+1; #end

#local N[0][0][0]=N[0][0][0]+N[1][0][0]+N[2][0][0]+N[3][0][0]+N[4][0][0];
#local N[1][0][0]=N[0][0][0];
#local N[2][0][0]=N[0][0][0];
#local N[3][0][0]=N[0][0][0];
#local N[4][0][0]=N[0][0][0];

#local N[0][SC][SC*2]=N[0][SC][SC*2]+N[1][SC][SC*2]+N[2][SC][SC*2]+N[3][SC][SC*2]+N[4][SC][SC*2];
#local N[1][SC][SC*2]=N[0][SC][SC*2];
#local N[2][SC][SC*2]=N[0][SC][SC*2];
#local N[3][SC][SC*2]=N[0][SC][SC*2];
#local N[4][SC][SC*2]=N[0][SC][SC*2];

#local N[0][SC][0]=N[0][SC][0]+N[1][0][SC];
#local N[1][0][SC]=N[0][SC][0];
#local N[1][SC][0]=N[1][SC][0]+N[2][0][SC];
#local N[2][0][SC]=N[1][SC][0];
#local N[2][SC][0]=N[2][SC][0]+N[3][0][SC];
#local N[3][0][SC]=N[2][SC][0];
#local N[3][SC][0]=N[3][SC][0]+N[4][0][SC];
#local N[4][0][SC]=N[3][SC][0];
#local N[4][SC][0]=N[4][SC][0]+N[0][0][SC];
#local N[0][0][SC]=N[4][SC][0];

#local N[0][SC][SC]=N[0][SC][SC]+N[1][0][SC+SC];
#local N[1][0][SC+SC]=N[0][SC][SC];
#local N[1][SC][SC]=N[1][SC][SC]+N[2][0][SC+SC];
#local N[2][0][SC+SC]=N[1][SC][SC];
#local N[2][SC][SC]=N[2][SC][SC]+N[3][0][SC+SC];
#local N[3][0][SC+SC]=N[2][SC][SC];
#local N[3][SC][SC]=N[3][SC][SC]+N[4][0][SC+SC];
#local N[4][0][SC+SC]=N[3][SC][SC];
#local N[4][SC][SC]=N[4][SC][SC]+N[0][0][SC+SC];
#local N[0][0][SC+SC]=N[4][SC][SC];

#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<SC)
#local J=0; #while (J<SC+SC)
#local K=0; #while (K<5)

#if (Smooth=false)
  triangle { P[K][I][J], P[K][I][J+1], P[K][I+1][J] }
  triangle { P[K][I+1][J+1], P[K][I+1][J], P[K][I][J+1] }
#else
  smooth_triangle { P[K][I][J], N[K][I][J],
    P[K][I+1][J], N[K][I+1][J], P[K][I][J+1], N[K][I][J+1] }
  smooth_triangle { P[K][I+1][J+1], N[K][I+1][J+1],
    P[K][I][J+1], N[K][I][J+1], P[K][I+1][J], N[K][I+1][J] }
#end

#local K=K+1; #end
#local J=J+1; #end
#local I=I+1; #end
}