astro.py

from collections import namedtuple
import numpy as np
d2r, r2d = np.pi/180.0, 180.0/np.pi

# Most of this is based around Meeus's Astronomical Algorithms, since it
# presents reasonably good approximations of all the quantities we require in a
# clear fashion.  Reluctant to go all out and use VSOP87 unless it can be shown
# to make a significant difference to the resulting accuracy of harmonic
# analysis.

#Convert a sexagesimal angle into decimal degrees
def s2d(degrees, arcmins = 0, arcsecs = 0, mas = 0, muas = 0):
	return (
			degrees
			+ (arcmins /  60.0)
			+ (arcsecs / (60.0*60.0))
			+ (mas	   / (60.0*60.0*1e3))
			+ (muas    / (60.0*60.0*1e6))
	)

#Evaluate a polynomial at argument
def polynomial(coefficients, argument):
	return sum([c * (argument ** i) for i,c in enumerate(coefficients)])

#Evaluate the first derivative of a polynomial at argument
def d_polynomial(coefficients, argument):
	return sum([c * i * (argument ** (i-1)) for i,c in enumerate(coefficients)])

#Meeus formula 11.1
def T(t):
	return (JD(t) - 2451545.0)/36525

#Meeus formula 7.1
def JD(t):
	Y, M = t.year, t.month
	D = (
		t.day
		+ t.hour / (24.0)
		+ t.minute / (24.0*60.0)
		+ t.second / (24.0*60.0*60.0)
		+ t.microsecond / (24.0 * 60.0 * 60.0 * 1e6)
	)
	if M <= 2:
		Y = Y - 1
		M = M + 12
	A = np.floor(Y / 100.0)
	B = 2 - A + np.floor(A / 4.0)
	return np.floor(365.25*(Y+4716)) + np.floor(30.6001*(M+1)) + D + B - 1524.5

#Meeus formula 21.3
terrestrial_obliquity_coefficients = (
	s2d(23,26,21.448),
	-s2d(0,0,4680.93),
	-s2d(0,0,1.55),
	s2d(0,0,1999.25),
	-s2d(0,0,51.38),
	-s2d(0,0,249.67),
	-s2d(0,0,39.05),
	s2d(0,0,7.12),
	s2d(0,0,27.87),
	s2d(0,0,5.79),
	s2d(0,0,2.45)
)

#Adjust these coefficients for parameter T rather than U
terrestrial_obliquity_coefficients = [
	c * (1e-2) ** i for i,c in enumerate(terrestrial_obliquity_coefficients)
]

#Not entirely sure about this interpretation, but this is the difference
#between Meeus formulae 24.2 and 24.3 and seems to work
solar_perigee_coefficients = (
	280.46645 - 357.52910,
	36000.76932 - 35999.05030,
	0.0003032 + 0.0001559,
	0.00000048
)

#Meeus formula 24.2
solar_longitude_coefficients = (
	280.46645,
	36000.76983,
	0.0003032
)

#This value is taken from JPL Horizon and is essentially constant
lunar_inclination_coefficients = (
	5.145,
)

#Meeus formula 45.1
lunar_longitude_coefficients = (
	218.3164591,
	481267.88134236,
	-0.0013268,
	1/538841.0
	-1/65194000.0
)

#Meeus formula 45.7
lunar_node_coefficients = (
	125.0445550,
	-1934.1361849,
	0.0020762,
	1/467410.0,
	-1/60616000.0
)

#Meeus, unnumbered formula directly preceded by 45.7
lunar_perigee_coefficients = (
	83.3532430,
	4069.0137111,
	-0.0103238,
	-1/80053.0,
	1/18999000.0
)

#Now follow some useful auxiliary values, we won't need their speed.
#See notes on Table 6 in Schureman for I, nu, xi, nu', 2nu''
def _I(N, i, omega):
	N, i, omega = d2r * N, d2r*i, d2r*omega
	cosI = np.cos(i)*np.cos(omega)-np.sin(i)*np.sin(omega)*np.cos(N)
	return r2d*np.arccos(cosI)

def _xi(N, i, omega):
	N, i, omega = d2r * N, d2r*i, d2r*omega
	e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N)
	e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N)
	e1, e2 = np.arctan(e1), np.arctan(e2)
	e1, e2 = e1 - 0.5*N, e2 - 0.5*N
	return -(e1 + e2)*r2d

def _nu(N, i, omega):
	N, i, omega = d2r * N, d2r*i, d2r*omega
	e1 = np.cos(0.5*(omega-i))/np.cos(0.5*(omega+i)) * np.tan(0.5*N)
	e2 = np.sin(0.5*(omega-i))/np.sin(0.5*(omega+i)) * np.tan(0.5*N)
	e1, e2 = np.arctan(e1), np.arctan(e2)
	e1, e2 = e1 - 0.5*N, e2 - 0.5*N
	return (e1 - e2)*r2d

#Schureman equation 224
#Can we be more precise than B "the solar coefficient" = 0.1681?
def _nup(N, i, omega):
	I = d2r * _I(N, i, omega)
	nu = d2r * _nu(N, i, omega)
	return r2d * np.arctan(np.sin(2*I)*np.sin(nu)/(np.sin(2*I)*np.cos(nu)+0.3347))

#Schureman equation 232
def _nupp(N, i, omega):
	I = d2r * _I(N, i, omega)
	nu = d2r * _nu(N, i, omega)
	tan2nupp = (np.sin(I)**2*np.sin(2*nu))/(np.sin(I)**2*np.cos(2*nu)+0.0727)
	return r2d * 0.5 * np.arctan(tan2nupp)

AstronomicalParameter = namedtuple('AstronomicalParameter', ['value', 'speed'])

def astro(t):
	a = {}
	#We can use polynomial fits from Meeus to obtain good approximations to
	#some astronomical values (and therefore speeds).
	polynomials = {
			's':     lunar_longitude_coefficients,
			'h':     solar_longitude_coefficients,
			'p':     lunar_perigee_coefficients,
			'N':     lunar_node_coefficients,
			'pp':    solar_perigee_coefficients,
			'90':    (90.0,),
			'omega': terrestrial_obliquity_coefficients,
			'i':     lunar_inclination_coefficients
	}
	#Polynomials are in T, that is Julian Centuries; we want our speeds to be
	#in the more convenient unit of degrees per hour.
	dT_dHour = 1 / (24 * 365.25 * 100)
	for name, coefficients in polynomials.items():
		a[name] = AstronomicalParameter(
				np.mod(polynomial(coefficients, T(t)), 360.0),
				d_polynomial(coefficients, T(t)) * dT_dHour
		)

	#Some other parameters defined by Schureman which are dependent on the
	#parameters N, i, omega for use in node factor calculations. We don't need
	#their speeds.
	args = list(each.value for each in [a['N'], a['i'], a['omega']])
	for name, function in {
		'I':    _I,
		'xi':   _xi,
		'nu':   _nu,
		'nup':  _nup,
		'nupp': _nupp
	}.items():
		a[name] = AstronomicalParameter(np.mod(function(*args), 360.0), None)

	#We don't work directly with the T (hours) parameter, instead our spanning
	#set for equilibrium arguments #is given by T+h-s, s, h, p, N, pp, 90.
	#This is in line with convention.
	hour = AstronomicalParameter((JD(t) - np.floor(JD(t))) * 360.0, 15.0)
	a['T+h-s'] = AstronomicalParameter(
		hour.value + a['h'].value - a['s'].value,
		hour.speed + a['h'].speed - a['s'].speed
	)
	#It is convenient to calculate Schureman's P here since several node
	#factors need it, although it could be argued that these
	#(along with I, xi, nu etc) belong somewhere else.
	a['P'] = AstronomicalParameter(
		np.mod(a['p'].value -a['xi'].value,360.0),
		None
	)
	return a