' program oplanets.bas April 13, 2017

' the ephemeris for the outer planets, except for Pluto, are based
' on the algorithms described in Low-Precision Formulae for Planetary
' Positions, T. C. Van Flandern and K. F. Pulkkinen, The Astrophysical
' Journal Supplement Series, 41:391-411, November 1979. To the
' precision of these algorithms (one arc minute) these coordinates can
' be considered to be true-of-date.

' the algorithm that computes the heliocentric position vector of Pluto
' relative to the ecliptic and equinox of J2000 is based on the method
' described in Chapter 36 of Astronomical Algorithms by Jean Meeus.

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

option default float

option base 1

const dtr = pi / 180.0, rtd = 180.0 / pi, pi2 = 2.0 * pi

dim pj(43), ps(43), pp(43), plona(43), plonb(43)

dim plata(43), platb(43), prva(43), prvb(43)

dim cmonth, cday, cyear

' read data for pluto

for i% = 1 to 43

read pj(i%), ps(i%), pp(i%)

read plona(i%), plonb(i%), plata(i%), platb(i%)

read prva(i%), prvb(i%)

next i%

' data for pluto - heliocentric longitude, latitude and radius vector

data 0, 0, 1, -19798886, 19848454, -5453098, -14974876, 66867334, 68955876
data 0, 0, 2, 897499, -4955707, 3527363, 1672673, -11826086, -333765
data 0, 0, 3, 610820, 1210521, -1050939, 327763, 1593657, -1439953
data 0, 0, 4, -341639, -189719, 178691, -291925, -18948, 482443
data 0, 0, 5, 129027, -34863, 18763, 100448, -66634, -85576
data 0, 0, 6, -38215, 31061, -30594, -25838, 30841, -5765
data 0, 1, -1, 20349, -9886, 4965, 11263, -6140, 22254
data 0, 1, 0, -4045, -4904, 310, -132, 4434, 4443
data 0, 1, 1, -5885, -3238, 2036, -947, -1518, 641
data 0, 1, 2, -3812, 3011, -2, -674, -5, 792
data 0, 1, 3, -601, 3468, -329, -563, 518, 518
data 0, 2, -2, 1237, 463, -64, 39, -13, -221
data 0, 2, -1, 1086, -911, -94, 210, 837, -494
data 0, 2, 0, 595, -1229, -8, -160, -281, 616
data 1, -1, 0, 2484, -485, -177, 259, 260, -395
data 1, -1, 1, 839, -1414, 17, 234, -191, -396
data 1, 0, -3, -964, 1059, 582, -285, -3218, 370
data 1, 0, -2, -2303, -1038, -298, 692, 8019, -7869
data 1, 0, -1, 7049, 747, 157, 201, 105, 45637
data 1, 0, 0, 1179, -358, 304, 825, 8623, 8444
data 1, 0, 1, 393, -63, -124, -29, -896, -801
data 1, 0, 2, 111, -268, 15, 8, 208, -122
data 1, 0, 3, -52, -154, 7, 15, -133, 65
data 1, 0, 4, -78, -30, 2, 2, -16, 1
data 1, 1, -3, -34, -26, 4, 2, -22, 7
data 1, 1, -2, -43, 1, 3, 0, -8, 16
data 1, 1, -1, -15, 21, 1, -1, 2, 9
data 1, 1, 0, -1, 15, 0, -2, 12, 5
data 1, 1, 1, 4, 7, 1, 0, 1, -3
data 1, 1, 3, 1, 5, 1, -1, 1, 0
data 2, 0, -6, 8, 3, -2, -3, 9, 5
data 2, 0, -5, -3, 6, 1, 2, 2, -1
data 2, 0, -4, 6, -13, -8, 2, 14, 10
data 2, 0, -3, 10, 22, 10, -7, -65, 12
data 2, 0, -2, -57, -32, 0, 21, 126, -233
data 2, 0, -1, 157, -46, 8, 5, 270, 1068
data 2, 0, 0, 12, -18, 13, 16, 254, 155
data 2, 0, 1, -4, 8, -2, -3, -26, -2
data 2, 0, 2, -5, 0, 0, 0, 7, 0
data 2, 0, 3, 3, 4, 0, 1, -11, 4
data 3, 0, -2, -1, -1, 0, 1, 4, -14
data 3, 0, -1, 6, -3, 0, 0, 18, 35
data 3, 0, 0, -1, -2, 0, 1, 13, 3

cmonth = 12
cday = 13
cyear = 1986

julian(cmonth, cday, cyear, jdate)

jupiter(jdate, plat, plon, pr)

print " "
print "julian day = ", str$(jdate, 0, 8)
print " "
print "plat_jupiter = ", str$(rtd * plat, 0, 4), " degrees"
print "plon_jupiter = ", str$(rtd * plon, 0, 4), " degrees"
print "pr_jupiter = ", str$(pr, 0, 8), " au"

print " "

' example from Astronomical Algorithms, page 247-251

jdate = 2448908.5

pluto(jdate, plat, plon, pr)

print " "
print "julian day = ", str$(jdate, 0, 8)
print " "
print "plat_pluto = ", str$(rtd * plat, 0, 6), " degrees"
print "plon_pluto = ", str$(rtd * plon, 0, 6), " degrees"
print "pr_pluto = ", str$(pr, 0, 8), " au"

print " "

end

'''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''

sub jupiter (jdate, plat, plon, pr)

' heliocentric coordinates of jupiter

' input

' jdate = julian day

' output

' plat = heliocentric latitude (radians)
' plon = heliocentric longitude (radians)
' pr = heliocentric distance (au)

''''''''''''''''''''''''''''''''''''

atr = pi / 648000.0

' fundamental time arguments

djd = jdate - 2451545.0

t = djd / 36525.0 + 1.0

l5 = r2r(0.089608 + 0.00023080893 * djd)
g5 = r2r(0.056531 + 0.00023080893 * djd)
g6 = r2r(0.882987 + 0.00009294371 * djd)
g7 = r2r(0.400589 + 0.00003269438 * djd)

pl = 19934 * sin(g5) + 5023 * t + 2511 + 1093 * cos(2 * g5 - 5 * g6)
pl = pl + 601 * sin(2 * g5) - 479 * sin(2 * g5 - 5 * g6)
pl = pl - 185 * sin(2 * g5 - 2 * g6) + 137 * sin(3 * g5 - 5 * g6)
pl = pl - 131 * sin(g5 - 2 * g6) + 79 * cos(g5 - g6)
pl = pl - 76 * cos(2 * g5 - 2 * g6) - 74 * t * cos(g5)
pl = pl + 68 * t * sin(g5) + 66 * cos(2 * g5 - 3 * g6)
pl = pl + 63 * cos(3 * g5 - 5 * g6) + 53 * cos(g5 - 5 * g6)
pl = pl + 49 * sin(2 * g5 - 3 * g6) - 43 * t * sin(2 * g5 - 5 * g6)
pl = pl - 37 * cos(g5) + 25 * sin(2 * l5) + 25 * sin(3 * g5)
pl = pl - 23 * sin(g5 - 5 * g6) - 19 * t * cos(2 * g5 - 5 * g6)
pl = pl + 17 * cos(2 * g5 - 4 * g6) + 17 * cos(3 * g5 - 3 * g6)
pl = pl - 14 * sin(g5 - g6) - 13 * sin(3 * g5 - 4 * g6)
pl = pl - 9 * cos(2 * l5) + 9 * cos(g6) - 9 * sin(g6)
pl = pl - 9 * sin(3 * g5 - 2 * g6) + 9 * sin(4 * g5 - 5 * g6)
pl = pl + 9 * sin(2 * g5 - 6 * g6 + 3 * g7) - 8 * cos(4 * g5 - 10 * g6)
pl = pl + 7 * cos(3 * g5 - 4 * g6) - 7 * cos(g5 - 3 * g6)
pl = pl - 7 * sin(4 * g5 - 10 * g6) - 7 * sin(g5 - 3 * g6)
pl = pl + 6 * cos(4 * g5 - 5 * g6) - 6 * sin(3 * g5 - 3 * g6)
pl = pl + 5 * cos(2 * g6) - 4 * sin(4 * g5 - 4 * g6)
pl = pl - 4 * cos(3 * g6) + 4 * cos(2 * g5 - g6)
pl = pl - 4 * cos(3 * g5 - 2 * g6) - 4 * t * cos(2 * g5)
pl = pl + 3 * t * sin(2 * g5) + 3 * cos(5 * g6)
pl = pl + 3 * cos(5 * g5 - 10 * g6) + 3 * sin(2 * g6)
pl = pl - 2 * sin(2 * l5 - g5) + 2 * sin(2 * l5 + g5)
pl = pl - 2 * t * sin(3 * g5 - 5 * g6) - 2 * t * sin(g5 - 5 * g6)

plon = modulo(l5 + atr * pl)

plat = -4692 * cos(g5) + 259 * sin(g5) + 227 - 227 * cos(2 * g5)
plat = plat + 30 * t * sin(g5) + 21 * t * cos(g5)
plat = plat + 16 * sin(3 * g5 - 5 * g6) - 13 * sin(g5 - 5 * g6)
plat = plat - 12 * cos(3 * g5) + 12 * sin(2 * g5)
plat = plat + 7 * cos(3 * g5 - 5 * g6) - 5 * cos(g5 - 5 * g6)

plat = atr * plat

' heliocentric distance (au)

pr = 5.20883 - 0.25122 * cos(g5) - 0.00604 * cos(2 * g5)
pr = pr + 0.0026 * cos(2 * g5 - 2 * g6) - 0.0017 * cos(3 * g5 - 5 * g6)
pr = pr - 0.00106 * sin(2 * g5 - 2 * g6) - 0.00091 * t * sin(g5)
pr = pr - 0.00084 * t * cos(g5) + .00069 * sin(2 * g5 - 3 * g6)
pr = pr - 0.00067 * sin(g5 - 5 * g6) + 0.00066 * sin(3 * g5 - 5 * g6)
pr = pr + 0.00063 * sin(g5 - g6) - 0.00051 * cos(2 * g5 - 3 * g6)
pr = pr - 0.00046 * sin(g5) - 0.00029 * cos(g5 - 5 * g6)
pr = pr + 0.00027 * cos(g5 - 2 * g6) - 0.00022 * cos(3 * g5)
pr = pr - 0.00021 * sin(2 * g5 - 5 * g6)

end sub

''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''

sub saturn (jdate, plat, plon, pr)

' heliocentric coordinates of saturn

' input

' jdate = julian day

' output

' plat = heliocentric latitude (radians)
' plon = heliocentric longitude (radians)
' pr = heliocentric distance (au)

''''''''''''''''''''''''''''''''''''

atr = pi / 648000.0

' fundamental time arguments

djd = jdate - 2451545.0

t = djd / 36525.0 + 1.0

g5 = r2r(0.056531 + 0.00023080893 * djd)
l6 = r2r(0.133295 + 0.00009294371 * djd)
g6 = r2r(0.882987 + 0.00009294371 * djd)
g7 = r2r(0.400589 + 0.00003269438 * djd)

pl = 23045 * sin(g6) + 5014 * t - 2689 * cos(2 * g5 - 5 * g6) + 2507
pl = pl + 1177 * sin(2 * g5 - 5 * g6) - 826 * cos(2 * g5 - 4 * g6)
pl = pl + 802 * sin(2 * g6) + 425 * sin(g5 - 2 * g6)
pl = pl - 229 * t * cos(g6) - 153 * cos(2 * g5 - 6 * g6)
pl = pl - 142 * t * sin(g6) - 114 * cos(g6)
pl = pl + 101 * t * sin(2 * g5 - 5 * g6) - 70 * cos(2 * l6)
pl = pl + 67 * sin(2 * l6) + 66 * sin(2 * g5 - 6 * g6)
pl = pl + 60 * t * cos(2 * g5 - 5 * g6) + 41 * sin(g5 - 3 * g6)
pl = pl + 39 * sin(3 * g6) + 31 * sin(g5 - g6)
pl = pl + 31 * sin(2 * g5 - 2 * g6) - 29 * cos(2 * g5 - 3 * g6)
pl = pl - 28 * sin(2 * g5 - 6 * g6 + 3 * g7) + 28 * cos(g5 - 3 * g6)
pl = pl + 22 * t * sin(2 * g5 - 4 * g6) - 22 * sin(g6 - 3 * g7)
pl = pl + 20 * sin(2 * g5 - 3 * g6) + 20 * cos(4 * g5 - 10 * g6)
pl = pl + 19 * cos(2 * g6 - 3 * g7) + 19 * sin(4 * g5 - 10 * g6)
pl = pl - 17 * t * cos(2 * g6) - 16 * cos(g6 - 3 * g7)
pl = pl - 12 * sin(2 * g5 - 4 * g6) + 12 * cos(g5)
pl = pl - 12 * sin(2 * g6 - 2 * g7) - 11 * t * sin(2 * g6)
pl = pl - 11 * cos(2 * g5 - 7 * g6) + 10 * sin(2 * g6 - 3 * g7)
pl = pl + 10 * cos(2 * g5 - 2 * g6) + 9 * sin(4 * g5 - 9 * g6)
pl = pl - 8 * sin(g6 - 2 * g7) - 8 * cos(2 * l6 + g6)
pl = pl + 8 * cos(2 * l6 - g6) + 8 * cos(g6 - g7)
pl = pl - 8 * sin(2 * l6 - g6) + 7 * sin(2 * l6 + g6)
pl = pl - 7 * cos(g5 - 2 * g6) - 7 * cos(2 * g6)
pl = pl - 6 * t * sin(4 * g5 - 10 * g6) + 6 * cos(4 * g5 - 10 * g6) * t
pl = pl + 6 * t * sin(2 * g5 - 6 * g6) - 5 * sin(3 * g5 - 7 * g6)
pl = pl - 5 * cos(3 * g5 - 3 * g6) - 5 * cos(2 * g6 - 2 * g7)
pl = pl + 5 * sin(3 * g5 - 4 * g6) + 5 * sin(2 * g5 - 7 * g6)
pl = pl + 4 * sin(3 * g5 - 3 * g6) + 4 * sin(3 * g5 - 5 * g6)
pl = pl + 4 * cos(g5 - 2 * g6) * t + 3 * t * cos(2 * g5 - 4 * g6)
pl = pl + 3 * cos(2 * g5 - 6 * g6 + 3 * g7) - 3 * t * sin(2 * l6)
pl = pl + 3 * cos(2 * g5 - 6 * g6) * t - 3 * t * cos(2 * l6)
pl = pl + 3 * cos(3 * g5 - 7 * g6) + 3 * cos(4 * g5 - 9 * g6)
pl = pl + 3 * sin(3 * g5 - 6 * g6) + 3 * sin(2 * g5 - g6)
pl = pl + 3 * sin(g5 - 4 * g6) + 2 * cos(3 * g6 - 3 * g7)
pl = pl + 2 * t * sin(g5 - 2 * g6) + 2 * sin(4 * g6)
pl = pl - 2 * cos(3 * g5 - 4 * g6) - 2 * cos(2 * g5 - g6)
pl = pl - 2 * sin(2 * g5 - 7 * g6 + 3 * g7) + 2 * cos(g5 - 4 * g6)
pl = pl + 2 * cos(4 * g5 - 11 * g6) - 2 * sin(g6 - g7)

plon = modulo(l6 + atr * pl)

plat = 8297 * sin(g6) - 3346 * cos(g6) + 462 * sin(2 * g6)
plat = plat - 189 * cos(2 * g6) + 185 + 79 * t * cos(g6)
plat = plat - 71 * cos(2 * g5 - 4 * g6) + 46 * sin(2 * g5 - 6 * g6)
plat = plat - 45 * cos(2 * g5 - 6 * g6) + 29 * sin(3 * g6)
plat = plat - 20 * cos(2 * g5 - 3 * g6) + 18 * t * sin(g6)
plat = plat - 14 * cos(2 * g5 - 5 * g6) - 11 * cos(3 * g6) - 10 * t
plat = plat + 9 * sin(g5 - 3 * g6) + 8 * sin(g5 - g6)
plat = plat - 6 * sin(2 * g5 - 3 * g6) + 5 * sin(2 * g5 - 7 * g6)
plat = plat - 5 * cos(2 * g5 - 7 * g6) + 4 * sin(2 * g5 - 5 * g6)
plat = plat - 4 * sin(2 * g6) * t - 3 * cos(g5 - g6)
plat = plat + 3 * cos(g5 - 3 * g6) + 3 * t * sin(2 * g5 - 4 * g6)
plat = plat + 3 * sin(g5 - 2 * g6) + 2 * sin(4 * g6)
plat = plat - 2 * cos(2 * g5 - 2 * g6)

plat = atr * plat

pr = 9.55774 - 0.53252 * cos(g6) - 0.01878 * sin(2 * g5 - 4 * g6)
pr = pr - 0.01482 * cos(2 * g6) + 0.00817 * sin(g5 - g6)
pr = pr - 0.00539 * cos(g5 - 2 * g6) - 0.00524 * t * sin(g6)
pr = pr + 0.00349 * sin(2 * g5 - 5 * g6) + 0.00347 * sin(2 * g5 - 6 * g6)
pr = pr + 0.00328 * t * cos(g6) - 0.00225 * sin(g6)
pr = pr + 0.00149 * cos(2 * g5 - 6 * g6) - 0.00126 * cos(2 * g5 - 2 * g6)
pr = pr + 0.00104 * cos(g5 - g6) + 0.00101 * cos(2 * g5 - 5 * g6)
pr = pr + 0.00098 * cos(g5 - 3 * g6) - 0.00073 * cos(2 * g5 - 3 * g6)
pr = pr - 0.00062 * cos(3 * g6) + 0.00042 * sin(2 * g6 - 3 * g7)
pr = pr + 0.00041 * sin(2 * g5 - 2 * g6) - 0.0004 * sin(g5 - 3 * g6)
pr = pr + 0.0004 * cos(2 * g5 - 4 * g6) - 0.00028 * t - 0.00023 * sin(g5)
pr = pr + 0.0002 * sin(2 * g5 - 7 * g6)

end sub

''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''

sub neptune(jdate, plat, plon, pr)

' heliocentric coordinates of neptune

' input

' jdate = julian day

' output

' plat = heliocentric latitude (radians)
' plon = heliocentric longitude (radians)
' pr = heliocentric distance (au)

''''''''''''''''''''''''''''''''''''

atr = pi / 648000.0

' fundamental time arguments

djd = jdate - 2451545.0

tjd = djd / 36525.0 + 1.0

g5 = r2r(0.056531 + 0.00023080893 * djd)
g6 = r2r(0.882987 + 0.00009294371 * djd)
l7 = r2r(0.870169 + 0.00003269438 * djd)
g7 = r2r(0.400589 + 0.00003269438 * djd)
l8 = r2r(0.846912 + 0.00001672092 * djd)
g8 = r2r(0.725368 + 0.00001672092 * djd)
f8 = r2r(0.480856 + 0.00001663715 * djd)

pl = 3523 * sin(g8) - 50 * sin(2 * f8)
pl = pl - 43 * tjd * cos(g8) + 29 * sin(g5 - g8)
pl = pl + 19 * sin(2 * g8) - 18 * cos(g5 - g8)
pl = pl + 13 * cos(g6 - g8) + 13 * sin(g6 - g8)
pl = pl - 9 * sin(2 * g7 - 3 * g8) + 9 * cos(2 * g7 - 2 * g8)
pl = pl - 5 * cos(2 * g7 - 3 * g8) - 4 * tjd * sin(g8)
pl = pl + 4 * cos(g7 - 2 * g8) + 4 * tjd * tjd * sin(g8)

plon = modulo(l8 + atr * pl)

pb = 6404 * sin(f8) + 55 * sin(g8 + f8)
pb = pb + 55 * sin(g8 - f8) - 33 * tjd * sin(f8)

plat = atr * pb

pr = 30.07175 - 0.25701 * cos(g8)
pr = pr - 0.00787 * cos(2 * l7 - g7 - 2 * l8)
pr = pr + 0.00409 * cos(g5 - g8) - 0.00314 * tjd * sin(g8)
pr = pr + 0.0025 * sin(g5 - g8) - 0.00194 * sin(g6 - g8)
pr = pr + 0.00185 * cos(g6 - g8)

end sub

'''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''

sub uranus(jdate, plat, plon, pr)

' heliocentric coordinates of uranus

' input

' jdate = julian day

' output

' plat = heliocentric latitude (radians)
' plon = heliocentric longitude (radians)
' pr = heliocentric distance (au)

''''''''''''''''''''''''''''''''''''

atr = pi / 648000.0

' fundamental time arguments

djd = jdate - 2451545.0

tjd = djd / 36525.0 + 1.0

g5 = r2r(0.056531 + 0.00023080893 * djd)
g6 = r2r(0.882987 + 0.00009294371 * djd)
l7 = r2r(0.870169 + 0.00003269438 * djd)
g7 = r2r(0.400589 + 0.00003269438 * djd)
f7 = r2r(0.664614 + 0.00003265562 * djd)
g8 = r2r(0.725368 + 0.00001672092 * djd)

pl = 19397 * sin(g7) + 570 * sin(2 * g7)
pl = pl - 536 * tjd * cos(g7) + 143 * sin(g6 - 2 * g7)
pl = pl + 110 * tjd * sin(g7) + 102 * sin(g6 - 3 * g7)
pl = pl + 76 * cos(g6 - 3 * g7) - 49 * sin(g5 - g7)
pl = pl + 32 * tjd * tjd - 30 * tjd * cos(2 * g7)
pl = pl + 29 * sin(2 * g5 - 6 * g6 + 3 * g7)
pl = pl + 29 * cos(2 * g7 - 2 * g8) - 28 * cos(g7 - g8)
pl = pl + 23 * sin(3 * g7) - 21 * cos(g5 - g7)
pl = pl + 20 * sin(g7 - g8) + 20 * cos(g6 - 2 * g7)
pl = pl - 19 * cos(g6 - g7) + 17 * sin(2 * g7 - 3 * g8)
pl = pl + 14 * sin(3 * g7 - 3 * g8) + 13 * sin(g6 - g7)
pl = pl - 12 * tjd * tjd * cos(g7) - 12 * cos(g7)
pl = pl + 10 * sin(2 * g7 - 2 * g8) - 9 * sin(2 * f7)
pl = pl - 9 * tjd * tjd * sin(g7) + 9 * cos(2 * g7 - 3 * g8)
pl = pl + 8 * tjd * cos(g6 - 2 * g7) + 7 * tjd * cos(g6 - 3 * g7)
pl = pl - 7 * tjd * sin(g6 - 3 * g7) + 7 * tjd * sin(2 * g7)
pl = pl + 6 * sin(2 * g5 - 6 * g6 + 2 * g7)
pl = pl + 6 * cos(2 * g5 - 6 * g6 + 2 * g7)
pl = pl + 5 * sin(g6 - 4 * g7) - 4 * sin(3 * g7 - 4 * g8)
pl = pl + 4 * cos(3 * g7 - 3 * g8) - 3 * cos(g8)
pl = pl - 2 * sin(g8)

plon = modulo(l7 + atr * pl)

pb = 2775 * sin(f7) + 131 * sin(g7 - f7) + 130 * sin(g7 + f7)

plat = atr * pb

pr = 19.21216 - 0.90154 * cos(g7) - 0.02488 * tjd * sin(g7)
pr = pr - 0.02121 * cos(2 * g7) - 0.00585 * cos(g6 - 2 * g7)
pr = pr - 0.00508 * tjd * cos(g7) - 0.00451 * cos(g5 - g7)
pr = pr + 0.00336 * sin(g6 - g7) + 0.00198 * sin(g5 - g7)
pr = pr + 0.00118 * cos(g6 - 3 * g7) + 0.00107 * sin(g6 - 2 * g7)
pr = pr - 0.00103 * tjd * sin(2 * g7)
pr = pr - 0.00081 * cos(3 * g7 - 3 * g8)

end sub

''''''''''''''''''''''''''
''''''''''''''''''''''''''

sub pluto(jdtdb, pl, pb, pr)

' heliocentric coordinates of pluto subroutine

' input

' jdtdb = TDB julian day

' output

' pl = heliocentric longitude (radians)
' pb = heliocentric latitude (radians)
' pr = heliocentric distance (au)

''''''''''''''''''''''''''''''''''

' fundamental time argument

tjd = (jdtdb - 2451545.0) / 36525.0

' planetary longitudes

j = dtr * (34.35 + 3034.9057 * tjd)

s = dtr * (50.08 + 1222.1138 * tjd)

p = dtr * (238.96 + 144.96 * tjd)

' sum series

lonsum = 0.0

latsum = 0.0

rvsum = 0.0

for i% = 1 to 43

alpha = pj(i%) * j + ps(i%) * s + pp(i%) * p

salpha = sin(alpha)

calpha = cos(alpha)

lonsum = lonsum + plona(i%) * salpha + plonb(i%) * calpha

latsum = latsum + plata(i%) * salpha + platb(i%) * calpha

rvsum = rvsum + prva(i%) * salpha + prvb(i%) * calpha

next i%

' heliocentric longitude (radians)

pl = dtr * (238.956785 + 144.96 * tjd + 0.000001 * lonsum)

' heliocentric latitude (radians)

pb = dtr * (-3.908202 + 0.000001 * latsum)

' heliocentric distance (au)

pr = 40.7247248 + 0.0000001 * rvsum

end sub

''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''

sub julian(month, day, year, jday)

' Gregorian date to julian day subroutine

' input

' month = calendar month
' day = calendar day
' year = calendar year (all four digits)

' output

' jday = julian day

' special notes

' (1) calendar year must include all digits

' (2) will report October 5, 1582 to October 14, 1582
' as invalid calendar dates and exit

'''''''''''''''''''''''''''''''''''''''''

local a, b, c, m, y

y = year

m = month

b = 0.0

c = 0.0

if (m <= 2.0) then

y = y - 1.0

m = m + 12.0

end if

if (y < 0.0) then c = -0.75

if (year < 1582.0) then

' null

elseif (year > 1582.0) then

a = fix(y / 100.0)

b = 2.0 - a + fix(a / 4.0)

elseif (month < 10.0) then

' null

elseif (month > 10.0) then

a = fix(y / 100.0)

b = 2.0 - a + fix(a / 4.0)

elseif (day <= 4.0) then

' null

elseif (day > 14.0) then

a = fix(y / 100.0)

b = 2.0 - a + fix(a / 4.0)

else

print "this date does not exist!!"

exit

end if

jday = fix(365.25 * y + c) + fix(30.6001 * (m + 1.0)) + day + b + 1720994.5

end sub

''''''''''''''''''''''
''''''''''''''''''''''

function r2r(x) as float

' revolutions to radians function

' input

' x = argument (revolutions 0 <= x <= 1)

' output

' r2r = equivalent x (radians 0 <= y <= 2 pi)

'''''''''''''''''''''''''''''''''''''''''''''''

r2r = pi2 * (x - fix(x))

end function

'''''''''''''''''''''''''
'''''''''''''''''''''''''

function modulo(x) as float

' modulo 2 pi function

''''''''''''''''''''''

local a

a = x - pi2 * fix(x / pi2)

if (a < 0.0) then

a = a + pi2

end if

modulo = a

end function