#include "math.inc"
#include "transforms.inc"

#local exagg_scale	= 1000;		// the factor to exaggerate the planets by
#local exagg_incln	= 1;		// the factor to exaggerate the planets by
#local exagg_eccty	= 1;		// the factor to exaggerate the planets by
#local markr_thick	= 1/200;
#local AU		= 149597870.691;


// ----------------------------------------

global_settings
{
	assumed_gamma 1.0
}
default {finish {ambient 1/2 diffuse 1}}
background {color rgb 1/2}


// ----------------------------------------

#local temp_camera_1 = camera
{
	orthographic
	location	z * -12
	direction	z
	up		y
	right		x * image_width/image_height
	rotate		x * 90
	rotate		y * 0
	scale		3
}

#local temp_light_1 = light_source
{
	0
	color rgb 1
	shadowless
}
#local temp_light_2 = light_source
{
	x+y+z
	color rgb 1
	parallel
}


// ----------------------------------------

union
{
//	object {temp_light_1}
	object {temp_light_2}
}

camera {temp_camera_1}

cylinder
{
	-x * 10, +x * 10, markr_thick
	pigment {color rgb x}
	no_shadow
}
cylinder
{
	-y * 10, +y * 10, markr_thick
	pigment {color rgb y}
	no_shadow
}
cylinder
{
	-z * 10, +z * 10, markr_thick
	pigment {color rgb z}
	no_shadow
}

plane
{
	y,0
	pigment
	{
		onion
		color_map
		{
			[0.00 color rgb 0	transmit 1/2]
			[0.01 color rgb 0	transmit 1/2]
			[0.01 color rgb (x + z)	transmit 1/2]
			[0.99 color rgb (x + z)	transmit 1/2]
			[0.99 color rgb 0	transmit 1/2]
			[1.00 color rgb 0	transmit 1/2]
		}
	}
//	no_shadow
}


// ----------------------------------------

#local SunMeanRadius	= 696000/AU;
#local SunObliquity	= 7.15;		// decimal part may be minutes

sphere
{
	0, SunMeanRadius * exagg_scale/100
	pigment {color rgb <1,1,0>}
	no_shadow
}

#macro hhmmssToDecimal(in_vec)
	#local mns = in_vec.z/60;
	#local hrs = (in_vec.y + mns)/60;
	#local day = (in_vec.x + hrs)/24;
	day
#end

#macro TimeToAngle(in_time,in_rate)
	in_time/in_rate * 360
#end

#macro m_KeplerPosition_trans(i, Omega, w, epsilon, a, P, tau, tol, t_clock)
	transform
	{
		translate a * m_orbitalPosition
		(
			epsilon,
			2 * pi * (t_clock-tau)/P,
			2 * pi * tol/P
		)
		rotate <0, -w, 0>
		rotate <-i, -Omega, 0>
	}
#end

#macro planet_pigment(in_color)
	pigment
	{
		radial
		color_map
		{
			[0/8 color rgb 1]
			[1/8 color rgb 1]
			[1/8 color rgb in_color]
			[3/8 color rgb in_color]
			[3/8 color rgb 0]
			[5/8 color rgb 0]
			[5/8 color rgb in_color]
			[7/8 color rgb in_color]
			[7/8 color rgb 1]
			[9/8 color rgb 1]
		}
	}
#end

#declare cylinder_pigment = pigment
{
	gradient y
	color_map
	{
		[0/2 color rgb 0]
		[1/2 color rgb 0]
		[1/2 color rgb 1]
		[2/2 color rgb 1]
	}
	translate	-y/2
	scale		1.0001
}


#include "fbruder_kepler.inc"


#local Today		= 2455376.500000000;
#local NumPlanets	= 4;
#local PlanetName	= array[NumPlanets]{            "Mercury",              "Venus",              "Earth",               "Mars"}	// string
#local MeanRadius	= array[NumPlanets]{           2440.00/AU,           6051.80/AU,           6371.01/AU,           3389.90/AU}	// AU
#local SemiMajorAxis	= array[NumPlanets]{3.870979296576140E-01,7.233286735377678E-01,9.992400161311368E-01,1.523739485637696E+00}	// AU
#local Eccentricity	= array[NumPlanets]{2.056242542315567E-01,6.810973251609722E-03,1.745257128512940E-02,9.331215215138063E-02}	// deg
#local AscendingNode	= array[NumPlanets]{4.831749863564502E+01,7.665036452914353E+01,1.251068496744354E+02,4.952430761299259E+01}	// deg
#local ArgOfPerihelion	= array[NumPlanets]{2.915316085559863E+01,5.488620905888749E+01,3.381688838185740E+02,2.865885963233656E+02}	// deg
#local Inclination	= array[NumPlanets]{7.004365876146918E+00,3.394489838966820E+00,2.652965379014084E-03,1.848898333793175E+00}	// deg
#local SidOrbitPeriod	= array[NumPlanets]{8.796900173442556E+01,2.246991436826062E+02,3.648400455253299E+02,6.870121845099741E+02}	// days
#local PigmentColor	= array[NumPlanets]{    <115,115,115>/255,    <222,199,165>/255,    <028,083,108>/255,    <135,100,060>/255}	// rgb color
#local DateOfPeriapsis	= array[NumPlanets]{ 2455372.937552663963, 2455333.428305923007, 2455200.636899325531, 2455630.005895495880}	// days
#local SidRotatPeriod	= array[NumPlanets]{              58.6462,            -243.0185,          23.93419/24,         24.622962/24}	// days, sidereal
#local Obliqueness	= array[NumPlanets]{                 2.11,                177.3,                23.45,                25.19}	// deg, make sure decimal parts are actual decimals, not minutes
#local TimeAtPeriapsis	= array[NumPlanets]{      <05,24,41.9530>,      <20,46,56.5818>,      <10,11,40.3137>,      <06,47,55.9578>}	// HH,MM,SS.FFFF


// Returns the position of a planet with given orbital elements at the current time (given by clock).
// Parameters:
// i     - inclination
// Omega - longitude of the ascending node, whereby the positive x direction is the vernal point
// w     - argument of perihelion
// epsilon - eccentricity
// a     - semi-major axis
// P     - orbital period
// tau   - time of periapsis passage
// tol   - tolerance of accuracy. This is an accuracy in time.
//         The returned position lies on the orbit and is passed at a time T with
//         (clock - tol) < T < (clock + tol).
//         A reasonable value for tol is clock_delta/2.
//         The number of iterations for solving Kepler's Equation will be time-independant. So the motion
//         will always be smooth.

// still need to know the starting angle of each planet WRT its obliqueness
union
{
	#local Iter = 0;
	#while (Iter < NumPlanets)
		#local Radius = MeanRadius[Iter] * exagg_scale;
		#local Planet_Pos = m_KeplerPosition
		(
			Inclination[Iter] * exagg_incln,
			AscendingNode[Iter],
			ArgOfPerihelion[Iter],
			Eccentricity[Iter] * exagg_eccty,
			SemiMajorAxis[Iter],
			SidOrbitPeriod[Iter],
			DateOfPeriapsis[Iter] - Today,
			1/2,
			clock
		);
		#local Marker_Pos1 = m_KeplerPosition
		(
			Inclination[Iter] * exagg_incln,
			AscendingNode[Iter],
			ArgOfPerihelion[Iter],
			Eccentricity[Iter] * exagg_eccty,
			SemiMajorAxis[Iter],
			1,
			0,
			1/2,
			0
		);
		#local Marker_Pos2 = m_KeplerPosition
		(
			Inclination[Iter] * exagg_incln,
			AscendingNode[Iter],
			ArgOfPerihelion[Iter],
			Eccentricity[Iter] * exagg_eccty,
			SemiMajorAxis[Iter],
			1,
			0,
			1/2,
			1/2
		);
		#local Marker_Pos3 = m_KeplerPosition
		(
			Inclination[Iter] * exagg_incln,
			AscendingNode[Iter],
			ArgOfPerihelion[Iter],
			Eccentricity[Iter] * exagg_eccty,
			SemiMajorAxis[Iter],
			1,
			0,
			1/2,
			1/4
		);
		#local Marker_Pos4 = m_KeplerPosition
		(
			Inclination[Iter] * exagg_incln,
			AscendingNode[Iter],
			ArgOfPerihelion[Iter],
			Eccentricity[Iter] * exagg_eccty,
			SemiMajorAxis[Iter],
			1,
			0,
			1/2,
			3/4
		);
		#local vector_a = vnormalize(Marker_Pos2-Marker_Pos1);
		#local vector_b = vnormalize(Marker_Pos4-Marker_Pos3);
		#local vector_c = VPerp_To_Plane(vector_b,vector_a);

		// maybe use "Solar sub-long & sub-lat" for obliqueness? dunno
		#local rotate_rate		= SidRotatPeriod[Iter];
		#local start_local_time		= hhmmssToDecimal(TimeAtPeriapsis[Iter]);
		#local start_angle		= start_local_time * 360;
		#local time_since_start		= Today - DateOfPeriapsis[Iter];
		#local angle_since_start	= TimeToAngle(time_since_start,rotate_rate);
		#local angle_now		= start_angle + angle_since_start;
		#debug concat("\nrotate_rate=",		str(rotate_rate,0,-1),"\n")
		#debug concat("\nstart_local_time=",	str(start_local_time,0,-1),"\n")
		#debug concat("\nstart_angle=",		str(start_angle,0,-1),"\n")
		#debug concat("\ntime_since_start=",	str(time_since_start,0,-1),"\n")
		#debug concat("\nangle_since_start=",	str(angle_since_start,0,-1),"\n")
		#debug concat("\nangle_now=",		str(angle_now,0,-1),"\n")
		union
		{
			sphere
			{
				0,1
				planet_pigment(PigmentColor[Iter])
			}
			cylinder
			{
				y*-2,y*+2,1/4
				pigment {cylinder_pigment scale 4}
			}
/*
			cylinder
			{
				-x*16, x*16, 1/8
				pigment {cylinder_pigment scale 32 rotate z * -90}
			}
*/
//			scale		Radius
			scale		1/20
			rotate		y * -angle_now			// should be negative, as that is the direction the planets rotate
			Point_At_Trans(vector_c)			// make sure no planets were flipped upside down
			translate	Planet_Pos
//			no_shadow
		}
		union
		{
//			cylinder {Marker_Pos1, Marker_Pos2, markr_thick}
//			cylinder {Marker_Pos3, Marker_Pos4, markr_thick}
			sphere_sweep
			{
				#local sweep_count = 0;
				#local sweep_total = 24;
				linear_spline
				sweep_total
				#while (sweep_count < sweep_total)				// change "<" to "<=" to close the gap
					m_KeplerPosition
					(
						Inclination[Iter] * exagg_incln,
						AscendingNode[Iter],
						ArgOfPerihelion[Iter],
						Eccentricity[Iter] * exagg_eccty,
						SemiMajorAxis[Iter],
						1,
						DateOfPeriapsis[Iter] - Today,
						1/2,
						sweep_count/sweep_total * 1
					),
					markr_thick
					#local sweep_count = sweep_count + 1;
				#end
	//			no_shadow
			}
			pigment {color rgb 1}
		}
		#local Iter = Iter + 1;
	#end
}