Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me.
|
Forum Index : Microcontroller and PC projects : MMX - outer planet ephemeris
Author | Message | ||||
cdeagle Senior Member Joined: 22/06/2014 Location: United StatesPosts: 261 |
As promised, this post is software that computes an ephemeris for the outer planets. The computer program implements a "quick and dirty" ephemeris for Jupiter, Saturn, Neptune and Uranus, and a more sophisticated algorithm for Pluto. Here's a PDF document that probably has more than you would ever want to know about ephemerides. The last word in the previous sentence is fun to try and pronounce. 2017-04-13_192335_ephemeris.pdf Here's a typical program example for Jupiter and Pluto. julian day = 2446777.50000000 plat_jupiter = -1.2668 degrees plon_jupiter = 356.3351 degrees pr_jupiter = 4.96509278 au julian day = 2448908.50000000 plat_pluto = 232.740093 degrees plon_pluto = 14.587687 degrees pr_pluto = 29.71138254 au Here's the MMBASIC source code listing. ' 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 |
||||
Frank N. Furter Guru Joined: 28/05/2012 Location: GermanyPosts: 813 |
David, you made my day! My dream is to build a projection planetarium and I need your code to calculate the position of the planets. Now I should be able to program a Orrery as first step! THANKS AGAIN! Frank |
||||
Print this page |