#declare ag = <0,-9.8,0>;		// Gravitational acceleration, m/s^2
#declare rhoa = 1.2; // Air density, kg/m^3
#declare rhos = 4542; // Lava density, kg/m^3
#declare deltat = 0.01; // Time differential
#declare alt = y*3500; // Initial launch position

#declare numShot = 10000;	// Number of blobs
#declare aRange = <65, 89, 3>; // azimuth Range
#declare tRange = <0, 160, 3>; // time range
#declare vRange = <150, 950, 4>; // velocity range
#declare rRange = <35, 65, 5>; // radius range
#declare rSeed = seed(4985); // Random number seed

#macro RangeNum(rng,sd,bias)
	/* Generates random number from range defined as vector
		rng.x -- range minimum
		rng.y -- range maximum
		rng.z -- number of samples, higher values make a bell curve
	*/
	#if (bias<.001 | bias>1000) #error "RangeNum: bias out of range. confine to 1e-3 to 1e3.\n" #end
	#local sep = rng.y-rng.x;
	#local i=0;
	#local trg=0;
	#while (i<rng.z) #local	trg = trg + rand(sd); #local i = i + 1; #end
	#local trg = pow(trg/rng.z,1/bias);
	#local ans = rng.x+sep*trg;
	ans
#end

#macro Volcano(gnd,oFile)
	#declare vPos = array[numShot]
	#declare vVel = array[numShot]
	#declare vRad = array[numShot]
	#declare vTim = array[numShot]
	#declare vImp = array[numShot]
	
	#local i = 0;
	#debug "Initializing... \n"
	#while (i<numShot)
		#declare vPos[i] = alt; // Initial position
		#declare vImp[i] = false; // True if surface impact detected 
		#declare vRad[i] = RangeNum(rRange,rSeed, .5); // Radius
		#local j = 0;
		#local tMax = RangeNum(tRange,rSeed, 1.2);
		#while (defined(vTim[j])) // Keep flight times sorted
			#if (vTim[j]>tMax)
				#local sw = tMax;
				#local tMax = vTim[j];
				#declare vTim[j] = sw;
			#end
			#local j = j + 1;
		#end
		#declare vTim[i] = tMax; // Flight time
		#local Azim = radians(RangeNum(aRange,rSeed,1));
		#local Lngt = RangeNum(<0, 2*pi, 1>,rSeed,1);
		#local iSpd = RangeNum(vRange,rSeed,1.3);
		#local iDir = <cos(Lngt)*cos(Azim), sin(Azim), sin(Lngt)*cos(Azim)>;
		#declare vVel[i] = iDir*iSpd;
		#local i = i + 1;
		#if (mod(i,25)=0) #debug concat("Particle #",str(i,0,0)," of ",str(numShot,0,0)," initiaized.\n") #end
	#end
	
	#debug "Launching particles... \n"
	#local Shot0 = 0; // Starting particle
	#local cTim = 0; // Current time
	#declare gNrm = <0,0,0>;
	#while (Shot0<numShot)
		#local cShot = Shot0;
		#declare cTim = cTim + deltat;
		#while (cShot<numShot)
			#if (vImp[cShot])
				// Impact has occurred, particle does not move
			#else
				// Move particle
				#local rd = vRad[cShot];
				#local spd = vlength(vVel[cShot]); // Speed, m/s
				#local dir = vnormalize(vVel[cShot]); // Direction
				#local rdFlat = 1+pow(spd,0.25)*.2; // Drag stretching factor
				#local Cd = 1.3/rdFlat/rdFlat; // Drag coefficient
				#local sMass = 4/3*pi*rd*rd*rd*rhos; // Mass, kg
				#local sCross = pi*rd*rd/rdFlat; // Cross section, m^2
				
				#local sAcc = ag - dir*0.5*Cd*rhoa*sCross*spd*spd/sMass;
				#local vo = vVel[cShot];
				#declare vVel[cShot] = vVel[cShot] + deltat*sAcc;
				#declare vPos[cShot] = vPos[cShot] + deltat*(vo+vVel[cShot])/2;
				#local gAlt = trace(gnd, vPos[cShot]+y*1e6, -y, gNrm);
				
				#if (gAlt.y>=vPos[cShot].y) // Impact, set flag
					#declare vImp[cShot] = true;
					#debug concat("Particle #",str(cShot,0,0)," impacted after ",str(cTim,0,3)," seconds.\n")
					#declare vPos[cShot] = gAlt;
					#declare vTim[cShot] = 0;
				#end
			#end	
			// Check for expired particle
			#if (vTim[cShot]<cTim)
				#if (mod(cShot,10)=0) #debug concat("Particle #",str(Shot0,0,0)," expired after ",str(cTim,0,3)," seconds.\n") #end
				#declare Shot0 = Shot0 + 1;
			#end
			#local cShot = cShot + 1;
		#end
	#end
	#local i=0;
	#fopen pFile oFile write
	#while (i<numShot)
		#if (vImp[i]) // On impact, blob is flattened disc
			#declare rdFlat = 0.5;
			#declare rdRot = <90,0,0>;
		#else // Flyng blob is stretched by part of the drag force
			#local spd = vlength(vVel[i]); // Speed, m/s
			#local dir = vnormalize(vVel[i]); // Direction
			#declare rdFlat = 1+pow(spd,0.25)*.2; // Drag stretching factor
			#declare rdRot = <degrees(asin(dir.y)), degrees(atan2(-dir.x,dir.z)), 0>;
		#end
		#local rdScale = <pow(rdFlat, -.5), pow(rdFlat, -.5), rdFlat>;
		
		#write( pFile, "sphere { \n  0,",vRad[i]," scale ",rdScale," rotate ",rdRot," translate ",vPos[i], "\n  hollow pigment { rgbt 1 } interior { intLava } }\n")
		#local i = i + 1;
		#if (mod(i,25)=0) #debug concat("Particle #",str(i,0,0)," of ",str(numShot,0,0)," in place.\n") #end
	#end
	#fclose pFile
#end
		