yuxing/lunar-calendar
1
0
Fork 0
mirror of https://github.com/infinet/lunar-calendar.git synced 2026-01-12 21:17:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

remove duplicated code from aa_full

Browse Source
This commit is contained in:
Chen Wei
2014-03-12 20:50:48 +08:00
parent 992686cc0c
commit 4da7b14e76
3 changed files with 48 additions and 799 deletions
Download Patch File Download Diff File Expand all files Collapse all files

216
aa.py
View File

@@ -1481,211 +1481,6 @@ CV = M_PHASE * DEG2RAD
C_V, CT_V, CTT_V = hsplit(CV, 3)
A_V, AT_V, ATT_V = hsplit(M_AMP, 3)
# for VSOP2010
# Mean Longituee J2000 (radian)
ci0 = [
4.402608631669e0, # Mercury
3.176134461576e0, # Venus
1.753470369433e0, # Earth-Moon Barycenter
6.203500014141e0, # Mars
4.091362210690e0, # Vesta
1.713743790353e0, # Iris
5.598651923117e0, # Bamberga
2.805135511956e0, # Ceres
2.326992146758e0, # Pallas
0.599546107035e0, # Jupiter
0.874018510107e0, # Saturn
5.548122539566e0, # Uranus
5.311897933164e0, # Neptune
0.e0, #
5.19846640063e0, # Moon (e)
1.62790513602e0, # Moon (F)
2.35555563875e0 # Moon (l)
]
# Mean Motions in longituee (radian/cy)
ci1 = [
0.2608790314068555e5, # Mercury
0.1021328554743445e5, # Venus
0.6283075850353215e4, # Earth-Moon Barycenter
0.3340612434145457e4, # Mars
0.1731170540074402e4, # Vesta
0.1704450784022772e4, # Iris
0.1428949097282629e4, # Bamberga
0.1364756486739947e4, # Ceres
0.1361923496417814e4, # Pallas
0.5296909615623250e3, # Jupiter
0.2132990086108488e3, # Saturn
0.7478165090307780e2, # Uranus
0.3813292073732270e2, # Neptune
0.3595362366859080e0, # Pluto (Mu from TOP2010)
0.777137714481804e5, # Moon (e)
0.843346615717837e5, # Moon (F)
0.832869142477147e5 # Moon (l)
]
# Planetary frequency in longituee
freqpla = [
0.2608790314068555e5, # Mecrcury
0.1021328554743445e5, # Venus
0.6283075850353215e4, # Earth-Moon Barycenter
0.3340612434145457e4, # Mars
0.5296909615623250e3, # Jupiter
0.2132990861084880e3, # Saturn
0.7478165903077800e2, # Uranus
0.3813292737322700e2, # Neptune
0.2533634111740826e2 # Pluto
]
def sunbyvsop2010(vsopobj, jde, ignorenutation=False):
''' calculate the apprent place of the Sun.
Arg:
jde as jde
Return:
geocentric longitude in radians, 0 - 2pi
'''
heliolong = vsopobj.lon(jde)
geolong = heliolong + PI
#if FK5:
#B = vsopLx(earth_B0, t)
#Lp = L - 13.97 * t - 0.031 * t * t
#deltaL = (-0.09033 + 0.03916 *
# (cos(math.radians(Lp)) + sin(math.radians(Lp))) *
# tan(math.radians(B))) / 3600.0
# appears -0.09033 is good enough
#deltaL = math.radians(-0.09033 / 3600.0)
#geolong += deltaL
# compensate nutation
if not ignorenutation:
geolong += nutation(jde)
# compensate light aberration, A&A p152, constant of aberration
# k = 20.49552"
t = (jde - J2000) / 36525.0
t2 = t * t
t3 = t2 * t
# mean anomaly of the Sun
M = 357.52910 + 35999.0503 * t - 0.0001559 * t2 - 0.00000048 * t3
# the eccentricity of Earth's orbit
e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2
# Sun's equation of center
C = ((1.9146 - 0.004817 * t - 0.000014 * t2) * sin(math.radians(M))
+ (0.019993 - 0.000101 * t) * sin(math.radians(2 * M))
+ 0.00029 * sin(math.radians(3 * M)))
# true anomaly
v = M + C
# Sun's distance from the Earth, in AU
#R = (1.000001018 * (1 - e * e)) / (1 + e * cos(math.radians(v)))
# finally, the aberration in radians
labbr = (math.radians(20.49552 / 3600.0) * (1 + e * cos(math.radians(v))) /
-1.000001018)
#print 'labbr', fmtdeg(math.degrees(labbr))
geolong += labbr
return normrad(geolong)
class VSOP2010():
''' compute mean longitude of the Earth-Moon Barycenter based on VSOP2010.
The coordinates are referred to the inertial frame defined by the dynamical
equinox and ecliptic J2000
TODO: the mean longitude returned is 1-2 degrees away from JPL.
'''
def __init__(self):
fp = open('/home/wei/tmp/astro/vsop2010/VSOP2010p3.dat')
headfmt = '9x,3i3,i7'
datafmt = '6x,4i3,1x,5i3,1x,4i4,1x,i6,1x,3i3,2a24'
term_count = 0
# iphi: 17 numerical coefficients a(i) (i=1,17)) (integer)
# ss: coefficient S
# cc: coefficient C
iphi = []
ss = []
cc = []
self.it_count = []
islon = False
delnext = False
for line in fp:
if term_count == 0:
# ip : planete index
# iv : variable index. 2 = mean longitude (radian)
# it : time power
# nt : number of terms in series
ip, iv, it, nt = fortran_readline(line, headfmt)
if iv == 2:
self.it_count.append(nt)
if it == 1:
# remove 1st term of 1st order time power
delnext = True
islon = True
else:
islon = False
term_count = nt
continue
if islon and not delnext:
tmp = fortran_readline(line, datafmt)
iphi.append(tmp[:17])
sa = tmp[17][:20] + 'e' + tmp[17][21:]
ca = tmp[18][:20] + 'e' + tmp[18][21:]
ss.append(float(sa))
cc.append(float(ca))
elif islon and delnext:
delnext = False
term_count -= 1
continue
term_count -= 1
self.it_count[1] -= 1
iphi = array(iphi)
a_ci0 = array(ci0)
a_ci1 = array(ci1)
self.aa = sum(iphi * a_ci0, axis=1)
self.bb = sum(iphi * a_ci1, axis=1)
self.ss = array(ss)
self.cc = array(cc)
def lon(self, jd, ignorenutation=False):
tn = 1
tc = (jd - J2000) / 365250.0
t = []
for it, count in enumerate(self.it_count):
tmp = [tn] * count
t.extend(tmp)
tn *= tc
t = array(t)
aa = self.aa
bb = self.bb
ss = self.ss
cc = self.cc
arg = aa + bb * tc
v = ne.evaluate('sum(t * (ss * sin(arg) + cc * cos(arg)))')
v += freqpla[2] * tc
#prec = precession(jd)
#v += prec
v = normrad(v)
#print 'debug: %f %f ' % (t[1], v)
return v
def vsopLx(vsopterms, t):
''' helper function for calculate VSOP87 '''
@@ -2008,17 +1803,18 @@ def lightabbr_high(jd):
t3 = t * t2
# mean anomaly of the Sun
M = 357.52910 + 35999.0503 * t - 0.0001559 * t2 - 0.00000048 * t3
M = (357.52910 + 35999.0503 * t - 0.0001559 * t2
- 0.00000048 * t3) * DEG2RAD
# the eccentricity of Earth's orbit
e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2
# Sun's equation of center
C = ((1.9146 - 0.004817 * t - 0.000014 * t2) * sin(math.radians(M))
+ (0.019993 - 0.000101 * t) * sin(math.radians(2 * M))
+ 0.00029 * sin(math.radians(3 * M)))
C = ((1.9146 - 0.004817 * t - 0.000014 * t2) * sin(M)
+ (0.019993 - 0.000101 * t) * sin(2 * M)
+ 0.00029 * sin(3 * M)) * DEG2RAD
# true anomaly
v = M + C
# Sun's distance from the Earth, in AU
R = (1.000001018 * (1 - e * e)) / (1 + e * cos(math.radians(v)))
R = (1.000001018 * (1 - e * e)) / (1 + e * cos(v))
res = -0.005775518 * R * var_lon * ASEC2RAD

595
aa_full.py
View File

@@ -22,111 +22,21 @@ __version__ = '0.0.3'
from numpy import *
import numexpr as ne
from aa_full_table import *
from aa import lightabbr_high
from aa import nutation
from aa import fmtdeg
from aa import g2jd, jd2g
J2000 = 2451545.0
SYNODIC_MONTH = 29.53
MOON_SPEED = 2 * pi / SYNODIC_MONTH
TROPICAL_YEAR = 365.24
SUN_SPEED = 2 * pi / TROPICAL_YEAR
ASEC2RAD = 4.848136811095359935899141e-6
DEG2RAD = 0.017453292519943295
ASEC360 = 1296000.0
PI = pi
TWOPI = 2 * pi
# post import process of LEA-406 tables
# horizontal split the numpy array
F0_V, F1_V, F2_V, F3_V, F4_V = hsplit(M_ARG, 5)
CV = M_PHASE * DEG2RAD
C_V, CT_V, CTT_V = hsplit(CV, 3)
A_V, AT_V, ATT_V = hsplit(M_AMP, 3)
# 77 terms IAU2000B Luni-Solar nutation table from USNO NOVAS c source
#
# Luni-Solar argument multipliers, coefficients, unit 1e-7 arcsec
# L L' F D Om longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin)
IAU2000BNutationTable = array([
[ 0, 0, 0, 0, 1, -172064161, -174666, 33386, 92052331, 9086, 15377 ],
[ 0, 0, 2, -2, 2, -13170906, -1675, -13696, 5730336, -3015, -4587 ],
[ 0, 0, 2, 0, 2, -2276413, -234, 2796, 978459, -485, 1374 ],
[ 0, 0, 0, 0, 2, 2074554, 207, -698, -897492, 470, -291 ],
[ 0, 1, 0, 0, 0, 1475877, -3633, 11817, 73871, -184, -1924 ],
[ 0, 1, 2, -2, 2, -516821, 1226, -524, 224386, -677, -174 ],
[ 1, 0, 0, 0, 0, 711159, 73, -872, -6750, 0, 358 ],
[ 0, 0, 2, 0, 1, -387298, -367, 380, 200728, 18, 318 ],
[ 1, 0, 2, 0, 2, -301461, -36, 816, 129025, -63, 367 ],
[ 0, -1, 2, -2, 2, 215829, -494, 111, -95929, 299, 132 ],
[ 0, 0, 2, -2, 1, 128227, 137, 181, -68982, -9, 39 ],
[ -1, 0, 2, 0, 2, 123457, 11, 19, -53311, 32, -4 ],
[ -1, 0, 0, 2, 0, 156994, 10, -168, -1235, 0, 82 ],
[ 1, 0, 0, 0, 1, 63110, 63, 27, -33228, 0, -9 ],
[ -1, 0, 0, 0, 1, -57976, -63, -189, 31429, 0, -75 ],
[ -1, 0, 2, 2, 2, -59641, -11, 149, 25543, -11, 66 ],
[ 1, 0, 2, 0, 1, -51613, -42, 129, 26366, 0, 78 ],
[ -2, 0, 2, 0, 1, 45893, 50, 31, -24236, -10, 20 ],
[ 0, 0, 0, 2, 0, 63384, 11, -150, -1220, 0, 29 ],
[ 0, 0, 2, 2, 2, -38571, -1, 158, 16452, -11, 68 ],
[ 0, -2, 2, -2, 2, 32481, 0, 0, -13870, 0, 0 ],
[ -2, 0, 0, 2, 0, -47722, 0, -18, 477, 0, -25 ],
[ 2, 0, 2, 0, 2, -31046, -1, 131, 13238, -11, 59 ],
[ 1, 0, 2, -2, 2, 28593, 0, -1, -12338, 10, -3 ],
[ -1, 0, 2, 0, 1, 20441, 21, 10, -10758, 0, -3 ],
[ 2, 0, 0, 0, 0, 29243, 0, -74, -609, 0, 13 ],
[ 0, 0, 2, 0, 0, 25887, 0, -66, -550, 0, 11 ],
[ 0, 1, 0, 0, 1, -14053, -25, 79, 8551, -2, -45 ],
[ -1, 0, 0, 2, 1, 15164, 10, 11, -8001, 0, -1 ],
[ 0, 2, 2, -2, 2, -15794, 72, -16, 6850, -42, -5 ],
[ 0, 0, -2, 2, 0, 21783, 0, 13, -167, 0, 13 ],
[ 1, 0, 0, -2, 1, -12873, -10, -37, 6953, 0, -14 ],
[ 0, -1, 0, 0, 1, -12654, 11, 63, 6415, 0, 26 ],
[ -1, 0, 2, 2, 1, -10204, 0, 25, 5222, 0, 15 ],
[ 0, 2, 0, 0, 0, 16707, -85, -10, 168, -1, 10 ],
[ 1, 0, 2, 2, 2, -7691, 0, 44, 3268, 0, 19 ],
[ -2, 0, 2, 0, 0, -11024, 0, -14, 104, 0, 2 ],
[ 0, 1, 2, 0, 2, 7566, -21, -11, -3250, 0, -5 ],
[ 0, 0, 2, 2, 1, -6637, -11, 25, 3353, 0, 14 ],
[ 0, -1, 2, 0, 2, -7141, 21, 8, 3070, 0, 4 ],
[ 0, 0, 0, 2, 1, -6302, -11, 2, 3272, 0, 4 ],
[ 1, 0, 2, -2, 1, 5800, 10, 2, -3045, 0, -1 ],
[ 2, 0, 2, -2, 2, 6443, 0, -7, -2768, 0, -4 ],
[ -2, 0, 0, 2, 1, -5774, -11, -15, 3041, 0, -5 ],
[ 2, 0, 2, 0, 1, -5350, 0, 21, 2695, 0, 12 ],
[ 0, -1, 2, -2, 1, -4752, -11, -3, 2719, 0, -3 ],
[ 0, 0, 0, -2, 1, -4940, -11, -21, 2720, 0, -9 ],
[ -1, -1, 0, 2, 0, 7350, 0, -8, -51, 0, 4 ],
[ 2, 0, 0, -2, 1, 4065, 0, 6, -2206, 0, 1 ],
[ 1, 0, 0, 2, 0, 6579, 0, -24, -199, 0, 2 ],
[ 0, 1, 2, -2, 1, 3579, 0, 5, -1900, 0, 1 ],
[ 1, -1, 0, 0, 0, 4725, 0, -6, -41, 0, 3 ],
[ -2, 0, 2, 0, 2, -3075, 0, -2, 1313, 0, -1 ],
[ 3, 0, 2, 0, 2, -2904, 0, 15, 1233, 0, 7 ],
[ 0, -1, 0, 2, 0, 4348, 0, -10, -81, 0, 2 ],
[ 1, -1, 2, 0, 2, -2878, 0, 8, 1232, 0, 4 ],
[ 0, 0, 0, 1, 0, -4230, 0, 5, -20, 0, -2 ],
[ -1, -1, 2, 2, 2, -2819, 0, 7, 1207, 0, 3 ],
[ -1, 0, 2, 0, 0, -4056, 0, 5, 40, 0, -2 ],
[ 0, -1, 2, 2, 2, -2647, 0, 11, 1129, 0, 5 ],
[ -2, 0, 0, 0, 1, -2294, 0, -10, 1266, 0, -4 ],
[ 1, 1, 2, 0, 2, 2481, 0, -7, -1062, 0, -3 ],
[ 2, 0, 0, 0, 1, 2179, 0, -2, -1129, 0, -2 ],
[ -1, 1, 0, 1, 0, 3276, 0, 1, -9, 0, 0 ],
[ 1, 1, 0, 0, 0, -3389, 0, 5, 35, 0, -2 ],
[ 1, 0, 2, 0, 0, 3339, 0, -13, -107, 0, 1 ],
[ -1, 0, 2, -2, 1, -1987, 0, -6, 1073, 0, -2 ],
[ 1, 0, 0, 0, 2, -1981, 0, 0, 854, 0, 0 ],
[ -1, 0, 0, 1, 0, 4026, 0, -353, -553, 0, -139 ],
[ 0, 0, 2, 1, 2, 1660, 0, -5, -710, 0, -2 ],
[ -1, 0, 2, 4, 2, -1521, 0, 9, 647, 0, 4 ],
[ -1, 1, 0, 1, 1, 1314, 0, 0, -700, 0, 0 ],
[ 0, -2, 2, -2, 1, -1283, 0, 0, 672, 0, 0 ],
[ 1, 0, 2, 2, 1, -1331, 0, 8, 663, 0, 4 ],
[ -2, 0, 2, 2, 2, 1383, 0, -2, -594, 0, -2 ],
[ -1, 0, 0, 0, 2, 1405, 0, 4, -610, 0, 2 ],
[ 1, 1, 2, -2, 2, 1290, 0, 0, -556, 0, 0 ]
])
def vsopLx(vsopterms, t):
''' helper function for calculate VSOP87 '''
@@ -135,6 +45,9 @@ def vsopLx(vsopterms, t):
return sum(lx)
# full VSOP87D tables
from aa_full_table import EAR_L0, EAR_L1, EAR_L2, EAR_L3, EAR_L4, EAR_L5
def vsop(jde, FK5=True):
''' Calculate ecliptical longitude of earth in heliocentric coordinates,
@@ -156,12 +69,12 @@ def vsop(jde, FK5=True):
'''
t = (jde - J2000) / 365250.0
L0 = vsopLx(earth_L0, t)
L1 = vsopLx(earth_L1, t)
L2 = vsopLx(earth_L2, t)
L3 = vsopLx(earth_L3, t)
L4 = vsopLx(earth_L4, t)
L5 = vsopLx(earth_L5, t)
L0 = vsopLx(EAR_L0, t)
L1 = vsopLx(EAR_L1, t)
L2 = vsopLx(EAR_L2, t)
L3 = vsopLx(EAR_L3, t)
L4 = vsopLx(EAR_L4, t)
L5 = vsopLx(EAR_L5, t)
lon = (L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5)))))
@@ -217,22 +130,6 @@ def npitopi(r):
return r
def fmtdeg(fdegree):
'''convert decimal degree to d m s format string'''
if abs(fdegree) > 360:
fdegree = math.fmod(fdegree, 360.0)
degree = math.fabs(fdegree)
tmp, deg = math.modf(degree)
minutes = tmp * 60
secs = math.modf(minutes)[0] * 60
sign =''
if fdegree < 0:
sign = '-'
res = '''%s%d%s%d'%.6f"''' % (sign, int(deg), u'\N{DEGREE SIGN}',
math.floor(minutes), secs)
return res
def f_solarangle(jd, r_angle):
''' Calculate the difference between target angle and solar geocentric
longitude at a given JDTT
@@ -283,7 +180,6 @@ def solarterm(year, angle):
# we searching for
est_vejd = g2jd(year, 3, 20.5)
x0 = est_vejd + angle * 360.0 / 365.24 # estimate
#x0 -= solarangle(x0, r) / SUN_SPEED # further closing
x1 = x0 + 0.5
return rootbysecand(f_solarangle, r, x0, x1, precision=ERROR)
@@ -334,9 +230,6 @@ def findnewmoons(start, count=15):
while nmcount < count:
b = newmoon(start)
if b != nm:
if nm > 0 and abs(nm - b) > (SYNODIC_MONTH + 1):
print 'last newmoon %s, found %s' % (fmtjde2ut(nm),
fmtjde2ut(b))
newmoons.append(b)
nm = b
nmcount += 1
@@ -350,7 +243,7 @@ def findnewmoons(start, count=15):
def apparentmoon(jde, ignorenutation=False):
''' calculate the apparent position of the Moon, it is an alias to the
lea406 function'''
return lea406(jde, ignorenutation)
return lea406_full(jde, ignorenutation)
def apparentsun(jde, ignorenutation=False):
@@ -374,453 +267,6 @@ def apparentsun(jde, ignorenutation=False):
return normrad(geolong)
# the higher accuracy light abberation table from A & A.
# the first item of row is the power of time, the 3rd and 4th item are the arg
# of sin, they have been converted into radians.
LIGHTABBR_TABLE = [
[ 0, 118.568, 1.527664004990, 6283.075850600876 ],
[ 0, 2.476, 1.484508993906, 12566.151699456424 ],
[ 0, 1.376, 0.486077687339, 77713.771468687730 ],
[ 0, 0.119, 1.276490181677, 7860.419392621536 ],
[ 0, 0.114, 5.885711004647, 5753.384884566103 ],
[ 0, 0.086, 3.884055717388, 11506.769769132206 ],
[ 0, 0.078, 2.841633387025, 161000.685738008644 ],
[ 0, 0.054, 1.441333038870, 18849.227550057301 ],
[ 0, 0.052, 2.993569534399, 3930.209696310768 ],
[ 0, 0.034, 0.529208263814, 71430.695618086858 ],
[ 0, 0.033, 2.091087703461, 5884.926847413107 ],
[ 0, 0.023, 4.320419446313, 5223.693920276658 ],
[ 0, 0.023, 5.674983441420, 5507.553238331428 ],
[ 0, 0.021, 2.707426294193, 11790.629088932305 ],
[ 1, 7.311, 5.819826570714, 6283.075850600876 ],
[ 1, 0.305, 5.776715192860, 12566.151699456424 ],
[ 1, 0.010, 5.733703298774, 18849.227550057301 ],
[ 2, 0.309, 4.214128894867, 6283.075850600876 ],
[ 2, 0.021, 3.578766215288, 12566.151699456424 ],
[ 2, 0.004, 5.198655163283, 77713.771468687730 ],
[ 3, 0.010, 2.700139544566, 6283.075850600876 ],
]
def lightabbr_high(jd):
'''compute light abberation based on A & A p156
the error will be less than 0.001"
'''
t = (jd - J2000) / 365250.0
# variation of the Sun's longitude
var_lon = 3548.330
for r in LIGHTABBR_TABLE:
var_lon += t ** r[0] * r[1] * sin(r[2] + r[3] * t)
t = (jd - J2000) / 36525.0
t2 = t * t
t3 = t * t2
# mean anomaly of the Sun
M = 357.52910 + 35999.0503 * t - 0.0001559 * t2 - 0.00000048 * t3
# the eccentricity of Earth's orbit
e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2
# Sun's equation of center
C = ((1.9146 - 0.004817 * t - 0.000014 * t2) * sin(math.radians(M))
+ (0.019993 - 0.000101 * t) * sin(math.radians(2 * M))
+ 0.00029 * sin(math.radians(3 * M)))
# true anomaly
v = M + C
# Sun's distance from the Earth, in AU
R = (1.000001018 * (1 - e * e)) / (1 + e * cos(math.radians(v)))
res = -0.005775518 * R * var_lon * ASEC2RAD
return res
def g2jd(y, m, d):
''' convert a Gregorian date to JD
from AA, p61
'''
if m <= 2:
y -= 1
m += 12
a = int(y / 100)
julian = False
if y < 1582:
julian = True
elif y == 1582:
if m < 10:
julian = True
if m == 10 and d <= 5:
julian = True
if m == 10 and d > 5 and d < 15:
return 2299160.5
if julian:
b = 0
else:
b = 2 - a + int(a / 4)
# 30.6001 is a hack Meeus suggested
return int(365.25 * (y + 4716)) + int(30.6001 * (m + 1)) + d + b - 1524.5
def td2jde(y, m, d):
''' convert Terrestrial (Dynamic) Time to JDE '''
return g2jd(y, m ,d)
def ut2jde(y, m, d):
''' convert UT to JDE, compensate the delta T '''
jd = g2jd(y, m, d)
deltat = deltaT(y, m) / 86400.0
return jd + deltat
def ut2jdut(y, m, d):
''' convert UT to JD, ignore delta T '''
return g2jd(y, m, d)
def jdut2ut(jd):
return jd2g(jd)
def jde2ut(jde, showtime=False):
g = jd2g(jde, showtime)
deltat = deltaT(g[0], g[1]) / 86400.0 # in day
return jd2g(jde - deltat, showtime)
def jde2td(jde):
return jd2g(jde)
def jd2g(jd):
''' convert JD to Gregorian date '''
jd += 0.5
z = int(jd)
f = fmod(jd, 1.0) # decimal part
if z < 2299161:
a = z
else:
alpha = int((z - 1867216.25) / 36524.25)
a = z + 1 + alpha - int(alpha / 4)
b = a + 1524
c = int((b - 122.1) / 365.25)
d = int(365.25 * c)
e = int((b - d) / 30.6001)
res_d = b - d - int(30.6001 * e) + f
if e < 14:
res_m = e - 1
else:
res_m = e - 13
if res_m > 2:
res_y = c - 4716
else:
res_y = c - 4715
return (res_y, res_m, res_d)
def jdptime(isodt, fmt, tz=0, ut=False):
'''
Args:
isodt: datetime string in ISO format
tz: a integer as timezone, e.g. -8 for UTC-8, 2 for UTC2
fmt: like strftime but currently only support three format
%y-%m-%d %H:%M:%S
%y-%m-%d %H:%M
%y-%m-%d
ut: convert to UTC(adjust delta T)
Return:
Julian Day
'''
if fmt == '%y-%m-%d':
isodate = isodt
isot = '00:00:00'
elif fmt == '%y-%m-%d %H:%M':
isodate, isot = isodt.split()
isot += ':00'
else:
isodate, isot = isodt.split()
y, m, d = [int(x) for x in isodate.split('-')]
#if isot:
H, M, S = [int(x) for x in isot.split(':')]
d += (H * 3600 + M * 60 + S) / 86400.0
return g2jd(y, m, d)
def jdftime(jde, fmt='%y-%m-%d %H:%M:%S', tz=0, ut=False):
''' format a Julian Day to ISO format datetime
Args:
jde: time in JDTT
tz: a integer as timezone, e.g. -8 for UTC-8, 2 for UTC2
fmt: like strftime but currently only support three format
%y-%m-%d %H:%M:%S
%y-%m-%d %H:%M
%y-%m-%d
ut: convert to UTC(adjust delta T)
Return:
time string in ISO format, e.g. 1984-01-01 23:59:59
'''
g = jd2g(jde)
deltat = deltaT(g[0], g[1]) if ut else 0
# convert JDE to seconds, then adjust deltat
utsec = jde * 86400.0 + tz * 3600 - deltat
jdut = utsec / 86400.0
# get time in seconds directly to minimize round error
secs = (utsec + 43200) % 86400
if fmt == '%y-%m-%d %H:%M':
secs = int(round(secs / 60.0) * 60.0)
else:
secs = int(secs)
if 86400 == secs:
jdut = int(jdut) + 0.5
secs = 0
y, m, d = jd2g(jdut)
d = int(d)
# use integer math hereafter
H = secs / 3600
ms = secs % 3600
M = ms / 60
S = ms % 60
# padding zero
d = '0%d' % d if d < 10 else str(d)
m = '0%d' % m if m < 10 else str(m)
H = '0%d' % H if H < 10 else str(H)
M = '0%d' % M if M < 10 else str(M)
S = '0%d' % S if S < 10 else str(S)
if fmt == '%y-%m-%d':
isodt = '%d-%s-%s' % (y, m, d)
elif fmt == '%y-%m-%d %H:%M':
isodt = '%d-%s-%s %s:%s' % (y, m, d, H, M)
else:
isodt = '%d-%s-%s %s:%s:%s' % (y, m, d, H, M, S)
return isodt
def deltaT(year, m ,d=0):
''' Polynomial Expressions for Delta T (ΔT) from nasa. Valid for -1999 to
+3000. http://eclipse.gsfc.nasa.gov/LEcat5/deltatpoly.html
Arg:
year: Gregorian year in integer
m: Gregorian month in integer
d: doesn't matter
Result:
ΔT in seconds
verfify against historical records from NASA
year history computed diff
-500 17190 17195.37 5.4
-400 15530 15523.84 6.2
-300 14080 14071.97 8.0
-200 12790 12786.59 3.4
-100 11640 11632.52 7.5
0 10580 10578.95 1.0
100 9600 9592.45 7.6
200 8640 8636.34 3.7
300 7680 7676.67 3.3
400 6700 6694.67 5.3
500 5710 5705.50 4.5
600 4740 4734.89 5.1
700 3810 3809.04 1.0
800 2960 2951.97 8.0
900 2200 2197.11 2.9
1000 1570 1571.65 1.7
1100 1090 1086.99 3.0
1200 740 735.10 4.9
1300 490 490.98 1.0
1400 320 321.09 1.1
1500 200 197.85 2.2
1600 120 119.55 0.5
1700 9 8.90 0.1
1750 13 13.44 0.4
1800 14 13.57 0.4
1850 7 7.16 0.2
1900 -3 -2.12 0.9
1950 29 29.26 0.3
1955 31 31.23 0.1
1960 33 33.31 0.1
1965 36 36.13 0.4
1970 40 40.66 0.5
1975 46 45.94 0.4
1980 50 50.93 0.4
1985 54 54.60 0.3
1990 57 57.20 0.3
1995 61 61.17 0.4
2000 64 64.00 0.2
2005 65 64.85 0.1
Also, JPL uses "last known leap-second is used over any future interval".
It causes the large error when compare apparent sun/moon position with JPL
Horizon
'''
y = year + (m - 0.5) / 12
if year < -500:
u = (year - 1820) / 100.0
return -20 + 32 * u ** 2
elif year < 500:
u = y / 100.0
return 10583.6 + u * (-1014.41 +
u * (33.78311 +
u * (-5.952053 +
u * (-0.1798452 +
u * (0.022174192 +
u * 0.0090316521)))))
elif year < 1600:
u = (y -1000) / 100.0
return 1574.2 + u * (-556.01 +
u * (71.23472 +
u * (0.319781 +
u * (-0.8503463 +
u * (-0.005050998 +
u * 0.0083572073)))))
elif year < 1700:
u = y -1600
return 120 + u * (-0.9808 +
u * (- 0.01532 +
u / 7129))
elif year < 1800:
u = y - 1700
return 8.83 + u * (0.1603 +
u * (-0.0059285 +
u * (0.00013336 +
u / -1174000)))
elif year < 1860:
u = y - 1800
return 13.72 + u * (-0.332447 +
u * (0.0068612 +
u * (0.0041116 +
u * (-0.00037436 +
u * (0.0000121272 +
u * (-0.0000001699 +
u * 0.000000000875))))))
elif year < 1900:
u = y - 1860
return 7.62 + u * (0.5737 +
u * (-0.251754 +
u * (0.01680668 +
u * (-0.0004473624 +
u / 233174))))
elif year < 1920:
u = y - 1900
return -2.79 + u * (1.494119 +
u * (-0.0598939 +
u * (0.0061966 +
u * -0.000197)))
elif year < 1941:
u = y - 1920
return 21.20 + u * (0.84493 +
u * (-0.076100 +
u * 0.0020936))
elif year < 1961:
u = y - 1950
return 29.07 + u * (0.407 +
u * (-1.0 / 233 +
u / 2547))
elif year < 1986:
u = y - 1975
return 45.45 + u * (1.067 +
u * (-1.0 / 260 +
u / -718))
elif year < 2005:
u = y -2000
return 63.86 + u * (0.3345 +
u * (-0.060374 +
u * (0.0017275 +
u * (0.000651814 +
u * 0.00002373599))))
'''
else:
# JPL uses "last known leap-second is used over any future interval" .
# It causes the large error when compare apparent sun/moon position
# with JPL Horizon
return 67.182963
'''
elif year < 2050:
u = y -2000
return 62.92 + u * (0.32217 + u * 0.005589)
elif year < 2150:
u = (y - 1820.0) / 100.0
return -20 + 32 * u ** 2 - 0.5628 * (2150 - y)
else:
u = (y - 1820.0) / 100.0
return -20 + 32 * u ** 2
def nutation(jde):
'''
Calculate nutation angles using the IAU 2000B model.
The table is from NOVAS(USNO) c source flle.
"...the 77-term, IAU-approved truncated nutation series, IAU 2000B, which
is accurate to about 0.001 arcsecond in the interval 1995-2050"
Arg:
jde as JDE
Return:
nutation of longitude in radians
'''
t = (jde - J2000) / 36525.0
#Mean anomaly of the Moon, in arcsec
L = 485868.249036 + t * 1717915923.2178
#Mean anomaly of the Sun.
Lp = 1287104.79305 + t * 129596581.0481
#Mean argument of the latitude of the Moon.
F = 335779.526232 + t * 1739527262.8478
#Mean elongation of the Moon from the Sun.
D = 1072260.70369 + t * 1602961601.2090
#Mean longitude of the ascending node of the Moon.
Om = 450160.398036 - t * 6962890.5431
m1, m2, m3, m4, m5, AA, BB, CC, EE, DD, FF = hsplit(IAU2000BNutationTable,
11)
args = (m1 * L + m2 * Lp + m3 * F + m4 * D + m5 * Om) * ASEC2RAD
lon = sum((AA + BB * t) * sin(args) + CC * cos(args))
# unit of longitude is 1.0e-7 arcsec, convert it to arcsec
lon *= 1.0e-7
# Constant account for the missing long-period planetary terms in the
# truncated nutation model, in arcsec
deplan = 0.000388
lon += deplan
lon *= ASEC2RAD # Convert from arcsec to radians
lon = fmod(lon, TWOPI)
return lon
#------------------------------------------------------------------------------
# LEA-406 Moon Solution
#
@@ -831,9 +277,16 @@ def nutation(jde):
# the tables M_AMP, M_PHASE, M_ARG are imported from aa_full_table
#------------------------------------------------------------------------------
FRM = [785939.924268, 1732564372.3047, -5.279, .006665, -5.522e-5]
from aa_full_table import M_ARG, M_AMP, M_PHASE
# post import process of LEA-406 tables, horizontal split the numpy array
F0_V, F1_V, F2_V, F3_V, F4_V = hsplit(M_ARG, 5)
CV = M_PHASE * DEG2RAD
C_V, CT_V, CTT_V = hsplit(CV, 3)
A_V, AT_V, ATT_V = hsplit(M_AMP, 3)
def lea406(jd, ignorenutation=False):
def lea406_full(jd, ignorenutation=False):
''' compute moon ecliptic longitude using lea406
numpy is used
'''
@@ -869,8 +322,8 @@ def lea406(jd, ignorenutation=False):
def main():
#jd = 2444239.5
jd = g2jd(1900, 1, 1)
for i in xrange(1):
l = normrad(lea406(jd))
for i in xrange(10):
l = normrad(lea406_full(jd))
#d = fmtdeg(math.degrees(npitopi(e -l )))
print jd, l, fmtdeg(math.degrees(l))
jd += 2000

34
aa_full_table.py
View File

@@ -31560,7 +31560,7 @@ M_PHASE = array([
# The D version of VSOP87 table, heliocentric spherical, coordinates referred
# to the mean equinox of the date
earth_L0 = array([
EAR_L0 = array([
[ 1.75347045673, 0, 0 ],
[ 0.03341656456, 4.66925680417, 6283.0758499914 ],
[ 0.00034894275, 4.62610241759, 12566.1516999828 ],
@@ -32122,7 +32122,7 @@ earth_L0 = array([
[ 0.00000000049, 2.44790934886, 13613.804277336 ],
])
earth_L1 = array([
EAR_L1 = array([
[ 6283.31966747491, 0, 0 ],
[ 0.00206058863, 2.67823455584, 6283.0758499914 ],
[ 0.0000430343, 2.63512650414, 12566.1516999828 ],
@@ -32466,7 +32466,7 @@ earth_L1 = array([
[ 0.00000000025, 5.71466092822, 25934.1243310894 ],
])
earth_L2 = array([
EAR_L2 = array([
[ 0.0005291887, 0, 0 ],
[ 0.00008719837, 1.07209665242, 6283.0758499914 ],
[ 0.00000309125, 0.86728818832, 12566.1516999828 ],
@@ -32611,7 +32611,7 @@ earth_L2 = array([
[ 0.0000000001, 4.93364992366, 12352.8526045448 ],
])
earth_L3 = array([
EAR_L3 = array([
[ 0.00000289226, 5.84384198723, 6283.0758499914 ],
[ 0.00000034955, 0, 0 ],
[ 0.00000016819, 5.48766912348, 12566.1516999828 ],
@@ -32636,7 +32636,7 @@ earth_L3 = array([
[ 0.00000000005, 4.28412873331, 6275.9623029906 ],
])
earth_L4 = array([
EAR_L4 = array([
[ 0.00000114084, 3.14159265359, 0 ],
[ 0.00000007717, 4.13446589358, 6283.0758499914 ],
[ 0.00000000765, 3.83803776214, 12566.1516999828 ],
@@ -32650,7 +32650,7 @@ earth_L4 = array([
[ 0.00000000002, 0.54912904658, 6438.4962494256 ],
])
earth_L5 = array([
EAR_L5 = array([
[ 0.00000000878, 3.14159265359, 0 ],
[ 0.00000000172, 2.7657906951, 6283.0758499914 ],
[ 0.0000000005, 2.01353298182, 155.4203994342 ],
@@ -32658,7 +32658,7 @@ earth_L5 = array([
[ 0.00000000005, 1.75600058765, 18849.2275499742 ],
])
earth_B0 = array([
EAR_B0 = array([
[ 0.0000027962, 3.19870156017, 84334.66158130829 ],
[ 0.00000101643, 5.42248619256, 5507.5532386674 ],
[ 0.00000080445, 3.88013204458, 5223.6939198022 ],
@@ -32845,7 +32845,7 @@ earth_B0 = array([
[ 0.00000000039, 3.1123991069, 96900.81328129109 ],
])
earth_B1 = array([
EAR_B1 = array([
[ 0.0000000903, 3.8972906189, 5507.5532386674 ],
[ 0.00000006177, 1.73038850355, 5223.6939198022 ],
[ 0.000000038, 5.24404145734, 2352.8661537718 ],
@@ -32947,7 +32947,7 @@ earth_B1 = array([
[ 0.00000000019, 0.85407021371, 14712.317116458 ],
])
earth_B2 = array([
EAR_B2 = array([
[ 0.00000001662, 1.62703209173, 84334.66158130829 ],
[ 0.00000000492, 2.41382223971, 1047.7473117547 ],
[ 0.00000000344, 2.24353004539, 5507.5532386674 ],
@@ -32999,7 +32999,7 @@ earth_B2 = array([
[ 0.00000000009, 5.94191743597, 7632.9432596502 ],
])
earth_B3 = array([
EAR_B3 = array([
[ 0.00000000011, 0.23877262399, 7860.4193924392 ],
[ 0.00000000009, 1.16069982609, 5507.5532386674 ],
[ 0.00000000008, 1.65357552925, 5884.9268465832 ],
@@ -33013,7 +33013,7 @@ earth_B3 = array([
[ 0.00000000007, 2.73399865247, 6309.3741697912 ],
])
earth_B4 = array([
EAR_B4 = array([
[ 0.00000000004, 0.79662198849, 6438.4962494256 ],
[ 0.00000000005, 0.84308705203, 1047.7473117547 ],
[ 0.00000000005, 0.05711572303, 84334.66158130829 ],
@@ -33021,7 +33021,7 @@ earth_B4 = array([
[ 0.00000000003, 2.89822201212, 6127.6554505572 ],
])
earth_R0 = array([
EAR_R0 = array([
[ 1.00013988799, 0, 0 ],
[ 0.01670699626, 3.09846350771, 6283.0758499914 ],
[ 0.00013956023, 3.0552460962, 12566.1516999828 ],
@@ -33550,7 +33550,7 @@ earth_R0 = array([
[ 0.0000000005, 6.15760345261, 78051.34191383338 ],
])
earth_R1 = array([
EAR_R1 = array([
[ 0.00103018608, 1.10748969588, 6283.0758499914 ],
[ 0.00001721238, 1.06442301418, 12566.1516999828 ],
[ 0.00000702215, 3.14159265359, 0 ],
@@ -33845,7 +33845,7 @@ earth_R1 = array([
[ 0.0000000002, 5.91915117116, 48739.859897083 ],
])
earth_R2 = array([
EAR_R2 = array([
[ 0.00004359385, 5.78455133738, 6283.0758499914 ],
[ 0.00000123633, 5.57934722157, 12566.1516999828 ],
[ 0.00000012341, 3.14159265359, 0 ],
@@ -33987,7 +33987,7 @@ earth_R2 = array([
[ 0.00000000009, 4.91488110218, 213.299095438 ],
])
earth_R3 = array([
EAR_R3 = array([
[ 0.00000144595, 4.27319435148, 6283.0758499914 ],
[ 0.00000006729, 3.91697608662, 12566.1516999828 ],
[ 0.00000000774, 0, 0 ],
@@ -34017,7 +34017,7 @@ earth_R3 = array([
[ 0.00000000005, 3.71102966917, 6290.1893969922 ],
])
earth_R4 = array([
EAR_R4 = array([
[ 0.00000003858, 2.56384387339, 6283.0758499914 ],
[ 0.00000000306, 2.2676950123, 12566.1516999828 ],
[ 0.00000000053, 3.44031471924, 5573.1428014331 ],
@@ -34030,7 +34030,7 @@ earth_R4 = array([
[ 0.00000000003, 1.28175749811, 6286.5989683404 ],
])
earth_R5 = array([
EAR_R5 = array([
[ 0.00000000086, 1.21579741687, 6283.0758499914 ],
[ 0.00000000012, 0.65617264033, 12566.1516999828 ],
[ 0.00000000001, 0.38068797142, 18849.2275499742 ],