Ephemeris
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © BullRider802

//@version=5

// @description TODO: add library description here
library("Ephemeris")

// ========================================================================================================================
//                                               PLANETARY KEPLER CONSTANTS
// ========================================================================================================================

//ELEMENTS @ J2000: a, e, i, mean longitude (L), longitude of perihelion, longitude of ascending node
//https://ssd.jpl.nasa.gov/?planet_pos
// a = semi major axis [0]
// e = eccentricity [1]
// i = orbit inclination [2]
// L = mean longitude [3]
// w = argument of perihelion [4]
// N = longitude of ascending node [5]

// Mercury Kepler Elements
// @returns Returns constant Mercury orbital elements used to calculate Kepler equations
export mercuryElements() =>
	mercuryElements = array.new<float>(6, 0.0)
	array.insert(mercuryElements, 0, 0.38709927)
	array.insert(mercuryElements, 1, 0.20563593)
	array.insert(mercuryElements, 2, 7.00497902)
	array.insert(mercuryElements, 3, 252.25032350)
	array.insert(mercuryElements, 4, 77.45779628)
	array.insert(mercuryElements, 5, 48.33076593)
	mercuryElements

// Mercury Kepler Rates
// @returns Returns constant Mercury orbital rates used to calculate Kepler equations
export mercuryRates() =>
	mercuryRates = array.new<float>(6, 0.0)
	array.insert(mercuryRates, 0, 0.00000037)
	array.insert(mercuryRates, 1, 0.00001906)
	array.insert(mercuryRates, 2, -0.00594749)
	array.insert(mercuryRates, 3, 149472.67411175)
	array.insert(mercuryRates, 4, 0.16047689)
	array.insert(mercuryRates, 5, -0.12534081)
	mercuryRates

// Venus Kepler Elements
// @returns Returns constant Venus orbital elements used to calculate Kepler equations
export venusElements() =>
	venusElements = array.new<float>(6, 0.0)
	array.insert(venusElements, 0, 0.72333566)
	array.insert(venusElements, 1, 0.00677672)
	array.insert(venusElements, 2, 3.39467605)
	array.insert(venusElements, 3, 181.97909950)
	array.insert(venusElements, 4, 131.60246718)
	array.insert(venusElements, 5, 76.67984255)
	venusElements

// Venus Kepler Rates
// @returns Returns constant Venus orbital rates used to calculate Kepler equations
export venusRates() =>
	venusRates = array.new<float>(6, 0.0)
	array.insert(venusRates, 0, 0.00000390)
	array.insert(venusRates, 1, -0.00004107)
	array.insert(venusRates, 2, -0.00078890)
	array.insert(venusRates, 3, 58517.81538729)
	array.insert(venusRates, 4, 0.00268329)
	array.insert(venusRates, 5, -0.27769418)
	venusRates

// Earth Kepler Elements
// @returns Returns constant Earth orbital elements used to calculate Kepler equations
export earthElements() =>
	earthElements = array.new<float>(6, 0.0)
	array.insert(earthElements, 0, 1.00000261)
	array.insert(earthElements, 1, 0.01671123)
	array.insert(earthElements, 2, -0.00001531)
	array.insert(earthElements, 3, 100.46457166)
	array.insert(earthElements, 4, 102.93768193)
	array.insert(earthElements, 5, 0.0)
	earthElements

// Earth Kepler Rates
// @returns Returns constant Earth orbital rates used to calculate Kepler equations
export earthRates() =>
	earthRates = array.new<float>(6, 0.0)
	array.insert(earthRates, 0, -0.00000562)
	array.insert(earthRates, 1, -0.00004392)
	array.insert(earthRates, 2, -0.01294668)
	array.insert(earthRates, 3, 35999.37244981)
	array.insert(earthRates, 4, 0.32327364)
	array.insert(earthRates, 5, 0.0)
	earthRates

// Mars Kepler Elements
// @returns Returns constant Mars orbital elements used to calculate Kepler equations
export marsElements() =>
	marsElements = array.new<float>(6, 0.0)
	array.insert(marsElements, 0, 1.52371034)
	array.insert(marsElements, 1, 0.09339410)
	array.insert(marsElements, 2, 1.84969142)
	array.insert(marsElements, 3, -4.55343205)
	array.insert(marsElements, 4, -23.94362959)
	array.insert(marsElements, 5, 49.55953891)
	marsElements

// @returns Returns constant Mars orbital rates used to calculate Kepler equations
export marsRates() =>
	marsRates = array.new<float>(6, 0.0)
	array.insert(marsRates, 0, 0.00001847)
	array.insert(marsRates, 1, 0.00007882)
	array.insert(marsRates, 2, 0.00813131)
	array.insert(marsRates, 3, 19140.30268499)
	array.insert(marsRates, 4, 0.44441088)
	array.insert(marsRates, 5, -0.29257343)
	marsRates

// @returns Returns constant Jupiter orbital elements used to calculate Kepler equations
export jupiterElements() =>
	jupiterElements = array.new<float>(6, 0.0)
	array.insert(jupiterElements, 0, 5.20288700)
	array.insert(jupiterElements, 1, 0.04838624)
	array.insert(jupiterElements, 2, 1.30439695)
	array.insert(jupiterElements, 3, 34.39644051)
	array.insert(jupiterElements, 4, 14.72847983)
	array.insert(jupiterElements, 5, 100.47390909)
	jupiterElements

// @returns Returns constant Jupiter orbital rates used to calculate Kepler equations
export jupiterRates() => 
	jupiterRates = array.new<float>(6, 0.0)
	array.insert(jupiterRates, 0, -0.00011607)
	array.insert(jupiterRates, 1, -0.00013253)
	array.insert(jupiterRates, 2, -0.00183714)
	array.insert(jupiterRates, 3, 3034.74612775)
	array.insert(jupiterRates, 4, 0.21252668)
	array.insert(jupiterRates, 5, 0.20469106)
	jupiterRates

// @returns Returns constant Saturn orbital elements used to calculate Kepler equations
export saturnElements() => 
	saturnElements = array.new<float>(6, 0.0)
	array.insert(saturnElements, 0, 9.53667594)
	array.insert(saturnElements, 1, 0.05386179)
	array.insert(saturnElements, 2, 2.48599187)
	array.insert(saturnElements, 3, 49.95424423)
	array.insert(saturnElements, 4, 92.59887831)
	array.insert(saturnElements, 5, 113.66242448)
	saturnElements

// @returns Returns constant Saturn orbital rates used to calculate Kepler equations
export saturnRates() =>
	saturnRates = array.new<float>(6, 0.0)
	array.insert(saturnRates, 0, -0.00125060)
	array.insert(saturnRates, 1, -0.00050991)
	array.insert(saturnRates, 2, 0.00193609)
	array.insert(saturnRates, 3, 1222.49362201)
	array.insert(saturnRates, 4, -0.41897216)
	array.insert(saturnRates, 5, -0.28867794)
	saturnRates

// @returns Returns constant Uranus orbital elements used to calculate Kepler equations
export uranusElements() => 
	uranusElements = array.new<float>(6, 0.0)
	array.insert(uranusElements, 0, 19.18916464)
	array.insert(uranusElements, 1, 0.04725744)
	array.insert(uranusElements, 2, 0.77263783)
	array.insert(uranusElements, 3, 313.23810451)
	array.insert(uranusElements, 4, 170.95427630)
	array.insert(uranusElements, 5, 74.01692503)
	uranusElements

// @returns Returns constant Uranus orbital rates used to calculate Kepler equations
export uranusRates() =>
	uranusRates = array.new<float>(6, 0.0)
	array.insert(uranusRates, 0, 0.00196176)
	array.insert(uranusRates, 1, -0.00004397)
	array.insert(uranusRates, 2, -0.00242939)
	array.insert(uranusRates, 3, 428.48202785)
	array.insert(uranusRates, 4, 0.40805281)
	array.insert(uranusRates, 5, 0.04240589)
	uranusRates


// @returns Returns constant Neptune orbital elements used to calculate Kepler equations
export neptuneElements() =>
	neptuneElements = array.new<float>(6, 0.0)
	array.insert(neptuneElements, 0, 30.06992276)
	array.insert(neptuneElements, 1, 0.00859048)
	array.insert(neptuneElements, 2, 1.77004347)
	array.insert(neptuneElements, 3, -55.12002969)
	array.insert(neptuneElements, 4, 44.96476227)
	array.insert(neptuneElements, 5, 131.78422574)
	neptuneElements

// @returns Returns constant Neptune orbital rates used to calculate Kepler equations
export neptuneRates() =>
	neptuneRates = array.new<float>(6, 0.0)
	array.insert(neptuneRates, 0, 0.00026291)
	array.insert(neptuneRates, 1, 0.00005105)
	array.insert(neptuneRates, 2, 0.00035372)
	array.insert(neptuneRates, 3, 218.45945325)
	array.insert(neptuneRates, 4, -0.32241464)
	array.insert(neptuneRates, 5, -0.00508664)
	neptuneRates

// ========================================================================================================================
//                                               CONSTANTS & HELPER FUNCTIONS
// ========================================================================================================================

// @function Normalize degrees to within [0, 360)
// @param x degrees to be normalized
// @returns Normalized degrees
export rev360(float x) =>
	float angle = x - math.floor(x / 360.0) * 360.0
	if angle < 0.0
		angle := angle + 360.0
	angle

// @function Scale angle in degrees
// @param x Angle in degrees
// @returns Scaled angle in degrees
export scaleAngle(float longitude, float magnitude, int harmonic) =>
	longitude*magnitude + 360*harmonic

// @function Constant Julian days per century
// @returns 36525
export julianCenturyInJulianDays() => 
	36525

// @function Julian date on J2000 epoch start (2000-01-01)
// @returns 2451545.0
export julianEpochJ2000() => 
	2451545.0

// @function Mean obliquity of the ecliptic on J2000 epoch start (2000-01-01)
// @returns 23.43928
export meanObliquityForJ2000() => 
	23.43928

// @function Convert calendar date to Julian date
// @param Year calendar year as integer (e.g. 2018)
// @param Month calendar month (January = 1, December = 12)
// @param Day calendar day of month (e.g. January valid days are 1-31)
// @param Hour valid values 0-23
// @param Minute valid values 0-60
export getJulianDate(int Year, int Month, int Day, int Hour, int Minute) =>
    y = Year
    m = Month
    d = Day
    h = Hour
    min = Minute
	//Adjust if B.C.
	if y < 0
		y := y + 1
	//Adjust if JAN or FEB
	if m == 1 or m == 2
		y := y - 1
		m := m + 12
	//Calculate A & B; ONLY if date is equal or after 1582-Oct-15
	float A = math.floor(y/100) //A
	float B = 2 - A + math.floor(A/4) //B
	//Added Minute Accuracy
	((math.floor(365.25*y)) + (math.floor(30.6001*(m+1))) + d + (0.04166666666666666666666666666667*h) + (0.000694444444444444*min) + 1720994.5 + B)

// @function Centuries since Julian Epoch 2000-01-01
// @param date Julian date to conver to Julian centuries
// @param epoch_start Julian date of epoch start (e.g. J2000 epoch = 2451545)
// @returns Julian date converted to Julian centuries
export julianCenturies(float date, float epoch_start) =>
    (date - epoch_start)/36525.0

// @function Calculate Julian centuries since epoch J2000 (2000-01-01)
// @param julianDate Julian Date in days
// @returns Julian centuries since epoch J2000 (2000-01-01)
export julianCenturiesSinceEpochJ2000(float julianDate) =>
	(julianDate - julianEpochJ2000()) / julianCenturyInJulianDays()

// ========================================================================================================================
//                                              UTILITY FUNCTIONS
// ========================================================================================================================

// @function Specialized arctan function
// @param y radians
// @param x radians
// @returns special arctan of y/x
export atan2(float y, float x) =>
	if x > 0
		math.atan(y/x)
	else if x < 0
		if y >= 0
			math.atan(y/x) + math.pi
		else
			math.atan(y/x) - math.pi
	else
		if y > 0
			math.pi / 2
		else if y < 0
			1.0 * math.pi / 2
		else
			runtime.error("X and Y cannot both be zero to compute atan2")
			na

//TAKEN FROM: http://www.jgiesen.de/kepler/kepler.html
//CREDIT: Juergen Giesen
//Used to solve for E
// @function Compute eccentricity of the anomaly
// @param ec Eccentricity of Orbit
// @param m_param Mean Anomaly ?
// @param dp Decimal places to round to 
// @returns Eccentricity of the Anomaly
export eccAnom(float ec, float m_param, int dp) =>
    m = m_param
	// arguments:
	// ec=eccentricity, m=mean anomaly,
	// dp=number of decimal places

	float pi = math.pi
	float K = pi / 180.0
	int maxIter = 30
	int i = 0
	float delta = math.pow(10,-dp)
	float E = na
	float F = na

	m := m / 360.0
	m := 2.0 * pi * (m - math.floor(m))

	if ec < 0.8
	    E := m
	else
	    E := pi

	F := E - ec * math.sin(m) - m

	while math.abs(F) > delta and i < maxIter
		E := E - F / (1.0 - ec * math.cos(E))
		F := E - ec * math.sin(E) - m
		i := i + 1

	E := E / K

	math.round(E * math.pow(10,dp)) / math.pow(10,dp)

// @function Compute planetary ephemeris (longtude relative to Earth or Sun) on a Julian date
// @param TGen Julian Date
// @param planetElementId All planet orbital elements in an array. This index references a specific planet's elements.
// @param planetRatesId All planet orbital rates in an array. This index references a specific planet's rates.
// @returns [xGen, yGen, zGen, rGen] X,Y,Z ecliptic rectangular coordinates and R radius from reference body.
export planetEphemerisCalc(float TGen, float[] planetElementId, float[] planetRatesId) =>
    float element_a = array.get(planetElementId, 0)
    float element_e = array.get(planetElementId, 1)
    float element_i = array.get(planetElementId, 2)
    float element_L = array.get(planetElementId, 3)
    float element_w = array.get(planetElementId, 4)
    float element_N = array.get(planetElementId, 5)

    float rate_a = array.get(planetRatesId, 0)
    float rate_e = array.get(planetRatesId, 1)
    float rate_i = array.get(planetRatesId, 2)
    float rate_L = array.get(planetRatesId, 3)
    float rate_w = array.get(planetRatesId, 4)
    float rate_N = array.get(planetRatesId, 5)
	//--------------------------------------------------------------------------------------------
	//1.
	//ORBIT SIZE
	//AU (CONSTANT = DOESN'T CHANGE)
	aGen = element_a + (rate_a * TGen)
	//2.
	//ORBIT SHAPE
	//ECCENTRICITY (CONSTANT = DOESN'T CHANGE)
	eGen = element_e + (rate_e * TGen)
	//--------------------------------------------------------------------------------------------
	//3.
	//ORBIT ORIENTATION
	//ORBITAL INCLINATION (CONSTANT = DOESN'T CHANGE)
	iGen = element_i + (rate_i * TGen)
	iGen := rev360(iGen) // iGen%360
	//4.
	//ORBIT ORIENTATION
	//LONG OF ASCENDING NODE (CONSTANT = DOESN'T CHANGE)
	NGen = element_N + (rate_N * TGen)
	NGen := rev360(NGen) // NGen%360
	//5.
	//ORBIT ORIENTATION
	//LONGITUDE OF THE PERIHELION
	wGen = element_w + (rate_w * TGen)
	wGen := rev360(wGen) // wGen%360
	//--------------------------------------------------------------------------------------------
	//6.
	//ORBIT POSITION
	//MEAN LONGITUDE (DYNAMIC = CHANGES OVER TIME)
	LGen = element_L + (rate_L * TGen)
	LGen := rev360(LGen) // LGen%360

	//MEAN ANOMALY --> Use this to determine Perihelion (0 degrees = Perihelion of planet)
	MGen = rev360(LGen - wGen)

	//ECCENTRIC ANOMALY
	EGen = eccAnom(eGen, MGen, 6)

	//ARGUMENT OF TRUE ANOMALY
	trueAnomalyArgGen = math.sqrt((1+eGen) / (1-eGen)) * (math.tan(math.toradians(EGen)/2))

	//TRUE ANOMALY (DYNAMIC = CHANGES OVER TIME)
	K = math.pi / 180.0 //Radian converter variable
	float nGen = na
	if trueAnomalyArgGen < 0
		nGen := 2 * (math.atan(trueAnomalyArgGen)/K+180) //ATAN = ARCTAN = INVERSE TAN
	else
		nGen := 2 * (math.atan(trueAnomalyArgGen)/K)
	//--------------------------------------------------------------------------------------------

	//CALCULATE RADIUS VECTOR
	rGen = (aGen * (1 - math.pow(eGen,2)))/(1 + (eGen * (math.cos(math.toradians(nGen)))))

	//TAKEN FROM: http://www.stargazing.net/kepler/ellipse.html
	//CREDIT: Keith Burnett
	//Used to determine Heliocentric Ecliptic Coordinates (this part works!)
	xGen = rGen * (math.cos(math.toradians(NGen)) * math.cos(math.toradians(nGen+wGen-NGen)) - math.sin(math.toradians(NGen)) * math.sin(math.toradians(nGen+wGen-NGen)) * math.cos(math.toradians(iGen)))
	yGen = rGen * (math.sin(math.toradians(NGen)) * math.cos(math.toradians(nGen+wGen-NGen)) + math.cos(math.toradians(NGen)) * math.sin(math.toradians(nGen+wGen-NGen)) * math.cos(math.toradians(iGen)))
	zGen = rGen * (math.sin(math.toradians(nGen+wGen-NGen)) * math.sin(math.toradians(iGen)))

	//Check
	rCheck = math.sqrt(math.pow(xGen,2) + math.pow(yGen,2) + math.pow(zGen,2))

	[xGen, yGen, zGen, rGen]

// @function Calculate right ascension and declination for a planet relative to Earth
// @param earthX Earth X ecliptic rectangular coordinate relative to Sun
// @param earthY Earth Y ecliptic rectangular coordinate relative to Sun
// @param earthZ Earth Z ecliptic rectangular coordinate relative to Sun
// @param planetX Planet X ecliptic rectangular coordinate relative to Sun
// @param planetY Planet Y ecliptic rectangular coordinate relative to Sun
// @param planetZ Planet Z ecliptic rectangular coordinate relative to Sun
// @returns [distPlanet, raPlanetDegrees, decPlanet] Planet geocentric orbital radius, geocentric right ascension, and geocentric declination
export calculateRightAscensionAndDeclination(float earthX, float earthY, float earthZ, float planetX, float planetY, float planetZ) =>
	//CALCULATE GEOCENTRIC ECLIPTIC COORDINATES
	//=========================================
	//TAKEN FROM: http://www.stargazing.net/kepler/ellipse.html#twig06
	//CREDIT: Keith Burnett
	xGeo = planetX - earthX
	yGeo = planetY - earthY
	zGeo = planetZ - earthZ

	//CALCULATE GEOCENTRIC EQUATORIAL COORDINATES
	//===========================================
	//TAKEN FROM: http://www.stargazing.net/kepler/ellipse.html#twig06
	//CREDIT: Keith Burnett
	meanObl = meanObliquityForJ2000()
	xEq = xGeo
	//Yq = Y' * cos(ec) - Z' * sin(ec)
	//NOTE: Hardcoded the earth's obliquity (i.e. axial tilt; angle between orbital plane & equatorial plane) based on what it was on January 1, 2000 (J2000)
	yEq = (math.cos(math.toradians(meanObl)) * yGeo) - (math.sin(math.toradians(meanObl)) * zGeo )
	//Zq = Y' * sin(ec) + Z' * cos(ec)
	//NOTE: Hardcoded the earth's obliquity (i.e. axial tilt; angle between orbital plane & equatorial plane) based on what it was on January 1, 2000 (J2000)
	zEq = (math.sin(math.toradians(meanObl)) * yGeo) + (math.cos(math.toradians(meanObl)) * zGeo )

	//CALCULATE RIGHT ASCENSION (R.A.)
	//================================
	//TAKEN FROM: http://www.stargazing.net/kepler/ellipse.html#twig06
	//CREDIT: Keith Burnett
	//
	//Calculate Right Ascension (R.A) --> Horizontal location in sky
	//360° = 24 hours (15° = 1 hour)
	//
	//alpha = math.atan(Yq/Xq)
	//If Xq is negative then add 180 degrees to alpha
	//If Xq is positive and Yq is negative then add 360 degrees to alpha
	float raPlanet = math.todegrees(math.atan(yEq/xEq))
	if xEq < 0
	    raPlanet := raPlanet + 180
	else if xEq > 0 and yEq < 0
	    raPlanet := raPlanet + 360
	raPlanetDegrees = raPlanet
	//Convert right ascension to hours
	raPlanet := raPlanet / 15.0

	//CALCULATE DECLINATION
	//=====================
	//TAKEN FROM: http://www.stargazing.net/kepler/ellipse.html#twig06
	//CREDIT: Keith Burnett
	//
	//Calculate Declination --> Vertical location in sky
	//delta = arctan( Zq / SQRT(Xq^2 + Yq^2))
	decPlanet = math.todegrees( math.atan((zEq / (math.sqrt(math.pow(xEq,2) + math.pow(yEq,2))))) )

	//CALCULATE DISTANCE
	//==================
	//Calculate Distance
	distPlanet = math.sqrt( math.pow(xEq,2) + math.pow(yEq,2) + math.pow(zEq,2) )

	//CONVERT R.A. FROM DEGREES TO HOURS, MINUTES & SECONDS
	//=====================================================
	//Degrees to hours
	raPlanetH = math.floor(raPlanet)
	//Degrees to minutes
	minfloat = (raPlanet - raPlanetH) * 60
	raPlanetM = math.floor(minfloat)
	//Degrees to seconds
	secfloat = (minfloat - raPlanetM) * 60.0
	raPlanetS = math.round(secfloat)

	//CONVERT DECLINATION FROM DEGREES TO DEGREES, MINUTES & SECONDS
	//==============================================================
	//Make value positive for correct calculation
	negative = false
	if decPlanet < 0
	    decPlanet := math.abs(decPlanet)
	    negative := true

	//Degrees to hours
	int decPlanetH = math.floor(decPlanet)
	//Degrees to minutes
	minfloat := (decPlanet - decPlanetH) * 60
	decPlanetM = math.floor(minfloat)
	//Degrees to seconds
	secfloat := (minfloat - decPlanetM) * 60
	decPlanetS = math.round(secfloat)

	//Correct degrees back to negative/positive
	if negative
		decPlanetH := -1 * decPlanetH
		decPlanet := -1 * decPlanet

	[distPlanet, raPlanetDegrees, decPlanet]

// ========================================================================================================================
//                                              HELIO EPHEMERIS
// ========================================================================================================================

// @function Compute Mercury heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Mercury heliocentric longitude on date
export mercuryHelio(float T) =>
	[xMercury, yMercury, zMercury, rMercury] = planetEphemerisCalc(T, mercuryElements(), mercuryRates())
	float mercuryHelio = na
	if yMercury != 0 and xMercury != 0
		mercuryHelio := rev360(math.todegrees(atan2(yMercury, xMercury)))
	mercuryHelio

// @function Compute Venus heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Venus heliocentric longitude on date
export venusHelio(float T) =>
	[xVenus, yVenus, zVenus, rVenus] = planetEphemerisCalc(T, venusElements(), venusRates())
	float venusHelio = na
	if yVenus != 0 and xVenus != 0
		venusHelio := rev360(math.todegrees(atan2(yVenus, xVenus)))
	venusHelio

// @function Compute Earth heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Earth heliocentric longitude on date
export earthHelio(float T) =>
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	float earthHelio = na
	if yEarth != 0 and xEarth != 0
		earthHelio := rev360(math.todegrees(atan2(yEarth, xEarth)))
	earthHelio

// @function Compute Mars heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Mars heliocentric longitude on date
export marsHelio(float T) =>
	[xMars, yMars, zMars, rMars] = planetEphemerisCalc(T, marsElements(), marsRates())
	float marsHelio = na
	if yMars != 0 and xMars != 0
		marsHelio := rev360(math.todegrees(atan2(yMars, xMars)))
	marsHelio

// @function Compute Jupiter heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Jupiter heliocentric longitude on date
export jupiterHelio(float T) =>
	[xJupiter, yJupiter, zJupiter, rJupiter] = planetEphemerisCalc(T, jupiterElements(), jupiterRates())
	float jupiterHelio = na
	if yJupiter != 0 and xJupiter != 0
		jupiterHelio := rev360(math.todegrees(atan2(yJupiter, xJupiter)))
	jupiterHelio

// @function Compute Saturn heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Saturn heliocentric longitude on date
export saturnHelio(float T) =>
	[xSaturn, ySaturn, zSaturn, rSaturn] = planetEphemerisCalc(T, saturnElements(), saturnRates())
	float saturnHelio = na
	if ySaturn != 0 and xSaturn != 0
		saturnHelio := rev360(math.todegrees(atan2(ySaturn, xSaturn)))
	saturnHelio

// @function Compute Neptune heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Neptune heliocentric longitude on date
export neptuneHelio(float T) =>
	[xNeptune, yNeptune, zNeptune, rNeptune] = planetEphemerisCalc(T, neptuneElements(), neptuneRates())
	float neptuneHelio = na
	if yNeptune != 0 and xNeptune != 0
		neptuneHelio := rev360(math.todegrees(atan2(yNeptune, xNeptune)))
	neptuneHelio

// @function Compute Uranus heliocentric longitude on date
// @param Julian Centuries since epoch J2000 (2000-01-01)
// @returns Uranus heliocentric longitude on date
export uranusHelio(float T) =>
	[xUranus, yUranus, zUranus, rUranus] = planetEphemerisCalc(T, uranusElements(), uranusRates())
	float uranusHelio = na
	if yUranus != 0 and xUranus != 0
		uranusHelio := rev360(math.todegrees(atan2(yUranus, xUranus)))
	uranusHelio

// ========================================================================================================================
//                                              GEO EPHEMERIS
// ========================================================================================================================

export sunGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Earth Celestial Coordinates
	[distEarth, earthGeo, decEarth] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xEarth, yEarth, zEarth)
	// Sun Geocentric Longitude
	sunGeo = rev360(earthGeo + 180)
	sunGeo

export mercuryGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Mercury Heliocentric Ecliptic Coordinates
	[xMercury, yMercury, zMercury, rMercury] = planetEphemerisCalc(T, mercuryElements(), mercuryRates())
	// Mercury Celestial Coordinates
	[distMercury, mercuryGeo, decMercury] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xMercury, yMercury, zMercury)
	mercuryGeo

export venusGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Venus Heliocentric Ecliptic Coordinates
	[xVenus, yVenus, zVenus, rVenus] = planetEphemerisCalc(T, venusElements(), venusRates())
	// Venus Celestial Coordinates
	[distVenus, venusGeo, decVenus] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xVenus, yVenus, zVenus)
	venusGeo

export marsGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Mars Heliocentric Ecliptic Coordinates
	[xMars, yMars, zMars, rMars] = planetEphemerisCalc(T, marsElements(), marsRates())
	// Mars Celestial Coordinates
	[distMars, marsGeo, decMars] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xMars, yMars, zMars)
	marsGeo

export jupiterGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Jupiter Heliocentric Ecliptic Coordinates
	[xJupiter, yJupiter, zJupiter, rJupiter] = planetEphemerisCalc(T, jupiterElements(), jupiterRates())
	// Jupiter Celestial Coordinates
	[distJupiter, jupiterGeo, decJupiter] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xJupiter, yJupiter, zJupiter)
	jupiterGeo

export saturnGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Saturn Heliocentric Ecliptic Coordinates
	[xSaturn, ySaturn, zSaturn, rSaturn] = planetEphemerisCalc(T, saturnElements(), saturnRates())
	// Saturn Celestial Coordinates
	[distSaturn, saturnGeo, decSaturn] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xSaturn, ySaturn, zSaturn)
	saturnGeo

export neptuneGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Neptune Heliocentric Ecliptic Coordinates
	[xNeptune, yNeptune, zNeptune, rNeptune] = planetEphemerisCalc(T, neptuneElements(), neptuneRates())
	// Neptune Celestial Coordinates
	[distNeptune, neptuneGeo, decNeptune] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xNeptune, yNeptune, zNeptune)
	neptuneGeo

export uranusGeo(float T) =>
	// Get Earth Heliocentric Ecliptic Coordinates
	// Need to calculate these values first, as you need them for all following calculations
	[xEarth, yEarth, zEarth, rEarth] = planetEphemerisCalc(T, earthElements(), earthRates())
	// Uranus Heliocentric Ecliptic Coordinates
	[xUranus, yUranus, zUranus, rUranus] = planetEphemerisCalc(T, uranusElements(), uranusRates())
	// Uranus Celestial Coordinates
	[distUranus, uranusGeo, decUranus] = calculateRightAscensionAndDeclination(xEarth, yEarth, zEarth, xUranus, yUranus, zUranus)
	uranusGeo

// ========================================================================================================================
//                                              MOON EPHEMERIS
// ========================================================================================================================

moonSumL() =>
	moonSumL = matrix.new<float>()
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,1,0,6288774)) //1
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,-1,0,1274027)) //2
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,0,0,658314)) //3
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,2,0,213618)) //4
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,1,0,0,-185116)) //5
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,0,2,-114332)) //6
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,-2,0,58793)) //7
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,-1,0,57066)) //8
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,1,0,53322)) //9
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,0,0,45758))//10
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,1,-1,0,-40923))//11
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,0,0,0,-34720))//12
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,1,1,0,-30383))//13
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,0,-2,15327))//14
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,1,2,-12528))//15
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,1,-2,10980))//16
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,0,-1,0,10675))//17
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,3,0,10034))//18
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,0,-2,0,8548))//19
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,1,-1,0,-7888))//20
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,1,0,0,-6766))//21
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,0,-1,0,-5163))//22
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,1,0,0,4987))//23
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,1,0,4036))//24
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,2,0,3994))//25
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,0,0,0,3861))//26
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,-3,0,3665))//27
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,1,-2,0,-2689))//28
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,-1,2,-2602))//29
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,-2,0,2390))//30
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,0,1,0,-2348))//31
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-2,0,0,2236))//32
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,1,2,0,-2120))//33
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,2,0,0,-2069))//34
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-2,-1,0,2048))//35
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,1,-2,-1773))//36
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,0,2,-1595))//37
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,-1,-1,0,1215))//38
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,2,2,-1110))//39
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(3,0,-1,0,-892))//40
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,1,1,0,-810))//41
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,-1,-2,0,759))//42
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,2,-1,0,-713))//43
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,2,-1,0,-700))//44
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,1,-2,0,691))//45
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,0,-2,596))//46
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,0,1,0,549))//47
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,4,0,537))//48
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,-1,0,0,520))//49
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,0,-2,0,-487))//50
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,1,0,-2,-399))//51
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,0,2,-2,-381))//52
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,1,1,0,351))//53
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(3,0,-2,0,-340))//54
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(4,0,-3,0,330))//55
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,-1,2,0,327))//56
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(0,2,1,0,-323))//57
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(1,1,-1,0,299))//58
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,3,0,294))//59
	matrix.add_row(moonSumL, matrix.rows(moonSumL), array.from(2,0,-1,-2,0))//60
	moonSumL

moonSumB() =>
	moonSumB = matrix.new<float>()
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,0,1,5128122)) //1
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,1,1,280602)) //2
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,1,-1,277693)) //3
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,0,-1,173237)) //4
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,-1,1,55413)) //5
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,-1,-1,46271)) //6
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,0,1,32573)) //7
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,2,1,17198)) //8
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,1,-1,9266)) //9
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,2,-1,8822)) //10
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,0,-1,8216)) //11
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,-2,-1,4324)) //12
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,1,1,4200)) //13
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,1,0,-1,-3359)) //14
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,-1,1,2463)) //15
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,0,1,2211)) //16
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,-1,-1,2065)) //17
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,-1,-1,-1870)) //18
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,-1,-1,1828)) //19
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,0,1,-1794)) //20
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,0,3,-1749)) //21
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,-1,1,-1565)) //22
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,0,0,1,-1491)) //23
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,1,1,-1475)) //24
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,1,-1,-1410)) //25
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,0,-1,-1344)) //26
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,0,0,-1,-1335)) //27
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,3,1,1107)) //28
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,0,-1,1021)) //29
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,-1,1,833)) //30
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,1,-3,777)) //31
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,-2,1,671)) //32
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,0,-3,607)) //33
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,2,-1,596)) //34
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,1,-1,491)) //35
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,-2,1,-451)) //36
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,3,-1,439)) //37
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,2,1,422)) //38
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,0,-3,-1,421)) //39
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,1,-1,1,-366)) //40
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,1,0,1,-351)) //41
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,0,1,331)) //42
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,1,1,315)) //43
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-2,0,-1,302)) //44
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,0,1,3,-283)) //45
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,1,1,-1,-229)) //46
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,1,0,-1,223)) //47
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,1,0,1,223)) //48
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,-2,-1,-220)) //49
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,1,-1,-1,-220)) //50
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,0,1,1,-185)) //51
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-1,-2,-1,181)) //52
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(0,1,2,1,-177)) //53
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,-2,-1,176)) //54
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,-1,-1,-1,166)) //55
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,0,1,-1,-164)) //56
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,0,1,-1,132)) //57
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(1,0,-1,-1,-119)) //58
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(4,-1,0,-1,115)) //59
	matrix.add_row(moonSumB, matrix.rows(moonSumB), array.from(2,-2,0,1,10)) //60
	moonSumB



moonSumNutation() =>
	sumNutation = matrix.new<float>()
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,0,0,1,-171996,-174.2,92025,8.9)) //1
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,0,2,2,-13187,-1.6,5736,-3.1)) //2
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,0,2,2,-2274,-0.2,977,-0.5)) //3
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,0,0,2,2062,0.2,-895,0.5)) //4
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,1,0,0,0,1426,-3.4,54,-0.1)) //5
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,0,0,712,0.1,-7,na)) //6
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,1,0,2,2,-517,1.2,224,-0.6)) //7
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,0,2,1,-386,-0.4,200,na)) //8
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,2,2,-301,na,129,-0.1)) //9
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,-1,0,2,2,217,-0.5,-95,0.3)) //10
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,1,0,0,-158,na,na,na)) //11
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,0,2,1,129,0.1,-70,na)) //12
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,-1,2,2,123,na,-53,na)) //13
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,0,0,0,63,na,na,na)) //14
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,0,1,63,0.1,-33,na)) //15
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,-1,2,2,-59,na,26,na)) //16
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,-1,0,1,-58,-0.1,32,na)) //17
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,2,1,-51,na,27,na)) //18
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,2,0,0,48,na,na,na)) //19
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,-2,2,1,46,na,-24,na)) //20
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,0,2,2,-38,na,16,na)) //21
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,2,2,2,-31,na,13,na)) //22
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,2,0,0,29,na,na,na)) //23
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,1,2,2,29,na,-12,na)) //24
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,0,2,0,26,na,na,na)) //25
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,0,2,0,-22,na,na,na)) //26
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,-1,2,1,21,na,-10,na)) //27
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,2,0,0,0,17,-0.1,na,na)) //28
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,-1,0,1,16,na,-8,na)) //29
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,2,0,2,2,-16,0.1,7,na)) //30
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,1,0,0,1,-15,na,9,na)) //31
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,1,0,1,-13,na,7,na)) //32
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,-1,0,0,1,-12,na,6,na)) //33
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,2,-2,0,11,na,na,na)) //34
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,-1,2,1,-10,na,5,na)) //35
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,1,2,2,-8,na,3,na)) //36
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,1,0,2,2,7,na,-3,na)) //37
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,1,1,0,0,-7,na,na,na)) //38
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,-1,0,2,2,-7,na,3,na)) //39
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,0,2,1,-7,na,3,na)) //40
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,1,0,0,6,na,na,na)) //41
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,2,2,2,6,na,-3,na)) //42
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,1,2,1,6,na,-3,na)) //43
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,-2,0,1,-6,na,3,na)) //44
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,0,0,0,1,-6,na,3,na)) //45
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,-1,1,0,0,5,na,na,na)) //46
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,-1,0,2,1,-5,na,3,na)) //47
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,0,0,1,-5,na,3,na)) //48
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,2,2,1,-5,na,3,na)) //49
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,0,2,0,1,4,na,na,na)) //50
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,1,0,2,1,4,na,na,na)) //51
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,-2,0,4,na,na,na)) //52
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-1,0,1,0,0,-4,na,na,na)) //53
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-2,1,0,0,0,-4,na,na,na)) //54
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(1,0,0,0,0,-4,na,na,na)) //55
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,1,2,0,3,na,na,na)) //56
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,-2,2,2,-3,na,na,na)) //57
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(-1,-1,1,0,0,-3,na,na,na)) //58
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,1,1,0,0,-3,na,na,na)) //59
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,-1,1,2,2,-3,na,na,na)) //60
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,-1,-1,2,2,-3,na,na,na)) //61
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(0,0,3,2,2,-3,na,na,na)) //62
	matrix.add_row(sumNutation, matrix.rows(sumNutation), array.from(2,-1,0,2,2,-3,na,na,na)) //63
	sumNutation


decimalToHMS(float decimal) =>
	hours = math.round(decimal)
	minutes = (decimal - math.floor(decimal)) * 60
	seconds = (minutes - math.floor(minutes)) * 60
	minutes := math.round(minutes)
	if minutes < 0
		minutes := math.abs(minutes)
	if seconds < 0
		seconds := math.abs(seconds)
	[hours, minutes, seconds]


export moonGeo(float T_JD) =>
	// Moon mean longitude
	L = rev360( 218.3164477 + (481267.88123421 * T_JD) - (0.0015786 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 538841) - (math.pow(T_JD,4) / 65194000) )
	// Moon mean elongation
	D = rev360( 297.8501921 + (445267.1114034 * T_JD) - (0.0018819 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 545868) - (math.pow(T_JD,4) / 113065000) )
	// Sun mean anomaly
	M_sun = rev360( 357.5291092 + (35999.0502909 * T_JD) - (0.0001536 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 24490000) )
	// Moon mean anomaly
	M_moon = rev360( 134.9633964 + (477198.8675055 * T_JD) + (0.0087414 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 69699) - (math.pow(T_JD,4) / 14712000) )
	// Moon argument of latitude
	F = rev360( 93.2720950 + (483202.0175233 * T_JD) - (0.0036539 * math.pow(T_JD,2)) - (math.pow(T_JD,3) / 3526000) + (math.pow(T_JD,4) / 863310000) )
	// Moon longitude of ascending node
	omega = rev360 ( 125.04452 - (1934.136261 * T_JD) + (0.0020708 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 450000) )

	// Eccentricity due to Venus
	A1 = rev360( 119.75 + (131.849 + T_JD) )
	// Eccentricity due to Jupiter
	A2 = rev360( 53.09 + (479264.29 * T_JD) )
	// ?
	A3 = rev360( 313.45 + (481266.484 * T_JD) )

	// Eccentricity of orbit
	E = rev360( 1 - (0.002516 * T_JD) - (0.0000074 * math.pow(T_JD,2)) )

	// RIGHT ASCENSION
	// Calculate sum of sigma l
	float sumL = 0
	moonSumL = moonSumL()
	for int i = 0 to matrix.rows(moonSumL) - 1
		// [0] = D, Moon mean elongation
		// [1] = M, Sun mean anomaly
		// [2] = M_prime, Moon mean anomaly
		// [3] = F, Moon arugment of latitude
		// [4] = Coefficient
		terms = matrix.row(moonSumL, i)
		termD = array.get(terms, 0) * D
		termMSun = array.get(terms, 1) * M_sun
		termMMoon = array.get(terms, 2) * M_moon
		termF = array.get(terms, 3) * F
		term = termD + termMSun + termMMoon + termF
		term := math.sin(math.toradians(term))
		// term M depends on eccentricity of Earth's orbit around Sun
		if array.get(terms, 1) != 0
			// in case of last term, the coefficient does not exist, so exlude it
			// multiply sum of fundamental arugs by coefficient (include E)
			if array.get(terms, 4) != 0
				term := array.get(terms, 4) * E * term
			else
				term := E * term
		else
			if array.get(terms, 4) != 0
				// multiply sum of fundamental arugs by coefficient (include E)
				term := array.get(terms, 4) * term
		sumL := sumL + term
	sumL := sumL + (3958 + math.sin(math.toradians(A1)))
	sumL := sumL + (1962 + math.sin(math.toradians(L - F)))
	sumL := sumL + (318 + math.sin(math.toradians(A2)))

	// DECLINATION
	float sumB = 0
	moonSumB = moonSumB()
	for int i = 0 to matrix.rows(moonSumB) - 1
		terms = matrix.row(moonSumB, i)
		termD = array.get(terms, 0) * D
		termMSun = array.get(terms, 1) * M_sun
		termMMoon = array.get(terms, 2) * M_moon
		termF = array.get(terms, 3) * F
		term = termD + termMSun + termMMoon + termF
		term := math.sin(math.toradians(term))
		if array.get(terms, 1) != 0
			if array.get(terms, 4) != 0
				term := array.get(terms, 4) * E * term
			else
				term := E * term
		else
			if array.get(terms, 4) != 0
				term := array.get(terms, 4) * term
		sumB := sumB + term
	sumB := sumB - (2235 * math.sin(math.toradians(L)))
	sumB := sumB + (382 * math.sin(math.toradians(A3)))
	sumB := sumB + (175 * math.sin(math.toradians(A1 - F)))
	sumB := sumB + (175 * math.sin(math.toradians(A1 + F)))
	sumB := sumB + (127 * math.sin(math.toradians(L - M_moon)))
	sumB := sumB - (115 * math.sin(math.toradians(L + M_moon)))

	moonRADeg = L + (sumL / 1000000)
	moonDecDeg = sumB / 1000000

	// Mean longitude of ascending node of Moon's mean orbit on the ecliptic
	MLAN = rev360( 125.04452 - (1934.136261 * T_JD) + (0.0020709 * math.pow(T_JD,2)) + (math.pow(T_JD,3) / 450000) )

	float sumNutation = 0
	float sumObl = 0
	moonSumNutation = moonSumNutation()
	for int i = 0 to matrix.rows(moonSumNutation) - 1
		//[0] = D (Moon Mean Elongation)
		//[1] = M (Sun Mean Anomaly)
		//[2] = M' (Moon Mean Anomaly)
		//[3] = F (Moon Argument of latitude)
		//[4] = Ω (Moon Mean Longitude of the ascending node)
		//[5] = Coefficient I
		//[6] = Coefficient II
		//[7] = e coefficient I
		//[8] = e coefficient II
		terms = matrix.row(moonSumNutation, i)
		termD = array.get(terms, 0) * D
		termMSun = array.get(terms, 1) * M_sun
		termMMoon = array.get(terms, 2) * M_moon
		termF = array.get(terms, 3) * F
		termPhi = array.get(terms, 4) * MLAN
		term = termD + termMSun + termMMoon + termF + termPhi
		termSin = math.sin(math.toradians(term))
		termCos = math.cos(math.toradians(term))
		if not na(array.get(terms, 6))
			termSin := (array.get(terms, 5) + (array.get(terms, 6) * T_JD)) * termSin
		else
			termSin := array.get(terms, 5) * termSin
		if not na(array.get(terms, 7))
			if not na(array.get(terms, 8))
				termCos := (array.get(terms, 7) + array.get(terms, 8) * T_JD) * termCos
			else
				termCos := array.get(terms, 7) * termCos
		sumNutation := sumNutation + termSin
		sumObl := sumObl + termCos
	sumNutation := sumNutation * 0.0001 / 3600
	sumObl := sumObl * 0.0001 / 3600

	appMoonRADeg = moonRADeg + sumNutation

	// Mean Obliquitiy of Ecliptic & True Obliquity of Ecliptic
	float earthObl0 = 23.43929 - (0.01300417 * T_JD) - (0.0000001638889 * math.pow(T_JD,2)) - (0.0000005036111 * math.pow(T_JD,3)) // E0
	float earthObl = earthObl0 + sumObl // true obliquity of ecliptic
	// earthObl0 := decimalToHMS(earthObl0) // convert decimal to H,M,S

	findY = (math.sin(math.toradians(appMoonRADeg)) * math.cos(math.toradians(earthObl))) - (math.tan(math.toradians(moonDecDeg)) * math.sin(math.toradians(earthObl)))
	findX = math.cos(math.toradians(appMoonRADeg))
	findAtan = math.todegrees(math.atan(findY/findX))
	// remove inverse tan ambiguity
	if findX < 0
		findAtan := findAtan + 180
	else if findX > 0 and findY < 0
		findAtan := findAtan + 360
	// calculate right ascension -> horizontal location in sky
	// convert Moon ecliptic longitude form degrees to H,M,S
	moonRADeg := findAtan
	moonRADeg


// ========================================================================================================================
//                                              COMPOSITE SINE WAVE
// ========================================================================================================================

// @function Mercury orbital period in Earth days
// @returns 87.9691
export mercuryOrbitalPeriod() => 
	87.9691
// @function Venus orbital period in Earth days
// @returns 224.701
export venusOrbitalPeriod() =>
	224.701
// @function Earth orbital period in Earth days
// @returns 365.256363004
export earthOrbitalPeriod() =>
	365.256363004
// @function Mars orbital period in Earth days
// @returns 686.980
export marsOrbitalPeriod() =>
	686.980
// @function Jupiter orbital period in Earth days
// @returns 4332.59
export jupiterOrbitalPeriod() =>
	4332.59
// @function Saturn orbital period in Earth days
// @returns 10759.22
export saturnOrbitalPeriod() =>
	10759.22
// @function Uranus orbital period in Earth days
// @returns 30688.5
export uranusOrbitalPeriod() =>
	30688.5
// @function Neptune orbital period in Earth days
// @returns 60195.0
export neptuneOrbitalPeriod() =>
	60195.0

// @returns Jupiter-Saturn composite sine wave period (period between conjunct dates) in Julian days
export jupiterSaturnCompositePeriod() =>
	getJulianDate(2020, 10, 5, 0, 0) - getJulianDate(1961, 6, 23, 0, 0)
// @returns Jupiter-Neptune composite sine wave period (period between conjunct dates) in Julian days
export jupiterNeptuneCompositePeriod() =>
	(getJulianDate(2020, 3, 23, 0, 0) - getJulianDate(1901, 8, 2, 0, 0)) / 10
// @returns Jupiter-Uranus composite sine wave period (period between conjunct dates) in Julian days
export jupiterUranusCompositePeriod() =>
	(getJulianDate(2020, 5, 11, 0, 0) - getJulianDate(1937, 4, 20, 0, 0)) / 7
// @returns Saturn-Neptune composite sine wave period (period between conjunct dates) in Julian days
export saturnNeptuneCompositePeriod() =>
	(getJulianDate(1996, 6, 7, 0, 0) - getJulianDate(1908, 1, 15, 0, 0)) / 3

// @returns Saturn-Uranus composite sine wave period (period between conjunct dates) in Julian days
export saturnUranusCompositePeriod() =>
	// Uses 1.5 period as the period because Bradley Cowen uses a 45 year cycle in his book... Why is this?
	// SU_h = getJulianDate(2008, 12, 4, 0, 0)
	// SU_g = getJulianDate(1997, 5, 6, 0, 0)
	// SU_f = getJulianDate(1981, 11, 23, 0, 0)
	// SU_e = getJulianDate(1965, 6, 28, 0, 0)
	// SU_d = getJulianDate(1952, 6, 2, 0, 0)
	// SU_c = getJulianDate(1938, 3, 14, 0, 0)
	// SU_b = getJulianDate(1920, 8, 26, 0, 0)
	// SU_a = getJulianDate(1908, 8, 31, 0, 0)
	// SUPeriod = ((SU_h-SU_e) + (SU_g-SU_d) + (SU_f-SU_c) + (SU_e-SU_b) + (SU_d-SU_a)) / 5
	SU_g = getJulianDate(1997, 5, 6, 0, 0)
	SU_e = getJulianDate(1965, 6, 28, 0, 0)
	SU_c = getJulianDate(1938, 3, 14, 0, 0)
	SU_a = getJulianDate(1908, 8, 31, 0, 0)
	SUPeriod_tops = (SU_g - SU_a) / 3

	SU_h = getJulianDate(2008, 12, 4, 0, 0)
	SU_f = getJulianDate(1981, 11, 23, 0, 0)
	SU_d = getJulianDate(1952, 6, 2, 0, 0)
	SU_b = getJulianDate(1920, 8, 26, 0, 0)
	SUPeriod_bottoms = (SU_h - SU_b) / 3

	SUPeriod = (SUPeriod_tops + SUPeriod_bottoms) / 2

// @function Convert heliocentric longitude of planet into a sine wave
// @param julianDate Julian Date in Julian centuries
// @param planetOrbitalPeriod Orbital period of planet in Earth days
// @param planetHelio Heliocentric longitude of planet in degrees
// @returns Sine of heliocentric longitude on a Julian date
export planetSineWave(float julianDateInCenturies, float planetOrbitalPeriod, float planetHelio) =>
	math.sin(( julianDateInCenturies * julianCenturyInJulianDays() / planetOrbitalPeriod ) * math.pi * 2 - math.toradians(planetHelio))