//based on code from
//http://www.javaworld.com/javaworld/jw-08-1998/jw-08-step_p.html
//left and right shift: http://www.janeg.ca/scjp/oper/shift.html

//macro takes density and roughness and seed and "falloff"
//Density is an integer from 1 on up, but 5 is a good
//test value, 6-9 probably a good value, but slow as hell
//Roughness ranges from 0 to 1, with zero meaning no random roughness.
//Roughnesses higher than 1 should work, but will produce roughnesses larger than the
//mountain peak.
//Seed is just the seed to the random number generator. Should allow you to
//repeat your mountains exactly, or choose new ones
//xPlateau is 0 for a mountain peak, otherwise is a mountain range extending left-right;
//the higher the value, the wider the range before it begins to fall off
//Plateau values can range from 0 to 1.
//Damp Types determine the falloff speed
#declare MT_PEAK = 0;
#declare MT_RANGE = .7;
#declare MT_LINEAR = 2;
#declare MT_SINE = 3;
#declare MT_SQUARE = 5;
#macro mountains(mr_Density, mr_Roughness, mr_Seed, mr_xPlateau, mr_DampType)
	//create the altitudes
	#local Divisions = pow(2,mr_Density);
	#local Terrain = array[Divisions+1][Divisions+1]
	#local RNG = seed(mr_Seed);

	#declare Terrain[0][0] = rand(RNG);
	#declare Terrain[0][Divisions] = rand(RNG);
	#declare Terrain[Divisions][Divisions] = rand(RNG);
	#declare Terrain[Divisions][0] = rand(RNG);

	#local Rough = mr_Roughness;
	#local I = 0;
	#while (I < mr_Density)
		#local Q = pow(2, I);
		#local R = pow(2, mr_Density - I);
		#local S = R / pow(2, 1);
		#local J = 0;
		#while (J < Divisions)
			#local K = 0;
			#while (K < Divisions)
				//diamond (j, k, r, rough);
				#if (R > 1)
					#local Half = R / 2;
					#local AVG = (Terrain[J][K] + Terrain[J+R][K] +
						Terrain[J+R][K+R] + Terrain[J][K+R])/4;
					#declare Terrain[J+Half][K+Half] = AVG + rand(RNG) * Rough;
				#end
				#declare K = K + R;
			#end
			#declare J = J + R;
		#end
		#if (S > 0)
			#local J = 0;
			#while (J <= Divisions)
				#declare K = mod(J + S, R); //this might be reversed
				#while (K <= Divisions)
					//square (j - s, k - s, r, rough);
					#local Xjs = J - S;
					#local Zks = K - S;
					#local Half = R / 2;
					#local AVG = 0;
					#local SUM = 0;
					#if (Xjs >= 0)
						#declare AVG = AVG + Terrain[Xjs][Zks+Half];
						#declare SUM = SUM + 1;
					#end
					#if (Zks >= 0)
						#declare AVG = AVG + Terrain[Xjs + Half][Zks];
						#declare SUM = SUM + 1;
					#end
					#if (Xjs + R <= Divisions)
						#declare AVG = AVG + Terrain[Xjs + R][Zks + Half];
						#declare SUM = SUM + 1;
					#end
					#if (Zks + R <= Divisions)
						#declare AVG = AVG + Terrain[Xjs + Half][Zks + R];
						#declare SUM = SUM + 1;
					#end
					#declare Terrain[Xjs + Half][Zks + Half] = AVG/SUM + rand(RNG)*Rough;
					#declare K = K + R;
				#end
				#declare J = J + S;
			#end
		#end
		#declare Rough = Rough*mr_Roughness;
		#declare I = I + 1;
	#end

	//damp the edges to zero and determine max height
	#local Max = 0;
	#local Ix = 0;
	#declare xCenter = Divisions/2;
	#declare zCenter = Divisions/2;
	#while (Ix <= Divisions)
		#declare xFromCenter = abs(Ix - xCenter)/xCenter;
		#if (xFromCenter <= mr_xPlateau)
			#declare xFromCenter = 0;
		#end
		#local Iz = 0;
		#while (Iz <= Divisions)
			#local zFromCenter = abs(Iz - zCenter)/zCenter;
			#local Height = 1-sqrt(zFromCenter*zFromCenter+xFromCenter*xFromCenter);
			#if (Height < 0)
				#declare Height = 0;
			#end
			#if (int(mr_DampType/MT_SINE) = mr_DampType/MT_SINE)
				#declare Height = (sin(radians(Height*180-90))+1)/2;
			#end
			#if (int(mr_DampType/MT_SQUARE) = mr_DampType/MT_SQUARE)
				#declare Height = Height*Height;
			#end
			#declare Terrain[Ix][Iz] = Terrain[Ix][Iz]*Height;
			
			//determine maximum
			#if (Terrain[Ix][Iz] > Max)
				#declare Max = Terrain[Ix][Iz];
			#end
			#declare Iz = Iz + 1;
		#end
		#declare Ix = Ix + 1;
	#end
	
	mesh {
		#declare mr_x = 0;
		#while (mr_x < Divisions)
			#declare mr_z = 0;
			#while (mr_z < Divisions)
				//draw the square
				#if (Terrain[mr_x][mr_z] = 0 & Terrain[mr_x+1][mr_z] = 0 & Terrain[mr_x+1][mr_z+1] = 0)
					//we don't need this triangle
				#else
					triangle {
						<mr_x, Terrain[mr_x][mr_z], mr_z>, 
						<mr_x+1, Terrain[mr_x+1][mr_z], mr_z>
						 <mr_x+1, Terrain[mr_x+1][mr_z+1], mr_z+1>
					}
				#end
				#if (Terrain[mr_x][mr_z] = 0 & Terrain[mr_x+1][mr_z+1] = 0 & Terrain[mr_x][mr_z+1] = 0)
					//we don't need this triangle
				#else
					triangle {
						<mr_x, Terrain[mr_x][mr_z], mr_z>,
						<mr_x+1, Terrain[mr_x+1][mr_z+1], mr_z+1>,
						<mr_x, Terrain[mr_x][mr_z+1], mr_z+1>
					}
				#end
				#declare mr_z = mr_z + 1;
			#end
			#declare mr_x = mr_x + 1;
		#end
		scale <1/Divisions, 1/Max, 1/Divisions>
	}
#end

#macro mountain_peak(mp_Density, mp_Roughness, mp_Seed)
	mountains(mp_Density, mp_Roughness, mp_Seed, MT_PEAK, MT_LINEAR)
	translate <-.5, 0, -.5>
#end

#macro foothill(mp_Density, mp_Roughness, mp_Seed)
	mountains(mp_Density, mp_Roughness, mp_Seed, MT_PEAK, MT_SINE*MT_SQUARE)
	translate <-.5, 0, -.5>
#end

#macro mountain_range(mp_Density, mp_Roughness, mp_Seed, mp_Length)
	mountains(mp_Density, mp_Roughness, mp_Seed, MT_RANGE, MT_LINEAR)
	translate <-.5, 0, -.5>
	scale <mp_Length, 1, 1>
#end


#declare groundGreen = texture {
	pigment {
		bumps
		turbulence .25
		color_map {
			[0 rgb <1, 1, 0>]
			[1 rgb <.1, .95, 0>]
		}
	}
	normal {
		bumps .25
		turbulence .25
	}
	finish {
		ambient .35
	}
}

#declare mountainGray = texture {
	pigment {
		color White*.8
	}
	
	normal {
		wrinkles 2
		scale .1
	}
}

#declare mountainSnow = texture {
	pigment {
		color White
	}
	
	normal {
		wrinkles 2
		scale .1
	}
	
	finish {
		ambient .4
	}
}

#declare mountainTexture = 	texture {
	gradient y
	turbulence .4
	texture_map {
		[0.0	groundGreen]
		[0.4	groundGreen]
		[0.6	mountainGray]
		[0.8	mountainGray]
		[0.9	mountainSnow]
		[1.0	mountainSnow]		
	}
}