.. a collection of basic astronomical functions for use in Euler .. program - written and tested on euler 1.36 comment ----------------------------------------------------- Astronomical functions 1999-08-19 Type 'help astro(RETURN)' for list of functions ----------------------------------------------------- endcomment function astro() ## Astronomical functions taken from ## Jean Meeus - 'Astronomical Algorithms' ## Montenbruck and Pfleger - 'Astronomy on the Personal Computer' ## Duffett-Smith - 'Practical astronomy with your calculator' ## other sources. ## For information on a function named fred, type 'help fred'. ## ## Report any problems or errors in these functions to ## Keith Burnett (keith@xylem.demon.co.uk) ## Latest version of this package is available from ## http://www.xylem.demon.co.uk/kepler/euler.html ## ## day jday gst nutation ## hmercury hvenus hearth hmars hjupiter hsaturn ## mercury venus sun mars jupiter saturn ## gmer gven gmar gjup ## gmoon amoon moon tmoon librate ## equatorial apparent mean altaz raltaz ## cart sph ## table ## ddegrees dsin dcos dtan dasin dacos datan datan2 ## brum reekie settle ## return 1; endfunction function ddegrees (d, m, s) ## converts degrees minutes and seconds to decimal degrees return d + m/60 + s/3600; endfunction .. Note that the built in mod function will give answers correct to .. 8 decimal places for angles up to 36,000,000 degrees - falls off .. for larger angles .. Note that Euler uses the same logical operators as C, so || is OR .. && is AND, ! is NOT. 'a NOT equal to b' is 'a != b', 'a equal to b' is .. 'a == b' and so on. function dsin (x) ## returns sine of x degrees return sin(pi()/180*mod(x, 360)); endfunction function dcos (x) ## returns cosine of x degrees return cos(pi()/180*mod(x, 360)); endfunction function dtan (x) ## returns tangent of x degrees return tan(pi()/180*mod(x, 360)); endfunction function dasin(x) ## returns angle in degrees corresponding to x return 180/pi()*asin(x) endfunction function dacos(x) ## returns angle in degrees corresponding to x return 180/pi()*acos(x) endfunction function datan(x) ## returns angle in degrees corresponding to x return 180/pi()*atan(x) endfunction function datan2(y, x) ## returns the angle in degrees within the correct ## quadrant given the coordinates y, x on the unit ## circle a = datan(y / x); if x < 0; a = a + 180; endif; if y < 0 && x > 0; a = a + 360; endif; return a; endfunction function day(y, m, d, h, min) ## returns the days since J2000.0 number assuming the Gregorian calendar ## after 1582 Oct 4th given the year, month, day, hour and minute ## This method is modified from Duffett-Smith Calculator page 7 ## Negative years CE are numbered according to the astronomical ## convention - i.e. calendrical year 1 BC is year zero in this function greg = y*10000 + m*100 + d; if m == 1 || m == 2; y = y - 1; m = m + 12; endif; if greg > 15821004; a = floor(y/100); b = 2 - a + floor(a/4); else; b = 0; endif; c = floor(365.25 * y); d1 = floor(30.6001 * (m + 1)); return b + c + d1 -730550.5 + d + (h+min/60)/24; endfunction function jday(y, m, d, h, min) ## returns the julian day number assuming the Gregorian calendar ## after 1582 Oct 4th given the year, month, day, hour and minute ## This method is taken from Duffett-Smith Calculator page 7 greg = y*10000 + m*100 + d; if m == 1 || m == 2; y = y - 1; m = m + 12; endif; if greg > 15821004; a = floor(y/100); b = 2 - a + floor(a/4); else; b = 0; endif; c = floor(365.25 * y); d1 = floor(30.6001 * (m + 1)); return b + c + d1 + 1720994.5 + d + (h+min/60)/24; endfunction function gst(day) ## returns the greenwich siderial time corresponding to the number of days ## before J2000.0 (Meeus Ch 11). ## Answer is given in degrees (360 = siderial day) t = day/36525; st = mod(280.46061837 + 360.98564736629 * day + 0.000387933 * t* t - t*t*t/38710000, 360); if st < 0; st = st + 360; endif; return st; endfunction function gmoon(day) ## returns the mean geocentric longitude, latitude and distance of the Moon ## given the instant in days since J2000.0 based on the truncated ELP-2000/82 ## theory in Meeus' book C45 - claimed good to 10 arcsec in longitude, ## 4 arcsec in latitude. .. Coefficients for mean longitude and radius vector of moon lcoef = [ 0 , 0 , 1 , 0 , 6288774 , -20905355; 2 , 0 , -1 , 0 , 1274027 , -3699111; 2 , 0 , 0 , 0 , 658314 , -2955968; 0 , 0 , 2 , 0 , 213618 , -569925; 0 , 1 , 0 , 0 , -185116 , 48888; 0 , 0 , 0 , 2 , -114332 , -3149; 2 , 0 , -2 , 0 , 58793 , 246158; 2 , -1 , -1 , 0 , 57066 , -152138; 2 , 0 , 1 , 0 , 53322 , -170733; 2 , -1 , 0 , 0 , 45758 , -204586; 0 , 1 , -1 , 0 , -40923 , -129620; 1 , 0 , 0 , 0 , -34720 , 108743; 0 , 1 , 1 , 0 , -30383 , 104755; 2 , 0 , 0 , -2 , 15327 , 10321; 0 , 0 , 1 , 2 , -12528 , 0; 0 , 0 , 1 , -2 , 10980 , 79661; 4 , 0 , -1 , 0 , 10675 , -34782; 0 , 0 , 3 , 0 , 10034 , -23210; 4 , 0 , -2 , 0 , 8548 , -21636; 2 , 1 , -1 , 0 , -7888 , 24208; 2 , 1 , 0 , 0 , -6766 , 30824; 1 , 0 , -1 , 0 , -5163 , -8379; 1 , 1 , 0 , 0 , 4987 , -16675; 2 , -1 , 1 , 0 , 4036 , -12831; 2 , 0 , 2 , 0 , 3994 , -10445; 4 , 0 , 0 , 0 , 3861 , -11650; 2 , 0 , -3 , 0 , 3665 , 14403; 0 , 1 , -2 , 0 , -2689 , -7003; 2 , 0 , -1 , 2 , -2602 , 0; 2 , -1 , -2 , 0 , 2390 , 10056; 1 , 0 , 1 , 0 , -2348 , 6322; 2 , -2 , 0 , 0 , 2236 , -9884; 0 , 1 , 2 , 0 , -2120 , 5751; 0 , 2 , 0 , 0 , -2069 , 0; 2 , -2 , -1 , 0 , 2048 , -4950; 2 , 0 , 1 , -2 , -1773 , 4130; 2 , 0 , 0 , 2 , -1595 , 0; 4 , -1 , -1 , 0 , 1215 , -3958; 0 , 0 , 2 , 2 , -1110 , 0; 3 , 0 , -1 , 0 , -892 , 3258; 2 , 1 , 1 , 0 , -810 , 2616; 4 , -1 , -2 , 0 , 759 , -1897; 0 , 2 , -1 , 0 , -713 , -2117; 2 , 2 , -1 , 0 , -700 , 2354; 2 , 1 , -2 , 0 , 691 , 0; 2 , -1 , 0 , -2 , 596 , 0; 4 , 0 , 1 , 0 , 549 , -1423; 0 , 0 , 4 , 0 , 537 , -1117; 4 , -1 , 0 , 0 , 520 , -1571; 1 , 0 , -2 , 0 , -487 , -1739; 2 , 1 , 0 , -2 , -399 , 0; 0 , 0 , 2 , -2 , -381 , -4421; 1 , 1 , 1 , 0 , 351 , 0; 3 , 0 , -2 , 0 , -340 , 0; 4 , 0 , -3 , 0 , 330 , 0; 2 , -1 , 2 , 0 , 327 , 0; 0 , 2 , 1 , 0 , -323 , 1165; 1 , 1 , -1 , 0 , 299 , 0; 2 , 0 , 3 , 0 , 294 , 0; 2 , 0 , -1 , -2 , 0 , 8752 ]; .. Coeficients for mean latitude of Moon bcoef = [ 0 , 0 , 0 , 1 , 5128122; 0 , 0 , 1 , 1 , 280602; 0 , 0 , 1 , -1 , 277693; 2 , 0 , 0 , -1 , 173237; 2 , 0 , -1 , 1 , 55413; 2 , 0 , -1 , -1 , 46271; 2 , 0 , 0 , 1 , 32573; 0 , 0 , 2 , 1 , 17198; 2 , 0 , 1 , -1 , 9266; 0 , 0 , 2 , -1 , 8822; 2 , -1 , 0 , -1 , 8216; 2 , 0 , -2 , -1 , 4324; 2 , 0 , 1 , 1 , 4200; 2 , 1 , 0 , -1 , -3359; 2 , -1 , -1 , 1 , 2463; 2 , -1 , 0 , 1 , 2211; 2 , -1 , -1 , -1 , 2065; 0 , 1 , -1 , -1 , -1870; 4 , 0 , -1 , -1 , 1828; 0 , 1 , 0 , 1 , -1794; 0 , 0 , 0 , 3 , -1749; 0 , 1 , -1 , 1 , -1565; 1 , 0 , 0 , 1 , -1491; 0 , 1 , 1 , 1 , -1475; 0 , 1 , 1 , -1 , -1410; 0 , 1 , 0 , -1 , -1344; 1 , 0 , 0 , -1 , -1335; 0 , 0 , 3 , 1 , 1107; 4 , 0 , 0 , -1 , 1021; 4 , 0 , -1 , 1 , 833; 0 , 0 , 1 , -3 , 777; 4 , 0 , -2 , 1 , 671; 2 , 0 , 0 , -3 , 607; 2 , 0 , 2 , -1 , 596; 2 , -1 , 1 , -1 , 491; 2 , 0 , -2 , 1 , -451; 0 , 0 , 3 , -1 , 439; 2 , 0 , 2 , 1 , 422; 2 , 0 , -3 , -1 , 421; 2 , 1 , -1 , 1 , -366; 2 , 1 , 0 , 1 , -351; 4 , 0 , 0 , 1 , 331; 2 , -1 , 1 , 1 , 315; 2 , -2 , 0 , -1 , 302; 0 , 0 , 1 , 3 , -283; 2 , 1 , 1 , -1 , -229; 1 , 1 , 0 , -1 , 223; 1 , 1 , 0 , 1 , 223; 0 , 1 , -2 , -1 , -220; 2 , 1 , -1 , -1 , -220; 1 , 0 , 1 , 1 , -185; 2 , -1 , -2 , -1 , 181; 0 , 1 , 2 , 1 , -177; 4 , 0 , -2 , -1 , 176; 4 , -1 , -1 , -1 , 166; 1 , 0 , 1 , -1 , -164; 4 , 0 , 1 , -1 , 132; 1 , 0 , -1 , -1 , -119; 4 , -1 , 0 , -1 , 115; 2 , -2 , 0 , 1 , 107]; ..Calculate the arguments, eccentricity and additive corrections t = day/36525; t2 = t*t; t3 = t2*t; t4 = t2*t2; .. Mean longitude including light travel time and referred to equinox of date l1 = mod(218.3164591 + 481267.88134236 *t - 0.0013268 *t2 + t3/538841 - t4/65194000, 360); .. Mean elongation d = mod(297.8502042 + 445267.1115168 *t - 0.0016300 *t2 + t3/545868 - t4/113065000, 360); .. Sun's mean anomaly m = mod(357.5291092 + 35999.0502909 *t - 0.0001536 *t2 + t3/24490000, 360); .. Moon's mean anomaly m1 = mod(134.9634114 + 477198.8676313 *t + 0.0089970 *t2 + t3/69699 - t4/14712000, 360); ..Moon's argument of latitude (mean distance from ascending node) f = mod(93.2720993 + 483202.0175273 *t - 0.0034029 *t2 - t3/3526000 + t4/863310000, 360); ..further arguments a1 = mod(119.75 + 131.849 * t,360); a2 = mod( 53.09 + 479264.290 * t,360); a3 = mod(313.45 + 481266.484 * t,360); .. Eccentricity of earth's orbit round Sun (affects terms with M or 2M) e = 1 - 0.002516 * t - 0.0000074 *t2; e2 = e * e; .. Use matrix algebra to sum the series (slight differences compared with .. for loop summation, around 0.1 arcsec .. Longitude and radius vector series tr = lcoef'; .. Condition coefficents with eccentricity correction for inc = 1 to cols(tr); if abs(tr[2, inc]) == 1; tr[5, inc] = e*tr[5, inc]; tr[6, inc] = e*tr[6, inc]; endif; if abs(tr[2, inc]) == 2; tr[5, inc] = e2*tr[5, inc]; tr[6, inc] = e2*tr[6, inc]; endif; end; .. Form the term vectors and sum the series arg = mod(tr[1]*d + tr[2]*m + tr[3]*m1 + tr[4]*f, 360); long = sum(tr[5] * dsin(arg)); r = sum(tr[6] * dcos(arg)); .. Latitude series tr = bcoef'; for inc = 1 to cols(tr); if abs(tr[2, inc]) == 1; tr[5, inc] = e*tr[5, inc]; endif; if abs(tr[2, inc]) == 2; tr[5, inc] = e2*tr[5, inc]; endif; end; arg = mod(tr[1]*d + tr[2]*m + tr[3]*m1 + tr[4]*f, 360); lat = sum(tr[5] * dsin(arg)); .. additive terms for longitude long = long + 3958 * dsin(a1) + 1962 * dsin (l1 - f) + 318 * dsin(a2); lambda = l1 + long/1000000; if lambda < 0; lambda = lambda + 360; endif; .. additive terms for latitude lat = lat - 2235 * dsin(l1) + 382 * dsin(a3) + 175 * dsin(a1 - f); lat = lat + 175 * dsin(a1 + f) + 127 * dsin(l1 - m1) - 115 * dsin(l1 + m1); beta = lat/1000000; .. calculate radius vector in Km delta = 385000.56 + r/1000; return [lambda, beta, delta]; endfunction function nutation(day) ## returns the values of delta-phi and delta-epsilon in degrees ## for UT instant. See Meeus Ch21 - good to 0.5 arcsec in phi and ## 0.1 arcsec in epsilon ## Returns a vector with entries [dphi, deps, 0] for compatibility t = day/36525; t2 = t*t; t3 = t2*t; .. Arguments .. Mean longitude of the Sun l = mod(280.4665 + 36000.7698 * t, 360); .. Mean longitude of the Moon l1 = mod(218.3165 + 481267.8813 * t, 360); .. Longitude of the ascending node of the Moon's orbit on Ecliptic plane, measured from .. mean equinox of UT instant o = mod(125.04452 - 1934.136261 * t - 0.0020708 * t2 + t3/327270, 360); .. Series dphi = -17.20 * dsin(o) - 1.32 * dsin(2*l) - 0.23 * dsin(2*l1) + 0.21 * dsin(2 * o); dphi = dphi/3600; deps = 9.20 * dcos(o) + 0.57 * dcos(2*l) + 0.10 * dcos(2*l1) - 0.09 * dcos(2 * o); deps = deps/3600; return [dphi, deps, 0]; endfunction function amoon(day) ## returns the apparent geocentric longitude, latitude and distance of the Moon ## corrected for nutation using a low precision nutation theory meanmoon = gmoon(day); nutcor = nutation(day); nutcor[2] = 0; return meanmoon + nutcor; endfunction function table(position, day1, day2, inc) ## builds a table of positions given the position function name, ## a starting UT instant, a stopping UT instant and a time increment. ## Each row of the table is a position. ..start = time c = 1; moontab = zeros([(day2-day1)/inc+1, 3]); for day = day1 to (day2 + inc) step inc; moonnow = position(day); moontab(c, 1) = moonnow(1); moontab(c, 2) = moonnow(2); moontab(c, 3) = moonnow(3); c = c+1; end; ..finish = time return moontab endfunction function hearth(day) ## returns the heliocentric ecliptic latitude of the Earth given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC .. define coeficient tables for the constants l0, l1 and so on l0 = [175347046 , 0.0000000 , 0.0000000; 3341656 , 4.6692568 , 6283.0758500; 34894 , 4.6261024 , 12566.1517000; 3497 , 2.7441178 , 5753.3848849; 3418 , 2.8288658 , 3.5231183; 3136 , 3.6276704 , 77713.7714681; 2676 , 4.4180835 , 7860.4193924; 2343 , 6.1351621 , 3930.2096962; 1324 , 0.7424634 , 11506.7697698; 1273 , 2.0370966 , 529.6909651; 1199 , 1.1096295 , 1577.3435424; 990 , 5.2326807 , 5884.9268466; 902 , 2.0450545 , 26.2983198; 857 , 3.5084915 , 398.1490034; 780 , 1.1788268 , 5223.6939198; 753 , 2.5333905 , 5507.5532387; 505 , 4.5829260 , 18849.2275500; 492 , 4.2050571 , 775.5226113; 357 , 2.9195411 , 0.0673103; 317 , 5.8490195 , 11790.6290887; 284 , 1.8986924 , 796.2980068; 271 , 0.3148626 , 10977.0788047; 243 , 0.3448145 , 5486.7778432; 206 , 4.8064663 , 2544.3144199; 205 , 1.8695377 , 5573.1428014; 202 , 2.4576779 , 6069.7767546; 156 , 0.8330608 , 213.2990954; 132 , 3.4111829 , 2942.4634233; 126 , 1.0829546 , 20.7753955; 115 , 0.6454491 , 0.9803211; 103 , 0.6359985 , 4694.0029547; 102 , 0.9756928 , 15720.8387849; 102 , 4.2667980 , 7.1135470; 99 , 6.2099293 , 2146.1654165; 98 , 0.6810134 , 155.4203994; 86 , 5.9832263 , 161000.6857377; 85 , 1.2987076 , 6275.9623030; 85 , 3.6708009 , 71430.6956181; 80 , 1.8079129 , 17260.1546547; 79 , 3.0369746 , 12036.4607349; 75 , 1.7550891 , 5088.6288398; 74 , 3.5031941 , 3154.6870849; 74 , 4.6792663 , 801.8209311; 70 , 0.8329762 , 9437.7629349; 62 , 3.9776391 , 8827.3902699; 61 , 1.8183989 , 7084.8967811; 57 , 2.7843046 , 6286.5989683; 56 , 4.3869487 , 14143.4952424; 56 , 3.4700606 , 6279.5527316; 52 , 0.1891495 , 12139.5535091; 52 , 1.3328274 , 1748.0164131; 51 , 0.2830683 , 5856.4776591; 49 , 0.4873501 , 1194.4470102; 41 , 5.3681759 , 8429.2412665; 41 , 2.3985094 , 19651.0484811; 39 , 6.1683302 , 10447.3878396; 37 , 6.0413386 , 10213.2855462; 37 , 2.5695748 , 1059.3819302; 36 , 1.7087581 , 2352.8661538; 36 , 1.7759689 , 6812.7668151; 33 , 0.5931028 , 17789.8456198; 30 , 0.4429446 , 83996.8473181; 30 , 2.7397512 , 1349.8674097; 25 , 3.1647089 , 4690.4798364 ]; l1 = [628331966747 0 0; 206059 2.678234558 6283.07585; 4303 2.635122335 12566.1517; 425 1.59046982 3.523118349; 119 5.795557656 26.2983198; 109 2.966310107 1577.343542; 93 2.592111095 18849.22755; 72 1.138405812 529.6909651; 68 1.874533003 398.1490034; 67 4.40932832 5507.553239; 59 2.888157906 5223.69392; 56 2.1747174 155.4203994; 45 0.397995029 796.2980068; 36 0.468754372 775.5226113; 29 2.647322546 7.113547001; 21 5.341382751 0.980321068; 19 1.84628376 5486.777843; 19 4.968551795 213.2990954; 17 2.991167606 6275.962303; 16 0.032165873 2544.31442; 16 1.430493013 2146.165416; 15 1.204697937 10977.0788; 12 2.834322821 1748.016413; 12 3.25805082 5088.62884; 12 5.273797604 1194.44701; 12 2.075020801 4694.002955; 11 0.76614723 553.5694028; 10 1.302634234 6286.598968; 10 4.239258653 1349.86741; 9 2.69956827 242.728604; 9 5.64476086 951.7184063; 8 5.300561729 2352.866154; 6 2.65034514 9437.762935; 6 4.666337263 4690.479836 ]; l2 = [52919 0 0; 8720 1.0721 6283.0758; 309 0.867 12566.152; 27 0.05 3.52; 16 5.19 26.3; 16 3.68 155.42; 10 0.76 18849.23; 9 2.06 77713.77; 7 0.83 775.52; 5 4.66 1577.34; 4 1.03 7.11; 4 3.44 5573.14; 3 5.14 796.3; 3 6.05 5507.55; 3 1.19 242.73; 3 6.12 529.69; 3 0.31 398.15; 3 2.28 553.57; 2 4.38 5223.69; 2 3.75 0.98 ]; l3 = [289 5.844 6283.076; 35 0 0; 17 5.49 12566.15; 3 5.2 155.42; 1 4.72 3.52; 1 5.3 18849.23; 1 5.97 242.73 ]; l4 = [114 3.142 0; 8 4.13 6283.08; 1 3.84 12556.15 ]; l5 = [1 3.14 0 ]; .. latitude terms - Meeus truncates most of these b0 = [ 280 3.19870156 84334.66158; 102 5.422486193 5507.553239; 80 3.880132045 5223.69392; 44 3.704446898 2352.866154; 32 4.000263698 1577.343542 ]; b1 = [ 9 3.90 5507.55; 6 1.73 5223.69 ]; .. Radius vector terms r0 = [ 100013989 0 0; 1670700 3.098463508 6283.07585; 13956 3.055246096 12566.1517; 3084 5.198466744 77713.77147; 1628 1.17387749 5753.384885; 1576 2.846852458 7860.419392; 925 5.452922341 11506.76977; 542 4.564091498 3930.209696; 472 3.661000221 5884.926847; 346 0.963686177 5507.553239; 329 5.899836465 5223.69392; 307 0.298671395 5573.142801; 243 4.273495362 11790.62909; 212 5.847145403 1577.343542; 186 5.021944472 10977.0788; 175 3.011936365 18849.22755; 110 5.055106363 5486.777843; 98 0.886813113 6069.776755; 86 5.689597783 15720.83878; 86 1.270837334 161000.6857; 65 0.272506138 17260.15465; 63 0.921771088 529.6909651; 57 2.01374292 83996.84732; 56 5.241597989 71430.69562; 49 3.245012404 2544.31442; 47 2.578050704 775.5226113; 45 5.537158073 9437.762935; 43 6.01110242 6275.962303; 39 5.360717382 4694.002955; 38 2.39255344 8827.39027; 37 0.829529223 19651.04848; 37 4.901075919 12139.55351; 36 1.67468059 12036.46073; 35 1.842706933 2942.463423; 33 0.243703001 7084.896781; 32 0.183682298 5088.62884; 32 1.777756421 398.1490034; 28 1.213448682 6286.598968; 28 1.899343309 6279.552732; 26 4.588968504 10447.38784 ]; r1 = [ 103019 1.107490 6283.075850; 1721 1.0644 12566.1517; 702 3.142 0; 32 1.02 18849.23; 31 2.84 5507.55; 25 1.32 5223.69; 18 1.42 1577.34; 10 5.91 10977.08; 9 1.42 6275.96; 9 0.27 5486.78 ]; r2 = [ 4359 5.7846 6283.0758; 124 5.579 12566.152; 12 3.14 0; 9 3.63 77713.77; 6 1.87 5573.14; 3 5.47 18849.23 ]; r3 = [ 145 4.273 6283.076; 7 3.92 12566.15 ]; r4 = [ 4 2.56 6283.076 ]; .. Now work out the sums of the terms t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); beta = (sb0 + sb1*t)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); r = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4)/100000000; return [lambda, beta, r]; endfunction function hvenus(day) ## returns the heliocentric ecliptic latitude of Venus given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC .. Longitude coefficents l0 = [ 317614667 0 0; 1353968 5.593133196 10213.28555; 89892 5.306500485 20426.57109; 5477 4.416306525 7860.419392; 3456 2.699644708 11790.62909; 2372 2.993775396 3930.209696; 1664 4.25018935 1577.343542; 1438 4.15745044 9683.594581; 1317 5.186682191 26.2983198; 1201 6.153571153 30639.85664; 769 0.816296159 9437.762935; 761 1.950147021 529.6909651; 708 1.064667072 775.5226113; 585 3.998398848 191.4482661; 500 4.123402101 15720.83878; 429 3.586428598 19367.18916; 327 5.677365837 5507.553239; 326 4.590564731 10404.73381; 232 3.162510571 9153.903616; 180 4.653379156 1109.378552; 155 5.570438889 19651.04848; 128 4.226044937 20.77539549; 128 0.962098227 5661.332049; 106 1.537211913 801.8209311 ]; l1 = [ 1021352943053 0 0; 95708 2.46424449 10213.28555; 14445 0.516245647 20426.57109; 213 1.795479294 30639.85664; 152 6.106352824 1577.343542; 174 2.655358794 26.2983198; 82 5.702341337 191.4482661; 70 2.68136035 9437.762935; 52 3.600130877 775.5226113; 38 1.03379038 529.6909651; 30 1.250563224 5507.553239; 25 6.106647929 10404.73381 ]; l2 = [ 54127 0 0; 3891 0.3451436 10213.28555; 1338 2.020112861 20426.57109; 24 2.04592119 26.2983198; 19 3.535273715 30639.85664; 10 3.971302211 775.5226113; 7 1.519625934 1577.343542; 6 0.999267579 191.4482661 ]; l3 = [ 136 4.80389021 10213.28555; 78 3.668763716 20426.57109; 26 0 0 ]; l4 = [ 114, 3.1416, 0; 3, 5.21, 20426.57; 2, 2.51, 10213.29 ]; l5 = [ 1, 3.14, 0 ]; .. latitude coefficients; b0 = [ 5923638 0.2670278 10213.2855462; 40108 1.14737 20426.57109; 32815 3.14159 0; 1011 1.0895 30639.8566; 149 6.254 18073.705; 138 0.860 1577.344; 130 3.672 9437.763; 120 3.705 2352.866; 108 4.539 22003.915 ]; b1 = [ 513348 1.803643 10213.285546; 4380 3.3862 20426.5711; 199 0 0; 197 2.530 30639.857 ]; b2 = [ 22378 3.38509 10213.28555; 282 0 0; 173 5.256 20426.571; 27 3.87 30639.86 ]; b3 = [ 647 4.992 10213.286; 20 3.14 0; 6 0.77 20426.57; 3 5.44 30639.86 ]; b4 = [ 14 0.32 10213.29 ]; .. Radius vector coefficients r0 = [ 72334821 0 0; 489824 4.021518 10213.285546; 1658 4.9021 20426.5711; 1632 2.8455 7860.4194; 1378 1.1285 11790.6291; 498 2.587 9683.595; 374 1.423 3930.210; 264 5.529 9437.763; 237 2.551 15720.839; 222 2.013 19367.189; 126 2.768 1577.344; 119 3.020 10404.734 ]; r1 = [ 34551 0.89199 10213.28555; 234 1.772 20426.571; 234 3.142 0 ]; r2 = [ 1407 5.0637 10213.2855; 16 5.47 20426.57; 13 0 0 ]; r3 = [ 50 3.22 10213.29 ]; r4 = [ 1 0.92 10213.29 ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4)/100000000; return [lambda, beta, delta]; endfunction function zzterm(a, t) ## calculates the value of a periodic term in the VSOP82 analytical theory ## for the position of the planets - called by the planet position functions tpi = 2*pi(); tr = a'; vec1 = tr[1] * cos(mod(tr[2] + tr[3] * t, tpi)); total = sum(vec1); return total; endfunction function hjupiter(day) ## returns the heliocentric ecliptic latitude of Jupiter given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## Accuracy of order 4 arcsec in longitude .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC l0 = [ 59954691 0 0; 9695899 5.061917931 529.6909651; 573610 1.44406206 7.113547001; 306389 5.4173473 1059.38193; 97178 4.142647088 632.7837393; 72903 3.640429093 522.5774181; 64264 3.411451852 103.0927742; 39806 2.293767449 419.4846439; 38858 1.272317249 316.3918697; 27965 1.784545895 536.8045121; 13590 5.774810316 1589.072895; 8769 3.630003244 949.175609; 8246 3.582279617 206.1855484; 7368 5.081011256 735.8765135; 6263 0.024976437 213.2990954; 6114 4.513195317 1162.474704; 5305 4.186250535 1052.268383; 5305 1.306712368 14.227094; 4905 1.320846317 110.2063212; 4647 4.699581095 3.932153263; 3045 4.316759603 426.5981909; 2610 1.566675949 846.0828348; 2028 1.063765474 3.181393738; 1921 0.971689288 639.8972863; 1765 2.141480778 1066.495477; 1723 3.880360089 1265.567479; 1633 3.582010898 515.4638711; 1432 4.296836903 625.6701923; 973 4.097649571 95.97922722; 884 2.437014261 412.3710969; 733 6.085341132 838.9692878; 731 3.80591234 1581.959348; 709 1.292725737 742.9900605; 692 6.133682229 2118.76386; 614 4.108534968 1478.866574; 582 4.539677176 309.2783227; 495 3.755674614 323.5054167; 441 2.958184609 454.9093665; 417 1.035544302 2.447680555; 390 4.897161059 1692.16567; 376 4.702991248 1368.660253; 341 5.714525258 533.6231184; 330 4.740498195 0.04818411; 262 1.87652461 0.963207847; 261 0.820472464 380.127768; 257 3.724107242 199.0720014; 244 5.220208789 728.7629665; 235 1.226939081 909.8187331; 220 1.65115016 543.9180591; 207 1.854616666 525.7588118; 202 1.806845742 1375.7738; 197 5.29252149 1155.361157; 175 3.729665548 942.062062; 175 3.226349034 1898.351218; 175 5.909735053 956.289156; 158 4.364839218 1795.258444; 151 3.906250226 74.78159857; 149 4.377451043 1685.052123; 141 3.135683579 491.5579295; 138 1.317979208 1169.588251; 131 4.168679455 1045.154836; 117 2.500221409 1596.186442; 117 3.38920921 0.521264862; 106 4.554397982 526.5095714 ]; l1 = [ 52993480757 0 0; 489741 4.220666899 529.6909651; 228919 6.02647464 7.113547001; 27655 4.572659568 1059.38193; 20721 5.459389363 522.5774181; 12106 0.16985765 536.8045121; 6068 4.42419502 103.0927742; 5434 3.984783826 419.4846439; 4238 5.890093513 14.227094; 2212 5.267714466 206.1855484; 1746 4.926693785 1589.072895; 1296 5.551327651 3.181393738; 1173 5.856473044 1052.268383; 1163 0.514508953 3.932153263; 1099 5.307049816 515.4638711; 1007 0.464783986 735.8765135; 1004 3.150403018 426.5981909; 848 5.758058505 110.2063212; 827 4.803120157 213.2990954; 816 0.586430549 1066.495477; 725 5.518274715 639.8972863; 568 5.988670495 625.6701923; 474 4.132452692 412.3710969; 413 5.736528913 95.97922722; 345 4.241595654 632.7837393; 336 3.73248749 1162.474704; 234 4.034699703 949.175609; 234 6.243022266 309.2783227; 199 1.504584428 838.9692878; 195 2.218790109 323.5054167; 187 6.086205659 742.9900605; 184 6.279635888 543.9180591; 171 5.416559838 199.0720014; 131 0.626433774 728.7629665; 115 0.680190502 846.0828348; 115 5.286416991 2118.76386; 108 4.492827601 956.289156; 80 5.824124003 1045.154836; 72 5.341626503 942.062062; 70 5.972634503 532.8723588; 67 5.733651265 21.340641; 66 0.129241914 526.5095714; 65 6.088034903 1581.959348; 59 0.58626971 1155.361157; 58 0.994530873 1596.186442; 57 5.968513048 1169.588251; 57 1.411984388 533.6231184; 55 5.428063837 10.29494074; 52 5.726614484 117.3198682; 52 0.229812991 1368.660253; 50 6.080751478 525.7588118; 47 3.626118432 1478.866574; 47 0.511440732 1265.567479; 40 4.161580136 1692.16567; 34 0.099139049 302.1647757; 33 5.035966895 220.4126424; 32 5.374925307 508.3503241; 29 5.422088971 1272.681026; 29 3.359272415 4.665866446; 29 0.759079097 88.86568022; 25 1.607230634 831.8557407 ]; l2 = [ 47234 4.321483236 7.113547001; 38966 0 0; 30629 2.930214402 529.6909651; 3189 1.055046156 522.5774181; 2729 4.845454814 536.8045121; 2723 3.414115266 1059.38193; 1721 4.187343852 14.227094; 383 5.767907144 419.4846439; 378 0.760489649 515.4638711; 367 6.055091204 103.0927742; 337 3.786443842 3.181393738; 308 0.693566541 206.1855484; 218 3.813891914 1589.072895; 199 5.339964434 1066.495477; 197 2.483564021 3.932153263; 156 1.406424265 1052.268383; 146 3.813731968 639.8972863; 142 1.63435169 426.5981909; 130 5.837388725 412.3710969; 117 1.414354626 625.6701923; 97 4.033834279 110.2063212; 91 1.10630629 95.97922722; 87 2.522351748 632.7837393; 79 4.637261313 543.9180591; 72 2.2171667 735.8765135; 58 0.832163174 199.0720014; 57 3.122920599 213.2990954; 49 1.672837916 309.2783227; 40 4.024854447 21.340641; 40 0.624169458 323.5054167; 36 2.32581247 728.7629665; 29 3.608383278 10.29494074; 28 3.239920137 838.9692878; 26 4.501182983 742.9900605; 26 2.512406239 1162.474704; 25 1.218681107 1045.154836; 24 3.005321393 956.289156; 19 4.290286447 532.8723588; 18 0.809539416 508.3503241; 17 4.200019777 2118.76386; 17 1.834021466 526.5095714; 15 5.810379869 1596.186442; 15 0.681741655 942.062062; 15 3.999896226 117.3198682; 14 5.951695685 316.3918697; 14 1.80336678 302.1647757; 13 2.518566436 88.86568022; 13 4.368562324 1169.588251; 11 4.435866346 525.7588118; 10 1.715631611 1581.959348; 9 2.176845635 1155.361157; 9 3.294527833 220.4126424; 9 3.319244936 831.8557407; 8 5.756722284 846.0828348; 8 2.709555168 533.6231184; 7 2.175600933 1265.567479; 6 0.499398635 949.175609 ]; l3 = [ 6502 2.598628805 7.113547001; 1357 1.346358864 529.6909651; 471 2.475039779 14.227094; 417 3.244512432 536.8045121; 353 2.97360159 522.5774181; 155 2.075655858 1059.38193; 87 2.514315843 515.4638711; 44 0 0; 34 3.826337945 1066.495477; 28 2.447547561 206.1855484; 24 1.276671723 412.3710969; 23 2.982313268 543.9180591; 20 2.10099934 639.8972863; 20 1.40255939 419.4846439; 19 1.593684035 103.0927742; 17 2.302146812 21.340641; 17 2.598214607 1589.072895; 16 3.145211173 625.6701923; 16 3.360301263 1052.268383; 13 2.759738922 95.97922722; 13 2.538622443 199.0720014; 13 6.265781104 426.5981909; 9 1.763349607 10.29494074; 9 2.265632563 110.2063212; 7 3.425664333 309.2783227; 7 4.038695629 728.7629665; 6 2.520964177 508.3503241; 5 2.911846871 1045.154836; 5 5.251961535 323.5054167; 4 4.302902612 88.86568022; 4 3.523813616 302.1647757; 4 4.091253151 735.8765135; 3 1.431759913 956.289156; 3 4.358175077 1596.186442; 3 1.252765908 213.2990954; 3 5.015058398 838.9692878; 3 2.237856733 117.3198682; 2 2.896624092 742.9900605; 2 2.355818712 942.062062 ]; l4 = [ 669 0.852824211 7.113547001; 114 3.141592654 0; 100 0.742589478 14.227094; 50 1.653462082 536.8045121; 44 5.820263866 529.6909651; 32 4.858299867 522.5774181; 15 4.290616358 515.4638711; 9 0.714785207 1059.38193; 5 1.295022594 543.9180591; 4 2.317155166 1066.495477; 4 0.483267975 21.340641; 3 3.002455427 412.3710969; 2 0.398589402 639.8972863; 2 4.259256203 199.0720014; 2 4.905362073 625.6701923; 2 4.261475808 206.1855484; 1 5.255469557 1052.268383; 1 4.716146338 95.97922722; 1 1.286045712 1589.072895 ]; l5 = [ 50 5.26 7.11; 16 5.25 14.23; 4 0.01 536.80; 2 1.10 522.58; 1 3.14 0 ]; .. latitude series b0 = [ 2268616 3.558526067 529.6909651; 110090 0 0; 109972 3.908093474 1059.38193; 8101 3.605095734 522.5774181; 6438 0.306271214 536.8045121; 6044 4.258831088 1589.072895; 1107 2.985344219 1162.474704; 944 1.675222884 426.5981909; 942 2.936190724 1052.268383; 894 1.754474299 7.113547001; 836 5.178819732 103.0927742; 767 2.154735941 632.7837393; 684 3.678087701 213.2990954; 629 0.643432823 1066.495477; 559 0.013548305 846.0828348; 532 2.703059544 110.2063212; 464 1.173372492 949.175609; 431 2.608250005 419.4846439; 351 4.610629907 2118.76386; 132 4.778169907 742.9900605; 123 3.349681814 1692.16567; 116 1.38688232 323.5054167; 115 5.048922954 316.3918697; 104 3.701038381 515.4638711; 103 2.318789996 1478.866574; 102 3.152937854 1581.959348 ]; b1 = [ 177352 5.701664885 529.6909651; 3230 5.779416193 1059.38193; 3081 5.474642965 522.5774181; 2212 4.734774802 536.8045121; 1694 3.141592654 0; 346 4.745951741 1052.268383; 234 5.188560999 1066.495477; 196 6.185542866 7.113547001; 150 3.927212261 1589.072895; 114 3.438972718 632.7837393; 97 2.914263041 949.175609; 82 5.076660975 1162.474704; 77 2.505221887 103.0927742; 77 0.612889814 419.4846439; 74 5.499582922 515.4638711; 61 5.447400844 213.2990954; 50 3.947996166 735.8765135; 46 0.538503609 110.2063212; 45 1.895166452 846.0828348; 37 4.698283928 543.9180591; 36 6.109525788 316.3918697; 32 4.92 1581.96 ]; b2 = [ 8094 1.463228437 529.6909651; 813 3.141592654 0; 742 0.95691639 522.5774181; 399 2.898886664 536.8045121; 342 1.446837897 1059.38193; 74 0.407246759 1052.268383; 46 3.480368958 1066.495477; 30 1.925041713 1589.072895; 29 0.990888318 515.4638711; 23 4.271240524 7.113547001; 14 2.922423873 543.9180591; 12 5.221689325 632.7837393; 11 4.880242225 949.175609; 6 6.210891084 1045.154836 ]; b3 = [ 252 3.381 529.691; 122 2.733 522.577; 49 1.04 536.80; 11 2.31 1052.27; 8 2.77 515.46; 7 4.25 1059.38; 6 1.78 1066.50; 4 1.13 543.92; 3 3.14 0 ]; b4 = [ 15 4.53 522.58; 5 4.47 529.69; 4 5.44 536.80; 3 0 0; 2 4.52 515.46; 1 4.20 1052.27 ]; b5 = [ 1 0.09 522.58 ]; .. radius vector r0 = [ 520887429 0 0; 25209327 3.4910864 529.6909651; 610600 3.841153656 1059.38193; 282029 2.574198799 632.7837393; 187647 2.075903801 522.5774181; 86793 0.710010906 419.4846439; 72063 0.214656947 536.8045121; 65517 5.979958508 316.3918697; 30135 2.161320584 949.175609; 29135 1.677592437 103.0927742; 23947 0.274578549 7.113547001; 23453 3.540231473 735.8765135; 22284 4.193627735 1589.072895; 13033 2.960430557 1162.474704; 12749 2.715501029 1052.268383; 9703 1.906695724 206.1855484; 9161 4.413526189 213.2990954; 7895 2.479075514 426.5981909; 7058 2.181847531 1265.567479; 6138 6.264175425 846.0828348; 5477 5.657293252 639.8972863; 4170 2.016050339 515.4638711; 4137 2.722199797 625.6701923; 3503 0.565312974 1066.495477; 2617 2.009939671 1581.959348; 2500 4.551820559 838.9692878; 2128 6.127514618 742.9900605; 1912 0.856219274 412.3710969; 1611 3.088677893 1368.660253; 1479 2.680261914 1478.866574; 1231 1.890429797 323.5054167; 1217 1.80171561 110.2063212; 1015 1.386732377 454.9093665; 999 2.872089401 309.2783227; 961 4.548769898 2118.76386; 886 4.147859485 533.6231184; 821 1.593425344 1898.351218; 812 5.940918991 909.8187331; 777 3.676969547 728.7629665; 727 3.988246864 1155.361157; 655 2.790656042 1685.052123; 654 3.381507753 1692.16567; 621 4.82284339 956.289156; 615 2.276249156 942.062062; 562 0.080959872 543.9180591; 542 0.283602664 525.7588118 ]; r1 = [ 1271802 2.649375111 529.6909651; 61662 3.00076251 1059.38193; 53444 3.897176442 522.5774181; 41390 0 0; 31185 4.882766635 536.8045121; 11847 2.413295882 419.4846439; 9166 4.759794086 7.113547001; 3404 3.34688538 1589.072895; 3203 5.210832855 735.8765135; 3176 2.792979871 103.0927742; 2806 3.742236936 515.4638711; 2677 4.330528787 1052.268383; 2600 3.634351016 206.1855484; 2412 1.469473083 426.5981909; 2101 3.927626823 639.8972863; 1646 5.309535109 1066.495477; 1641 4.416286698 625.6701923; 1050 3.16113623 213.2990954; 1025 2.55432643 412.3710969; 806 2.677508014 632.7837393; 741 2.170946306 1162.474704; 677 6.249534798 838.9692878; 567 4.576554147 742.9900605; 485 2.468827932 949.175609; 469 4.709734635 543.9180591; 445 0.402811814 323.5054167; 416 5.368360182 728.7629665; 402 4.605288415 309.2783227; 347 4.681488087 14.227094; 338 3.167819511 956.289156; 261 5.342903061 846.0828348; 247 3.923138235 942.062062; 220 4.84210965 1368.660253; 203 5.599954254 1155.361157; 200 4.438888144 1045.154836; 197 3.705514614 2118.76386; 196 3.758775871 199.0720014; 184 4.265267697 95.97922722; 180 4.401654912 532.8723588; 170 4.846474889 526.5095714; 146 6.129583655 533.6231184; 133 1.322457359 110.2063212; 132 4.511879508 525.7588118 ]; r2 = [ 79645 1.358658966 529.6909651; 8252 5.777739354 522.5774181; 7030 3.274769658 536.8045121; 5314 1.838351097 1059.38193; 1861 2.976821394 7.113547001; 964 5.48031822 515.4638711; 836 4.198898817 419.4846439; 498 3.141592654 0; 427 2.227531018 639.8972863; 406 3.782507304 1066.495477; 377 2.242483529 1589.072895; 363 5.367618473 206.1855484; 342 6.099229693 1052.268383; 339 6.12690864 625.6701923; 333 0.003289612 426.5981909; 280 4.261625558 412.3710969; 257 0.96295365 632.7837393; 230 0.705307662 735.8765135; 201 3.068506234 543.9180591; 200 4.428841653 103.0927742; 139 2.932356716 14.227094; 114 0.787139113 728.7629665; 95 1.704980411 838.9692878; 86 5.14434752 323.5054167; 83 0.058348735 309.2783227; 80 2.981223618 742.9900605; 75 1.604951959 956.289156; 70 1.509883575 213.2990954; 67 5.473071781 199.0720014; 62 6.101378899 1045.154836; 56 0.955348105 1162.474704; 52 5.584356256 942.062062; 50 2.720631623 532.8723588; 45 5.524456214 508.3503241; 44 0.271181526 526.5095714; 40 5.945665062 95.97922722 ]; r3 = [ 3519 6.058006338 529.6909651; 1073 1.673213458 536.8045121; 916 1.413296761 522.5774181; 342 0.522965427 1059.38193; 255 1.196254735 7.113547001; 222 0.952252262 515.4638711; 90 3.141592654 0; 69 2.268852823 1066.495477; 58 1.413897453 543.9180591; 58 0.525801176 639.8972863; 51 5.980163647 412.3710969; 47 1.57864238 625.6701923; 43 6.116896091 419.4846439; 37 1.182627623 14.227094; 34 1.66671707 1052.268383; 34 0.847849779 206.1855484; 31 1.042902459 1589.072895; 30 4.63236245 426.5981909; 21 2.500712438 728.7629665; 15 0.891369984 199.0720014; 14 0.960401971 508.3503241; 13 1.502337886 1045.154836; 12 2.609526145 735.8765135; 12 3.555135101 323.5054167; 11 1.790414376 309.2783227; 11 6.278451127 956.289156; 10 6.260168595 103.0927742; 9 3.451268125 838.9692878 ]; r4 = [ 129 0.084193096 536.8045121; 113 4.248588558 529.6909651; 83 3.297549094 522.5774181; 38 2.733266111 515.4638711; 27 5.691425886 7.113547001; 18 5.400125369 1059.38193; 13 6.015604161 543.9180591; 9 0.768139465 1066.495477; 8 5.682280657 14.227094; 7 1.427512921 412.3710969; 6 5.122869325 639.8972863; 5 3.335019473 625.6701923; 3 3.403348051 1052.268383; 3 4.16090413 728.7629665; 3 2.898020351 426.5981909 ]; r5 = [ 11 4.75 536.80; 4 5.92 522.58; 2 5.57 515.46; 2 4.30 543.92; 2 3.69 7.11; 2 4.13 1059.38; 2 5.49 1066.50 ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); sr5 = zzterm(r5, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4 + sr5*t5)/100000000; return [lambda, beta, delta]; endfunction function hmars(day) ## returns the heliocentric ecliptic latitude of Mars given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## More terms than Meeus included as an experiment .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC l0 = [ 620347711.583000000000 0.000000000000 0.000000000000; 18656368.100000000000 5.050371003030 3340.612426699800; 1108216.792000000000 5.400998369580 6681.224853399600; 91798.394000000000 5.754787451110 10021.837280099400; 27744.987000000000 5.970495129420 3.523118349000; 12315.897000000000 0.849560812380 2810.921461605200; 10610.230000000000 2.939585249730 2281.230496510600; 8926.772000000000 4.156978459390 0.017253652200; 8715.688000000000 6.110051597920 13362.449706799200; 7774.867000000000 3.339686550740 5621.842923210400; 6797.552000000000 0.364622436260 398.149003408200; 4161.101000000000 0.228149753300 2942.463423291600; 3575.079000000000 1.661865401410 2544.314419883400; 3075.250000000000 0.856965970820 191.448266111600; 2937.543000000000 6.078937114080 0.067310302800; 2628.122000000000 0.648061435700 3337.089308350800; 2579.842000000000 0.029967061970 3344.135545048800; 2389.420000000000 5.038964013490 796.298006816400; 1798.808000000000 0.656340268440 529.690965094600; 1546.408000000000 2.915796333920 1751.539531416000; 1528.140000000000 1.149793062280 6151.533888305000; 1286.232000000000 3.067959246260 2146.165416475200; 1264.356000000000 3.622750922310 5092.151958115800; 1024.907000000000 3.693342935550 8962.455349910200; 891.567000000000 0.182938990900 16703.062133499000; 858.760000000000 2.400937042040 2914.014235823800; 832.724000000000 4.494957534580 3340.629680352000; 832.718000000000 2.464185912820 3340.595173047600; 748.724000000000 3.822483994680 155.420399434200; 723.863000000000 0.674975658010 3738.761430108000; 712.899000000000 3.663360147880 1059.381930189200; 655.163000000000 0.488640751760 3127.313331261800; 635.557000000000 2.921827042750 8432.764384815600; 552.746000000000 4.474788630160 1748.016413067000; 550.472000000000 3.810012054080 0.980321068200; 472.164000000000 3.625478194100 1194.447010224600; 425.972000000000 0.553651381720 6283.075849991400; 415.132000000000 0.496623147740 213.299095438000; 312.141000000000 0.998533228430 6677.701735050600; 306.552000000000 0.380528629730 6684.747971748600; 302.377000000000 4.486181503210 3532.060692811400; 299.396000000000 2.783237056970 6254.626662523600; 293.199000000000 4.221312779140 20.775395492400; 283.600000000000 5.768854941230 3149.164160588200; 281.073000000000 5.881633729450 1349.867409658800; 274.035000000000 0.133725012110 3340.679737002600; 274.028000000000 0.542221418410 3340.545116397000; 238.857000000000 5.371554716720 4136.910433516200; 236.114000000000 5.755045155760 3333.498879699000; 231.185000000000 1.282406852940 3870.303391794400; 221.225000000000 3.504666722030 382.896532223200; 204.161000000000 2.821332661850 1221.848566321400; 193.126000000000 3.357151377450 3.590428651800; 188.639000000000 1.491030164860 9492.146315004800; 179.196000000000 1.005611125740 951.718406250600; 174.068000000000 2.413603325760 553.569402842400; 172.110000000000 0.439430417190 5486.777843175000; 160.011000000000 3.948547351920 4562.460993021200; 144.305000000000 1.418741934180 135.065080035400; 139.897000000000 3.325925161640 2700.715140385800; 138.245000000000 4.301451769150 7.113547000800; 130.993000000000 4.044917202640 12303.067776610000; 128.102000000000 2.208066510080 1592.596013632800; 128.062000000000 1.806656433320 5088.628839766800; 116.945000000000 3.128052822070 7903.073419721000; 113.486000000000 3.700707981230 1589.072895283800; 110.375000000000 1.051950796870 242.728603974000; 104.541000000000 0.785353820760 8827.390269874800; 100.090000000000 3.243437408610 11773.376811515400; 98.947000000000 4.845582947400 6681.242107051800; 98.946000000000 2.814811403710 6681.207599747400; 95.592000000000 0.539541811490 20043.674560198800; 86.931000000000 2.201867405230 11243.685846420800; 86.751000000000 1.020922215630 7079.373856807800; 84.187000000000 3.989707207300 4399.994356889000; 83.749000000000 3.202561309900 4690.479836358600; 75.034000000000 0.766434182520 6467.925757961600; 73.476000000000 2.184280125670 8429.241266466600; 72.091000000000 5.846721025250 5884.926846583200; 71.437000000000 2.803075500160 3185.192027265600; 68.984000000000 3.763997317880 6041.327567085600; 68.414000000000 2.738349144120 2288.344043511400; 66.706000000000 0.736306207660 3723.508958923000; 65.320000000000 2.681185975780 28.449187467800; 63.376000000000 0.912962407980 3553.911522137800; 63.314000000000 4.527714704700 426.598190876000; 61.683000000000 6.168315094190 2274.116949509800; 56.629000000000 5.062504102060 15.252471185000; 56.396000000000 1.687271503040 6872.673119511200; 55.909000000000 3.462608334950 263.083923372800; 55.488000000000 4.606254670200 4292.330832950400; 52.256000000000 0.899415313070 9623.688276691200; 51.678000000000 2.813074926820 3339.632105631600; 51.332000000000 4.148236365340 3341.592747768000; 48.542000000000 3.956704187190 4535.059436924400; 45.905000000000 0.287189814970 5614.729376209600; 45.829000000000 0.787842350620 1990.745017041000; 44.174000000000 3.195297367020 5628.956470211200; 41.939000000000 3.583264251150 8031.092263058400; 41.223000000000 6.020193299220 3894.181829542200; 40.671000000000 3.138326218290 9595.239089223400; 39.495000000000 5.632253921600 3097.883822725790; 38.790000000000 1.351984987950 10018.314161750400; 38.352000000000 5.828807074260 3191.049229565200; 38.206000000000 2.348359840630 162.466636132200; 38.107000000000 0.734019463200 10025.360398448400; 37.752000000000 4.154829552990 2803.807914604400; 37.135000000000 0.685081507740 2818.035008606000; 36.716000000000 2.637207751020 692.157601226800; 34.031000000000 2.595440825090 11769.853693166400; 33.626000000000 6.119924010520 6489.776587288000; 33.148000000000 1.140237700040 5.522924307400; 32.562000000000 0.484006593330 6681.292163702400; 32.561000000000 0.892503168880 6681.157543096800; 31.168000000000 3.981609129820 20.355319398800; 29.007000000000 2.427073856740 3319.837031207400; 28.686000000000 5.720554567340 7477.522860216000; 27.584000000000 1.596912030580 7210.915818494200; 27.540000000000 6.083899423370 6674.111306398800; 27.278000000000 4.556453281220 3361.387822192200; 26.357000000000 1.345326465740 3496.032826134000; 25.637000000000 0.249635234200 522.577418093800; 25.512000000000 3.432423528040 3443.705200918400; 25.380000000000 0.520931161120 10.636665349800; 24.554000000000 4.003231830880 11371.704689758200; 24.378000000000 0.969946964130 632.783739313200; 23.764000000000 1.840583772560 12832.758741704600; 23.079000000000 4.749902142230 3347.725973700600; 22.816000000000 3.526282121060 1648.446757197400; 22.662000000000 3.954463244170 4989.059183897200; 22.542000000000 5.648617034380 2388.894020449200; 22.274000000000 0.721061337210 266.607041721800; 21.530000000000 6.153887571770 3264.346355424200; 21.343000000000 4.282187578630 4032.770027926600; 21.202000000000 3.118244722840 2957.715894476600; 20.158000000000 3.671315049460 1758.653078416800; 20.093000000000 1.082474160650 7064.121385622800; 19.849000000000 2.376689207450 10713.994881326200; 17.709000000000 3.697423439740 3344.202855351600 ]; l1 = [ 334085627474.34200000000 0.00000000000 0.00000000000; 1458227.05100000000 3.60426053609 3340.61242669980; 164901.34300000000 3.92631250962 6681.22485339960; 19963.33800000000 4.26594061030 10021.83728009940; 3452.39900000000 4.73210386365 3.52311834900; 2485.48000000000 4.61277567318 13362.44970679920; 841.55100000000 4.45858256765 2281.23049651060; 537.56600000000 5.01589727492 398.14900340820; 521.04100000000 4.99422678175 3344.13554504880; 432.61400000000 2.56066402860 191.44826611160; 429.65600000000 5.31646162367 155.42039943420; 381.74700000000 3.53881289437 796.29800681640; 314.12900000000 4.96335266049 16703.06213349900; 282.80400000000 3.15967518204 2544.31441988340; 205.66400000000 4.56891455660 2146.16541647520; 168.80500000000 1.32894813366 3337.08930835080; 157.58700000000 4.18501035954 1751.53953141600; 133.68600000000 2.23325104196 0.98032106820; 133.56300000000 5.97421903927 1748.01641306700; 117.59100000000 6.02407213861 6151.53388830500; 116.56100000000 2.21347652545 1059.38193018920; 113.87600000000 2.12869455089 1194.44701022460; 113.59500000000 5.42803224317 3738.76143010800; 91.09800000000 1.09627836591 1349.86740965880; 85.34200000000 3.90854841008 553.56940284240; 83.30100000000 5.29636626272 6684.74797174860; 80.77600000000 4.42813405865 529.69096509460; 79.53100000000 2.24864266330 8962.45534991020; 72.94600000000 2.50189460554 951.71840625060; 72.50500000000 5.84208163240 242.72860397400; 71.48700000000 3.85636094435 2914.01423582380; 67.58200000000 5.02327686473 382.89653222320; 65.08900000000 1.01802439311 3340.59517304760; 65.08900000000 3.04879603978 3340.62968035200; 61.50800000000 4.15183159800 3149.16416058820; 56.52000000000 3.88813699320 4136.91043351620; 48.47700000000 4.87362121538 213.29909543800; 47.61300000000 1.18238046057 3333.49887969900; 46.58400000000 1.31452419914 3185.19202726560; 41.34300000000 0.71385375517 1592.59601363280; 40.27200000000 2.72542480614 7.11354700080; 40.05500000000 5.31611875491 20043.67456019880; 32.88600000000 5.41067411968 6283.07584999140; 28.24400000000 0.04534124888 9492.14631500480; 26.57900000000 3.88960724782 1221.84856632140; 26.55400000000 5.11271747607 2700.71514038580; 23.33500000000 6.16762213077 3532.06069281140; 22.79700000000 1.54504711003 2274.11694950980; 22.61200000000 0.83775884934 3097.88382272579; 22.43100000000 5.46592525433 20.35531939880; 22.29400000000 5.88516997273 3870.30339179440; 21.42500000000 4.97081508139 3340.67973700260; 21.41800000000 5.37934044204 3340.54511639700; 21.10400000000 3.52525428062 15.25247118500; 20.43100000000 2.36353950189 1589.07289528380; 20.18600000000 3.36375535766 5088.62883976680; 20.02900000000 4.73119428749 4690.47983635860; 19.96400000000 5.78652958398 7079.37385680780; 19.67500000000 2.57805423988 12303.06777661000; 19.46800000000 0.49216434489 6677.70173505060; 19.45500000000 2.53112676345 4399.99435688900; 18.50500000000 5.57863503922 1990.74501704100; 17.81100000000 6.12537931996 4292.33083295040; 16.59900000000 1.25519718278 3894.18182954220; 16.47200000000 2.60291845066 3341.59274776800; 15.38100000000 2.47009470350 4535.05943692440; 15.30700000000 2.26515985343 3723.50895892300; 15.00000000000 1.03464802434 2288.34404351140; 14.70500000000 3.36979890389 6681.24210705180; 14.70500000000 1.33902715586 6681.20759974740; 13.64400000000 1.97739547259 5614.72937620960; 13.53500000000 2.12334410410 5486.77784317500; 13.27500000000 3.42243595774 5621.84292321040; 13.01300000000 1.51424752315 5628.95647021120; 12.95000000000 5.61929676688 10025.36039844840; 12.68200000000 2.95022113262 3496.03282613400; 11.85400000000 5.47630686910 3553.91152213780; 11.85000000000 3.12701832949 426.59819087600; 11.76400000000 2.58551521265 8432.76438481560; 11.35300000000 6.23438193885 135.06508003540; 11.13200000000 5.84178807242 2803.80791460440; 10.95800000000 4.15771850822 2388.89402044920; 10.86700000000 5.28184140482 2818.03500860600; 10.47200000000 2.73581537999 2787.04302385740; 9.70800000000 4.52957217749 6489.77658728800; 8.84000000000 4.23294294197 7477.52286021600; 8.69500000000 4.43542512603 5092.15195811580; 8.65000000000 4.33256981793 3339.63210563160; 8.56200000000 3.16141186861 162.46663613220; 8.49000000000 1.91378007528 11773.37681151540; 8.43400000000 3.16292250873 3347.72597370060; 8.34400000000 2.18273563186 23.87843774780; 8.13300000000 1.61295625304 2957.71589447660; 8.03400000000 5.69983564288 6041.32756708560; 7.69600000000 5.71877332978 9623.68827669120; 7.37200000000 6.17831593269 3583.34103067380; 6.66400000000 5.07517838003 8031.09226305840 ]; l2 = [ 58015.79100000000 2.04979463279 3340.61242669980; 54187.64500000000 0.00000000000 0.00000000000; 13908.42600000000 2.45742359888 6681.22485339960; 2465.10400000000 2.80000020929 10021.83728009940; 398.37900000000 3.14118428289 13362.44970679920; 222.02200000000 3.19436080019 3.52311834900; 120.95700000000 0.54325292454 155.42039943420; 61.51700000000 3.48529427371 16703.06213349900; 53.63800000000 3.54191121461 3344.13554504880; 34.26800000000 6.00188499119 2281.23049651060; 31.66500000000 4.14015171788 191.44826611160; 29.83900000000 1.99870679845 796.29800681640; 23.16800000000 4.33403365928 242.72860397400; 21.65900000000 3.44532466378 398.14900340820; 20.37000000000 5.42191375400 553.56940284240; 16.22700000000 0.65678953303 0.98032106820; 16.04400000000 6.11000472441 2146.16541647520; 15.64800000000 1.22086121940 1748.01641306700; 14.92700000000 6.09541783564 3185.19202726560; 14.41600000000 4.01923812101 951.71840625060; 14.31700000000 2.61851897591 1349.86740965880; 13.35200000000 0.60189008414 1194.44701022460; 11.93400000000 3.86122163021 6684.74797174860; 11.26000000000 4.71822363671 2544.31441988340; 10.39600000000 0.25038714677 382.89653222320; 9.46800000000 0.68170713564 1059.38193018920; 9.22900000000 3.83209092321 20043.67456019880; 9.00500000000 3.88271826102 3738.76143010800; 7.50100000000 5.46498630412 1751.53953141600; 6.85900000000 2.57522504136 3149.16416058820; 6.68100000000 2.37843690339 4136.91043351620; 6.49700000000 5.47773072872 1592.59601363280; 6.31100000000 2.34104793674 3097.88382272579; 5.87000000000 1.14783576679 7.11354700080; 4.76400000000 2.89684755585 3333.49887969900; 4.64700000000 4.42957708526 6151.53388830500; 4.16600000000 3.68631477611 5614.72937620960; 4.04500000000 6.12493402657 5628.95647021120; 3.65300000000 4.06679068397 1990.74501704100; 3.61800000000 2.46868561769 529.69096509460; 3.27700000000 0.68101740787 8962.45534991020; 3.25300000000 2.79565340390 3894.18182954220; 3.09100000000 4.56861203364 3496.03282613400; 2.92100000000 5.41458945995 2914.01423582380; 2.92100000000 1.23050883841 2787.04302385740; 2.88800000000 3.41062353663 3337.08930835080; 2.78400000000 1.38911141844 4292.33083295040; 2.63000000000 4.67640929857 3583.34103067380; 2.62000000000 1.04061894134 3341.59274776800; 2.60200000000 2.64911714813 2388.89402044920; 2.59400000000 1.49510566123 3340.62968035200; 2.59300000000 5.74934234498 3340.59517304760; 2.41800000000 0.96341462666 4535.05943692440; 2.41600000000 1.04749173375 4399.99435688900; 2.38600000000 4.27072575550 7079.37385680780; 2.35700000000 4.84628239765 9492.14631500480; 2.34400000000 4.18104725028 10025.36039844840; 2.34400000000 0.01425578204 4690.47983635860; 2.19100000000 3.26449527357 213.29909543800; 2.18700000000 0.16036551231 6525.80445396540; 2.12600000000 0.48290227706 2700.71514038580; 1.83000000000 0.97181050149 1589.07289528380; 1.76300000000 2.51660121293 2810.92146160520; 1.75900000000 3.81170701376 3723.50895892300; 1.63300000000 1.10703599922 12303.06777661000; 1.62900000000 4.94267977718 1221.84856632140; 1.61700000000 4.95614491689 5088.62883976680; 1.51300000000 2.92614134711 640.87760738220; 1.50400000000 0.11031912519 2957.71589447660; 1.43100000000 2.97747408389 6489.77658728800; 1.40800000000 1.54395721611 3347.72597370060; 1.40100000000 3.85907867678 6283.07584999140; 1.39200000000 2.73498041122 7477.52286021600; 1.33800000000 5.29685392418 6677.70173505060; 1.23600000000 3.77245965590 2699.73481931760; 1.23400000000 1.88931735265 6681.24210705180; 1.23400000000 6.14168429036 6681.20759974740 ]; l3 = [ 1482.42300000000 0.44434694876 3340.61242669980; 662.09500000000 0.88469178686 6681.22485339960; 188.26800000000 1.28799982497 10021.83728009940; 41.47400000000 1.64850786997 13362.44970679920; 25.99400000000 0.00000000000 0.00000000000; 22.66100000000 2.05267665262 155.42039943420; 10.45400000000 1.58006906385 3.52311834900; 8.02400000000 1.99858757687 16703.06213349900; 4.90000000000 2.82452457966 242.72860397400; 3.78200000000 2.01914272515 3344.13554504880; 3.17600000000 4.59144897927 3185.19202726560; 3.13400000000 0.65044714325 553.56940284240; 1.68400000000 5.53835848782 951.71840625060; 1.51100000000 5.71795850828 191.44826611160; 1.44800000000 0.45869142895 796.29800681640; 1.44200000000 2.34368495577 20043.67456019880; 1.30200000000 5.36284013048 0.98032106820; 1.16900000000 4.14601161433 1349.86740965880; 1.13300000000 2.38180830662 6684.74797174860; 1.03700000000 1.76892750558 382.89653222320; 0.89400000000 5.33688328934 1194.44701022460; 0.80700000000 2.74798886181 1748.01641306700; 0.64700000000 3.17645475605 3583.34103067380; 0.64000000000 6.10665147849 3496.03282613400; 0.56700000000 5.85922384979 7.11354700080; 0.55800000000 1.85212342360 398.14900340820; 0.51900000000 4.93376176788 6525.80445396540; 0.50800000000 1.01139298015 3149.16416058820; 0.45200000000 5.98109989317 2787.04302385740; 0.40500000000 1.27295444059 2281.23049651060 ]; l4 = [ 113.96900000000 3.14159265359 0.00000000000; 28.72500000000 5.63662412043 6681.22485339960; 24.44700000000 5.13868481454 3340.61242669980; 11.18700000000 6.03161074431 10021.83728009940; 3.25200000000 0.13228350651 13362.44970679920; 3.19000000000 3.56267988299 155.42039943420; 0.78700000000 0.49340783377 16703.06213349900; 0.77600000000 1.31734531594 242.72860397400; 0.49400000000 3.06356214498 3185.19202726560; 0.37400000000 2.15785846355 553.56940284240; 0.33100000000 6.23159792887 3.52311834900; 0.19700000000 0.44350153983 3344.13554504880; 0.18100000000 0.81531283571 20043.67456019880; 0.16800000000 3.73509781785 3496.03282613400; 0.11500000000 1.66898531261 3583.34103067380; 0.09200000000 3.40530361815 6525.80445396540; 0.08600000000 0.79259553758 6684.74797174860; 0.06400000000 4.47443580658 2787.04302385740; 0.04500000000 5.17216217058 3097.88382272579; 0.04100000000 1.21875027733 23384.28698689860; 0.03600000000 5.53975653407 3149.16416058820 ]; l5 = [ 0.86800000000 3.14159265359 0.00000000000; 0.71000000000 4.04089996521 6681.22485339960; 0.51000000000 4.49214901625 10021.83728009940; 0.35700000000 5.07435505061 155.42039943420; 0.22300000000 3.51351884241 3340.61242669980 ]; .. latitude series b0 = [ 3197134.98600000000 3.76832042432 3340.61242669980; 298033.23400000000 4.10616996243 6681.22485339960; 289104.74200000000 0.00000000000 0.00000000000; 31365.53800000000 4.44651052853 10021.83728009940; 3484.10000000000 4.78812547889 13362.44970679920; 443.40100000000 5.02642620491 3344.13554504880; 442.99900000000 5.65233015876 3337.08930835080; 399.10900000000 5.13056814700 16703.06213349900; 292.50600000000 3.79290644595 2281.23049651060; 181.98200000000 6.13648011704 6151.53388830500; 163.15900000000 4.26399626634 529.69096509460; 159.67800000000 2.23194610246 1059.38193018920; 149.29700000000 2.16501209917 5621.84292321040; 142.68600000000 1.18215016110 3340.59517304760; 142.68500000000 3.21292180820 3340.62968035200; 139.32300000000 2.41796344238 8962.45534991020; 86.37700000000 5.74429648412 3738.76143010800; 83.27600000000 5.98866315739 6677.70173505060; 82.54400000000 5.36667872319 6684.74797174860; 73.64000000000 5.09187524843 398.14900340820; 72.66000000000 5.53775710437 6283.07584999140; 63.11100000000 0.73049113369 5884.92684658320; 62.33800000000 4.85071999184 2942.46342329160; 60.11600000000 3.67960808826 796.29800681640; 47.19900000000 4.52184736343 3149.16416058820; 46.95300000000 5.13486627234 3340.67973700260; 46.95100000000 5.54339723804 3340.54511639700; 46.63000000000 5.47361665459 20043.67456019880; 45.58800000000 2.13262507507 2810.92146160520; 41.26900000000 0.20003189001 9492.14631500480; 38.54000000000 4.08008443274 4136.91043351620; 33.06900000000 4.06581918329 1751.53953141600; 29.69400000000 5.92218297386 3532.06069281140 ]; b1 = [ 350068.84500000000 5.36847836211 3340.61242669980; 14116.03000000000 3.14159265359 0.00000000000; 9670.75500000000 5.47877786506 6681.22485339960; 1471.91800000000 3.20205766795 10021.83728009940; 425.86400000000 3.40843812875 13362.44970679920; 102.03900000000 0.77617286189 3337.08930835080; 78.84800000000 3.71768293865 16703.06213349900; 32.70800000000 3.45803723682 5621.84292321040; 26.17100000000 2.48293558065 2281.23049651060; 20.71200000000 1.44120802297 6151.53388830500; 18.29400000000 6.03102943125 529.69096509460; 16.97500000000 4.81115186866 3344.13554504880; 15.68000000000 3.93075566599 8962.45534991020; 15.62200000000 2.78241383265 3340.59517304760; 15.62200000000 4.81318636318 3340.62968035200; 14.26800000000 0.24640247719 2942.46342329160; 13.77100000000 1.67983063909 3532.06069281140; 13.06700000000 0.97324736181 6677.70173505060; 12.71100000000 4.04546734935 20043.67456019880; 12.49300000000 2.25620513522 5884.92684658320; 8.80000000000 0.34079528233 398.14900340820 ]; b2 = [ 16726.69000000000 0.60221392419 3340.61242669980; 4986.79900000000 3.14159265359 0.00000000000; 302.14100000000 5.55871276021 6681.22485339960; 25.76700000000 1.89662673499 13362.44970679920; 21.45200000000 0.91749968618 10021.83728009940; 11.82000000000 2.24240738700 3337.08930835080; 7.98500000000 2.24892866611 16703.06213349900; 2.96000000000 5.89425825808 3496.03282613400; 2.44500000000 5.18770525274 5621.84292321040; 1.77900000000 2.58759968520 20043.67456019880; 1.50100000000 3.18533003542 3532.06069281140; 1.42800000000 1.25238140580 2281.23049651060; 1.25900000000 4.80695172904 3185.19202726560; 1.10900000000 3.80982317372 5884.92684658320; 1.02900000000 2.35029907056 6677.70173505060 ]; b3 = [ 606.50600000000 1.98050633529 3340.61242669980; 42.61100000000 0.00000000000 0.00000000000; 13.65200000000 1.79588228800 6681.22485339960; 2.73000000000 3.45377082121 10021.83728009940; 0.92900000000 3.75226159072 3337.08930835080; 0.60700000000 0.10618486408 13362.44970679920; 0.61700000000 1.14471772765 3496.03282613400; 0.47900000000 0.70504966293 16703.06213349900; 0.18500000000 3.28778562029 3185.19202726560; 0.16900000000 0.29980532608 5621.84292321040 ]; b4 = [ 11.33400000000 3.45724352586 3340.61242669980; 13.36900000000 0.00000000000 0.00000000000; 0.74400000000 0.50445805257 6681.22485339960; 0.14800000000 1.05056602649 10021.83728009940; 0.10200000000 2.66185835593 3496.03282613400; 0.05300000000 5.27888218929 3337.08930835080 ]; b5 = [ 0.457 4.86794125358 3340.61242669980; 0.053 5.30547050586 6681.22485339960 ]; .. radius vector r0 = [ 153033488.27600000000 0.00000000000 0.00000000000; 14184953.15300000000 3.47971283519 3340.61242669980; 660776.35700000000 3.81783442097 6681.22485339960; 46179.11700000000 4.15595316284 10021.83728009940; 8109.73800000000 5.55958460165 2810.92146160520; 7485.31500000000 1.77238998069 5621.84292321040; 5523.19300000000 1.36436318880 2281.23049651060; 3825.16000000000 4.49407182408 13362.44970679920; 2484.38500000000 4.92545577893 2942.46342329160; 2306.53900000000 0.09081742493 2544.31441988340; 1999.39900000000 5.36059605227 3337.08930835080; 1960.19800000000 4.74249386323 3344.13554504880; 1167.11500000000 2.11261501155 5092.15195811580; 1102.82800000000 5.00908264160 398.14900340820; 992.25200000000 5.83862401067 6151.53388830500; 899.07700000000 4.40790433994 529.69096509460; 807.34800000000 2.10216647104 1059.38193018920; 797.91000000000 3.44839026172 796.29800681640; 740.98000000000 1.49906336892 2146.16541647520; 725.58300000000 1.24516913473 8432.76438481560; 692.34000000000 2.13378814785 8962.45534991020; 633.14400000000 0.89353285018 3340.59517304760; 633.14000000000 2.92430448169 3340.62968035200; 629.97600000000 1.28738135858 1751.53953141600; 574.35200000000 0.82896196337 2914.01423582380; 526.18700000000 5.38292276228 3738.76143010800; 472.77600000000 5.19850457873 3127.31333126180; 348.09500000000 4.83219198908 16703.06213349900; 283.70200000000 2.90692294913 3532.06069281140; 279.55200000000 5.25749247548 6283.07584999140; 275.50100000000 1.21767967781 6254.62666252360; 275.22400000000 2.90818883832 1748.01641306700; 269.89100000000 3.76394728622 5884.92684658320; 239.13300000000 2.03669896238 1194.44701022460; 233.82700000000 5.10546492529 5486.77784317500; 228.12800000000 3.25529020620 6872.67311951120; 223.19000000000 4.19861593779 3149.16416058820; 219.42800000000 5.58340248784 191.44826611160; 208.33600000000 4.84626442122 3340.67973700260; 208.33300000000 5.25476080773 3340.54511639700; 186.21300000000 5.69871555748 6677.70173505060; 182.68600000000 5.08062683355 6684.74797174860; 178.61300000000 4.18423025538 3333.49887969900; 175.99500000000 5.95341786369 3870.30339179440; 163.53400000000 3.79889068111 4136.91043351620; 144.28600000000 0.21296012258 5088.62883976680; 141.75900000000 2.47790321309 4562.46099302120; 133.12000000000 1.53910106710 7903.07341972100; 128.55500000000 5.49883294915 8827.39026987480; 118.78100000000 2.12178071222 1589.07289528380; 114.94100000000 4.31745088059 1349.86740965880; 111.53800000000 0.55339169625 11243.68584642080; 102.09600000000 6.18138550087 9492.14631500480; 86.65900000000 1.74988330093 2700.71514038580; 85.31200000000 1.61621097912 4690.47983635860; 84.47000000000 0.62274593110 1592.59601363280; 83.21200000000 0.61553380568 8429.24126646660; 82.49800000000 1.62227044590 11773.37681151540; 71.82600000000 2.47489899385 12303.06777661000; 68.59900000000 2.40197828418 4399.99435688900; 66.50900000000 2.21307705185 6041.32756708560; 63.64100000000 2.67334126661 426.59819087600; 62.01500000000 1.10065866221 1221.84856632140; 58.95900000000 3.26242666052 6681.24210705180; 58.95900000000 1.23165502899 6681.20759974740; 58.55900000000 4.72052787516 213.29909543800; 55.81100000000 1.23288325946 3185.19202726560; 55.68600000000 5.44686699242 3723.50895892300; 54.98900000000 5.72691385306 951.71840625060; 52.41800000000 3.02366828926 4292.33083295040; 51.56100000000 5.72326937712 7079.37385680780; 48.93900000000 5.61614696751 3553.91152213780; 45.41400000000 5.43290921705 6467.92575796160; 44.62900000000 2.01473640390 8031.09226305840; 44.29200000000 5.00341366850 5614.72937620960; 43.25600000000 1.03732072925 11769.85369316640; 42.44400000000 2.26551590902 155.42039943420; 42.19100000000 1.63253742760 5628.95647021120; 39.23700000000 1.24237122859 3339.63210563160; 38.95600000000 2.57760416009 3341.59274776800; 36.43500000000 4.43921812388 3894.18182954220; 35.98000000000 1.15966567007 2288.34404351140; 35.26500000000 5.49029710802 1990.74501704100; 33.62300000000 5.17029029766 20043.67456019880; 33.06500000000 0.85467740581 553.56940284240; 32.25900000000 2.38215172582 4535.05943692440; 31.97200000000 1.93970478412 382.89653222320; 31.94300000000 4.59258406791 2274.11694950980; 31.87000000000 4.37521442752 3.52311834900; 30.34500000000 2.44177670130 11371.70468975820; 29.35000000000 4.06034813442 3097.88382272579; 27.90400000000 4.25805969214 3191.04922956520; 27.54300000000 1.57668567401 9595.23908922340; 26.16600000000 5.58466944895 9623.68827669120; 25.15900000000 0.81355213242 10713.99488132620; 24.77200000000 5.38970742761 2818.03500860600; 24.73200000000 2.58034797703 2803.80791460440; 23.35900000000 6.01453778225 3496.03282613400; 22.07000000000 0.85747723964 3319.83703120740 ]; r1 = [ 1107433.34000000000 2.03250524950 3340.61242669980; 103175.88600000000 2.37071845682 6681.22485339960; 12877.20000000000 0.00000000000 0.00000000000; 10815.88000000000 2.70888093803 10021.83728009940; 1194.55000000000 3.04702182503 13362.44970679920; 438.57900000000 2.88835072628 2281.23049651060; 395.69800000000 3.42324611291 3344.13554504880; 182.57200000000 1.58428644001 2544.31441988340; 135.85000000000 3.38507017993 16703.06213349900; 128.36200000000 6.04343360441 3337.08930835080; 128.20400000000 0.62991220570 1059.38193018920; 127.06800000000 1.95389775740 796.29800681640; 118.44300000000 2.99761345074 2146.16541647520; 87.53700000000 3.42052758979 398.14900340820; 83.02600000000 3.85574986653 3738.76143010800; 75.59800000000 4.45101839349 6151.53388830500; 71.99900000000 2.76442180680 529.69096509460; 66.54200000000 2.54892602695 1751.53953141600; 66.43000000000 4.40597549957 1748.01641306700; 57.51800000000 0.54354327916 1194.44701022460; 54.31400000000 0.67750943459 8962.45534991020; 51.03500000000 3.72585409207 6684.74797174860; 49.42800000000 5.72959428364 3340.59517304760; 49.42400000000 1.47717922226 3340.62968035200; 48.31800000000 2.58061691301 3149.16416058820; 47.86300000000 2.28527896843 2914.01423582380; 38.95300000000 2.31900090554 4136.91043351620; 37.17600000000 5.81439911546 1349.86740965880; 36.38400000000 6.02728752344 3185.19202726560; 36.03600000000 5.89508336048 3333.49887969900; 31.11500000000 0.97832506960 191.44826611160; 27.24400000000 5.41367977087 1592.59601363280; 24.30000000000 3.75843924498 155.42039943420; 22.80400000000 1.74830773908 5088.62883976680; 22.32400000000 0.93932040730 951.71840625060; 21.70800000000 3.83571581352 6283.07584999140; 21.63100000000 4.56895741061 3532.06069281140; 21.30400000000 0.78049229782 1589.07289528380; 20.42800000000 3.13540712557 4690.47983635860; 18.23700000000 0.41328624131 5486.77784317500; 17.95600000000 4.21930481803 3870.30339179440; 16.85000000000 4.53690440252 4292.33083295040; 16.80300000000 5.54857987615 3097.88382272579; 16.52800000000 0.96752074938 4399.99435688900; 16.46300000000 3.53882915214 2700.71514038580; 16.25100000000 3.80760134974 3340.54511639700; 16.25100000000 3.39910907599 3340.67973700260; 16.17100000000 2.34870953554 553.56940284240; 15.75500000000 4.75736730681 9492.14631500480; 15.74600000000 3.72356090283 20043.67456019880; 14.69900000000 5.95325006816 3894.18182954220; 14.25900000000 3.99897353022 1990.74501704100; 13.16900000000 0.41461716663 5614.72937620960; 13.01000000000 5.14230107067 6677.70173505060; 12.49200000000 1.03211063742 3341.59274776800; 12.40800000000 6.23142869816 5628.95647021120; 11.27200000000 1.02375627844 12303.06777661000 ]; r2 = [ 44242.24700000000 0.47930603943 3340.61242669980; 8138.04200000000 0.86998398093 6681.22485339960; 1274.91500000000 1.22594050809 10021.83728009940; 187.38700000000 1.57298991982 13362.44970679920; 52.39600000000 3.14159265359 0.00000000000; 40.74400000000 1.97080175060 3344.13554504880; 26.61600000000 1.91665615762 16703.06213349900; 17.82500000000 4.43499505333 2281.23049651060; 11.71300000000 4.52510453730 3185.19202726560; 10.20900000000 5.39143469548 1059.38193018920; 9.95000000000 0.41870577185 796.29800681640; 9.23700000000 4.53579272961 2146.16541647520; 7.78500000000 5.93369079547 1748.01641306700; 7.29900000000 3.14218509183 2544.31441988340; 7.21700000000 2.29300859074 6684.74797174860; 6.80800000000 5.26702580055 155.42039943420; 6.74900000000 5.30194395749 1194.44701022460; 6.52800000000 2.30781369329 3738.76143010800; 5.84000000000 1.05191350362 1349.86740965880; 5.39100000000 1.00223256041 3149.16416058820; 4.69500000000 0.76880586144 3097.88382272579; 4.40600000000 2.45556303355 951.71840625060; 4.28600000000 3.89643520638 1592.59601363280 ]; r3 = [ 1113.10700000000 5.14987350142 3340.61242669980; 424.44600000000 5.61343766478 6681.22485339960; 100.04400000000 5.99726827028 10021.83728009940; 19.60600000000 0.07633062094 13362.44970679920; 4.69300000000 3.14159265359 0.00000000000; 3.47700000000 0.42951907576 16703.06213349900; 2.86900000000 0.44711842697 3344.13554504880; 2.42800000000 3.02115527957 3185.19202726560; 0.68800000000 0.80693359444 6684.74797174860; 0.57700000000 0.77853275120 20043.67456019880; 0.54000000000 3.86836515672 1059.38193018920; 0.48700000000 1.60862391345 3583.34103067380; 0.46800000000 4.52450786544 3496.03282613400; 0.39700000000 5.71967986581 3149.16416058820; 0.36200000000 4.42397903418 2787.04302385740 ]; r4 = [ 19.55200000000 3.58211650473 3340.61242669980; 16.32300000000 4.05116076923 6681.22485339960; 5.84800000000 4.46383962094 10021.83728009940; 1.53200000000 4.84374321619 13362.44970679920; 0.37500000000 1.50962233608 3185.19202726560; 0.33900000000 5.20684967613 16703.06213349900; 0.15100000000 5.16472931648 3344.13554504880; 0.14800000000 0.00000000000 0.00000000000; 0.12500000000 2.19233532803 3496.03282613400; 0.08700000000 0.10275067375 3583.34103067380 ]; r5 = [ 19.55200000000 2.47617204701 6681.22485339960; 16.32300000000 2.91510547706 10021.83728009940; 5.84800000000 1.76888962311 3340.61242669980; 1.53200000000 3.31378377179 13362.44970679920 ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); sr5 = zzterm(r5, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4 + sr5*t5)/100000000; return [lambda, beta, delta]; endfunction function hmercury(day) ## returns the heliocentric ecliptic latitude of Mercury given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## More terms than Meeus included as an experiment .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC l0 = [ 440250710.144 0.00000000000 0.0000000000; 40989414.976 1.48302034194 26087.9031415742; 5046294.199 4.47785489540 52175.8062831484; 855346.843 1.16520322351 78263.7094247225; 165590.362 4.11969163181 104351.6125662960; 34561.897 0.77930765817 130439.5157078700; 7583.476 3.71348400510 156527.4188494450; 3559.740 1.51202669419 1109.3785520934; 1726.012 0.35832239908 182615.3219910190; 1803.463 4.10333178410 5661.3320491522; 1364.682 4.59918318745 27197.2816936676; 1589.923 2.99510417815 25028.5212113850; 1017.332 0.88031439040 31749.2351907264; 714.182 1.54144865265 24978.5245894808; 643.759 5.30266110787 21535.9496445154; 404.200 3.28228847025 208703.2251325930; 352.441 5.24156297101 20426.5710924220; 343.313 5.76531885335 955.5997416086; 339.214 5.86327765000 25558.2121764796; 451.137 6.04989275289 51116.4243529592; 325.335 1.33674334780 53285.1848352418; 259.587 0.98732428184 4551.9534970588; 345.212 2.79211901539 15874.6175953632; 272.947 2.49451163975 529.6909650946; 234.830 0.26672118900 11322.6640983044; 238.793 0.11343953378 1059.3819301892; 264.336 3.91705094013 57837.1383323006; 216.645 0.65987207348 13521.7514415914; 183.359 2.62878670784 27043.5028831828; 175.965 4.53636829858 51066.4277310550; 181.629 2.43413502466 25661.3049506982; 208.995 2.09178234008 47623.8527860896; 172.643 2.45200164173 24498.8302462904; 142.316 3.36003948842 37410.5672398786; 137.942 0.29098447849 10213.2855462110; 118.233 2.78149786369 77204.3274945333; 96.860 6.20398202740 234791.1282741670; 125.219 3.72079804425 39609.6545831656; 86.819 2.64219349385 51646.1153180537 ]; l1 = [ 2608814706222.740 0.00000000000 0.0000000000; 1126007.832 6.21703970996 26087.9031415742; 303471.395 3.05565472363 52175.8062831484; 80538.452 6.10454743366 78263.7094247225; 21245.035 2.83531934452 104351.6125662960; 5592.094 5.82675673328 130439.5157078700; 1472.233 2.51845458395 156527.4188494450; 352.244 3.05238094403 1109.3785520934; 388.318 5.48039225891 182615.3219910190; 93.540 6.11791163931 27197.2816936676; 90.579 0.00045481669 24978.5245894808; 102.743 2.14879173777 208703.2251325930; 51.941 5.62107554052 5661.3320491522; 44.370 4.57348500464 25028.5212113850; 28.070 3.04195430989 51066.4277310550; 22.003 0.86475371243 955.5997416086; 27.295 5.09210138837 234791.1282741670; 20.425 3.71509622702 20426.5710924220 ]; l2 = [ 53049.845 0.00000000000 0.0000000000; 16903.658 4.69072300649 26087.9031415742; 7396.711 1.34735624669 52175.8062831484; 3018.297 4.45643539705 78263.7094247225; 1107.419 1.26226537554 104351.6125662960; 378.173 4.31998055900 130439.5157078700; 122.998 1.06868541052 156527.4188494450; 38.663 4.08011610182 182615.3219910190; 14.898 4.63343085810 1109.3785520934; 11.861 0.79187646439 208703.2251325930; 5.243 4.71799772791 24978.5245894808; 3.575 3.77317513032 234791.1282741670 ]; l3 = [ 188.077 0.03466830117 52175.8062831484; 142.152 3.12505452600 26087.9031415742; 96.877 3.00378171915 78263.7094247225; 43.669 6.01867965826 104351.6125662960; 35.395 0.00000000000 0.0000000000; 18.045 2.77538373991 130439.5157078700; 6.971 5.81808665742 156527.4188494450; 2.556 2.57014364454 182615.3219910190; 0.900 5.59308888939 208703.2251325930 ]; l4 = [ 114.078 3.14159265359 0.0000000000; 3.247 2.02848007619 26087.9031415742; 1.914 1.41731803758 78263.7094247225; 1.727 4.50137643801 52175.8062831484; 1.237 4.49970181057 104351.6125662960; 0.645 1.26591776986 130439.5157078700; 0.298 4.30600984981 156527.4188494450 ]; l5 = [ 0.877 3.14159265359 0.0000000000; 0.059 3.37513289692 52175.8062831484 ]; .. latitude series b0 = [ 11737528.962 1.98357498767 26087.9031415742; 2388076.996 5.03738959685 52175.8062831484; 1222839.532 3.14159265359 0.0000000000; 543251.810 1.79644363963 78263.7094247225; 129778.770 4.83232503961 104351.6125662960; 31866.927 1.58088495667 130439.5157078700; 7963.301 4.60972126348 156527.4188494450; 2014.189 1.35324164694 182615.3219910190; 513.953 4.37835409309 208703.2251325930; 207.674 4.91772564073 27197.2816936676; 208.584 2.02020294153 24978.5245894808; 132.013 1.11908492283 234791.1282741670; 100.454 5.65684734206 20426.5710924220; 121.395 1.81271752059 53285.1848352418; 91.566 2.28163128692 25028.5212113850; 99.214 0.09391887097 51116.4243529592 ]; b1 = [ 429151.362 3.50169780393 26087.9031415742; 146233.668 3.14159265359 0.0000000000; 22675.295 0.01515366880 52175.8062831484; 10894.981 0.48540174006 78263.7094247225; 6353.462 3.42943919982 104351.6125662960; 2495.743 0.16051210665 130439.5157078700; 859.585 3.18452433647 156527.4188494450; 277.503 6.21020774184 182615.3219910190; 86.233 2.95244391822 208703.2251325930; 26.133 5.97708962692 234791.1282741670; 27.696 0.29068938889 27197.2816936676; 12.831 3.37744320558 53285.1848352418; 12.720 0.53792661684 24978.5245894808; 7.781 2.71768609268 260879.0314157410 ]; b2 = [ 11830.934 4.79065585784 26087.9031415742; 1913.516 0.00000000000 0.0000000000; 1044.801 1.21216540536 52175.8062831484; 266.213 4.43418336532 78263.7094247225; 170.280 1.62255638714 104351.6125662960; 96.300 4.80023692017 130439.5157078700; 44.692 1.60758267772 156527.4188494450; 18.316 4.66904655377 182615.3219910190; 6.927 1.43404888930 208703.2251325930; 2.479 4.47495202955 234791.1282741670; 1.739 1.83080039600 27197.2816936676; 0.852 1.22749255198 260879.0314157410; 0.641 4.87358642253 53285.1848352418 ]; b3 = [ 235.423 0.35387524604 26087.9031415742; 160.537 0.00000000000 0.0000000000; 18.904 4.36275460261 52175.8062831484; 6.376 2.50715381439 78263.7094247225; 4.580 6.14257817571 104351.6125662960; 3.061 3.12497552681 130439.5157078700; 1.732 6.26642412058 156527.4188494450; 0.857 3.07673166705 182615.3219910190; 0.384 6.14815319932 208703.2251325930; 0.159 2.92437378320 234791.1282741670 ]; b4 = [ 4.276 1.74579932115 26087.9031415742; 1.023 3.14159265359 0.0000000000; 0.425 4.03419509143 52175.8062831484 ]; b5 = [ 0.106 3.94555784256 26087.90314157420 ]; .. radius vector r0 = [ 39528271.652 0.00000000000 0.0000000000; 7834131.817 6.19233722599 26087.9031415742; 795525.557 2.95989690096 52175.8062831484; 121281.763 6.01064153805 78263.7094247225; 21921.969 2.77820093975 104351.6125662960; 4354.065 5.82894543257 130439.5157078700; 918.228 2.59650562598 156527.4188494450; 260.033 3.02817753482 27197.2816936676; 289.955 1.42441936951 25028.5212113850; 201.855 5.64725040350 182615.3219910190; 201.499 5.59227724202 31749.2351907264; 141.980 6.25264202645 24978.5245894808; 100.144 3.73435608689 21535.9496445154; 77.561 3.66972526976 20426.5710924220; 63.277 4.29905918105 25558.2121764796; 62.951 4.76588899933 1059.3819301892; 66.754 2.52520309182 5661.3320491522 ]; r1 = [ 217347.739 4.65617158663 26087.9031415742; 44141.826 1.42385543975 52175.8062831484; 10094.479 4.47466326316 78263.7094247225; 2432.804 1.24226083435 104351.6125662960; 1624.367 0.00000000000 0.0000000000; 603.996 4.29303116561 130439.5157078700; 152.851 1.06060779810 156527.4188494450; 39.202 4.11136751416 182615.3219910190; 17.760 4.54424653085 27197.2816936676; 17.999 4.71193725810 24978.5245894808; 10.154 0.87893548494 208703.2251325930 ]; r2 = [ 3117.867 3.08231840296 26087.9031415742; 1245.396 6.15183317423 52175.8062831484; 424.822 2.92583352960 78263.7094247225; 136.130 5.97983925842 104351.6125662960; 42.175 2.74936980629 130439.5157078700; 21.759 3.14159265359 0.0000000000; 12.793 5.80143162209 156527.4188494450; 3.825 2.56993599584 182615.3219910190; 1.042 3.14648120079 24978.5245894808; 1.131 5.62142196970 208703.2251325930; 0.483 6.14307654520 27197.2816936676 ]; r3 = [ 32.676 1.67971635359 26087.9031415742; 24.166 4.63403168997 52175.8062831484; 12.133 1.38983781545 78263.7094247225; 5.140 4.43915386930 104351.6125662960 ]; r4 = [ 0.394 0.36735403840 26087.9031415742 ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4)/100000000; return [lambda, beta, delta]; endfunction function hsaturn(day) ## returns the heliocentric ecliptic latitude of Saturn given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## More terms than Meeus included as an experiment .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC to a resolution .. of 1 arcsec as predicted using Meeus truncation formula l0 = [ 87401354.029 0 0.0000000000; 11107659.780 3.96205090194 213.2990954380; 1414150.958 4.58581515873 7.1135470008; 398379.386 0.52112025957 206.1855484372; 350769.223 3.30329903015 426.5981908760; 206816.296 0.24658366938 103.0927742186; 79271.288 3.84007078530 220.4126424388; 23990.338 4.66976934860 110.2063212194; 16573.583 0.43719123541 419.4846438752; 14906.995 5.76903283845 316.3918696566; 15820.300 0.93808953760 632.7837393132; 14609.562 1.56518573691 3.9321532631; 13160.308 4.44891180176 14.2270940016; 15053.509 2.71670027883 639.8972863140; 13005.305 5.98119067061 11.0457002639; 10725.066 3.12939596466 202.2533951741; 5863.207 0.23657028777 529.6909650946; 5227.771 4.20783162380 3.1813937377; 6126.308 1.76328499656 277.0349937414; 5019.658 3.17787919533 433.7117378768; 4592.541 0.61976424374 199.0720014364; 4005.862 2.24479893937 63.7358983034; 2953.815 0.98280385206 95.9792272178; 3873.696 3.22282692566 138.5174968707; 2461.172 2.03163631205 735.8765135318; 3269.490 0.77491895787 949.1756089698; 1758.143 3.26580514774 522.5774180938; 1640.183 5.50504966218 846.0828347512; 1391.336 4.02331978116 323.5054166574; 1580.641 4.37266314120 309.2783226558; 1123.515 2.83726793572 415.5524906121; 1017.258 3.71698151814 227.5261894396; 848.643 3.19149825839 209.3669421749; 1087.237 4.18343232481 2.4476805548; 956.752 0.50740889886 1265.5674786264; 789.205 5.00745123149 0.9632078465; 686.965 1.74714407827 1052.2683831884; 654.470 1.59889331515 0.0481841098; 748.811 2.14398149298 853.1963817520; 633.980 2.29889903023 412.3710968744; 743.584 5.25276954625 224.3447957019; 852.677 3.42141350697 175.1660598002; 579.857 3.09259007048 74.7815985673; 624.904 0.97046831256 210.1177017003; 529.861 4.44938897119 117.3198682202; 542.643 1.51824320514 9.5612275556; 474.279 5.47527185987 742.9900605326; 448.542 1.28990416161 127.4717966068; 546.358 2.12678554211 350.3321196004; 478.054 2.96488054338 137.0330241624; 354.944 3.01286483030 838.9692877504; 451.827 1.04436664241 490.3340891794; 347.413 1.53928227764 340.7708920448; 343.475 0.24604039134 0.5212648618; 309.001 3.49486734909 216.4804891757; 322.185 0.96137456104 203.7378678824; 372.308 2.27819108625 217.2312487011; 321.543 2.57182354537 647.0108333148; 330.196 0.24715617844 1581.9593482830; 249.116 1.47010534421 1368.6602528450; 286.688 2.37043745859 351.8165923087; 220.225 4.20422424873 200.7689224658; 277.775 0.40020408926 211.8146227297; 204.500 6.01082206600 265.9892934775; 207.663 0.48349820488 1162.4747044078; 208.655 1.34516255304 625.6701923124; 182.454 5.49122292426 2.9207613068; 226.609 4.91003163138 12.5301729722; 207.659 1.28302218900 39.3568759152; 173.914 1.86305806814 0.7507595254; 184.690 3.50344404958 149.5631971346; 183.511 0.97254952728 4.1927856940; 146.068 6.23102544071 195.1398481733; 164.541 0.44005517520 5.4166259714; 147.526 1.53529320509 5.6290742925; 139.666 4.29450260069 21.3406410024; 131.283 4.06828961903 10.2949407385; 117.283 2.67920400584 1155.3611574070; 149.299 5.73594349789 52.6901980395; 122.373 1.97588777199 4.6658664460; 113.747 5.59427544714 1059.3819301892; 102.702 1.19748124058 1685.0521225016; 118.156 5.34072933900 554.0699874828; 109.275 3.43812715686 536.8045120954; 110.399 0.16604024090 1.4844727083; 124.969 6.27737805832 1898.3512179396; 89.949 5.80392934702 114.1384744825; 103.956 2.19210363069 88.8656802170; 112.437 1.10502663534 191.2076949102; 106.570 4.01156608514 956.2891559706; 91.430 1.87521577510 38.1330356378; 83.791 5.48810655641 0.1118745846; 83.461 2.28972767279 628.8515860501; 96.987 4.53666595763 302.1647756550; 100.631 4.96513666539 269.9214467406; 75.491 2.18045274099 728.7629665310; 96.330 2.83319189210 275.5505210331; 82.363 3.05469876064 440.8252848776; 73.888 5.08914205084 1375.7737998458; 71.633 5.10940743430 65.2203710117; 70.409 4.86846451411 0.2124483211; 69.760 3.71029022489 14.9778535270; 88.772 3.86334563977 278.5194664497; 68.090 0.73415460990 1478.8665740644; 66.501 0.02677580336 70.8494453042; 65.682 2.02165559602 142.4496501338; 75.765 1.61410487792 284.1485407422; 63.153 3.49493353034 479.2883889155; 62.539 2.58713611532 422.6660376129; 69.313 3.43979731402 515.4638710930; 79.021 4.45154941586 35.4247226521; 63.664 3.31749528708 62.2514255951; 52.939 5.51392725227 0.2606324309; 53.011 3.18480701697 8.0767548473; 54.492 2.45674090515 22.0914005278; 50.514 4.26749346978 99.1606209555; 55.170 0.96797446150 942.0620619690; 49.288 2.38641424063 1471.7530270636; 47.199 2.02515248245 312.1990839626; 61.080 1.50295092063 210.8514148832; 45.126 0.93109376473 2001.4439921582; 60.556 2.68715551585 388.4651552382; 43.452 2.52602011714 288.0806940053; 42.544 3.81793980322 330.6189636582; 39.915 5.71378652900 408.4389436113; 50.145 6.03164759907 2214.7430875962; 45.860 0.54229721801 212.3358875915; 54.165 0.78154835399 191.9584544356; 47.016 4.59934671151 437.6438911399; 42.362 1.90070070955 430.5303441391; 39.722 1.63259419913 1066.4954771900; 36.345 0.84756992711 213.3472795478; 35.468 4.18603772925 215.7467759928; 36.344 3.93295730315 213.2509113282; 38.005 0.31313803095 423.4167971383; 44.746 1.12488341174 6.1503391543; 37.902 1.19795851115 2.7083129857; 43.402 1.37363944007 563.6312150384; 43.764 3.93043802956 525.4981794006; 34.825 1.01566605408 203.0041546995; 31.755 1.69273634405 0.1600586944; 30.880 6.13525703832 417.0369633204; 36.388 6.00586032647 18.1592472647; 29.032 1.19660544505 404.5067903482; 32.812 0.53649479713 107.0249274817; 30.433 0.72335287989 222.8603229936; 32.644 0.81204701486 1795.2584437210; 37.769 3.69666903716 1272.6810256272; 27.679 1.45663979401 7.1617311106; 27.187 1.89731951902 1045.1548361876; 37.699 4.51997049537 24.3790223882; 34.885 4.46095761791 214.2623032845; 32.650 0.66372395761 692.5874843535; 30.324 5.30369950147 33.9402499438; 27.480 6.22702216249 1.2720243872; 26.657 4.56713198392 7.0653628910; 31.745 5.49798599565 56.6223513026; 28.050 5.64447420566 128.9562693151; 24.277 3.93966553574 414.0680179038; 32.017 5.22260660455 92.0470739547; 26.976 0.06705123981 205.2223405907; 22.974 3.65817751770 207.6700211455; 31.775 5.59198119173 6069.7767545534; 23.153 2.10054506119 1788.1448967202; 31.025 0.37190053329 703.6331846174; 29.376 0.14742155778 131.4039498699; 22.562 5.24009182383 212.7778305762; 26.185 5.41311252822 140.0019695790; 25.673 4.36038885283 32.2433289144; 20.392 2.82413909260 429.7795846137; 20.659 0.67091805084 2317.8358618148; 24.397 3.08740396398 145.6310438715; 23.735 2.54365387567 76.2660712756; 20.157 5.06708675157 617.8058857862; 23.307 3.97357729211 483.2205421786; 22.878 6.10452832642 177.8743727859; 22.978 3.20140795404 208.6332289920; 20.638 5.22128727027 6.5922821390; 21.446 0.72034565528 1258.4539316256; 18.034 6.11382719947 210.3783341312 ]; l1 = [ 21354295595.986 0.00000000000 0.0000000000; 1296855.005 1.82820544701 213.2990954380; 564347.566 2.88500136429 7.1135470008; 98323.030 1.08070061328 426.5981908760; 107678.770 2.27769911872 206.1855484372; 40254.586 2.04128257090 220.4126424388; 19941.734 1.27954662736 103.0927742186; 10511.706 2.74880392800 14.2270940016; 6939.233 0.40493079985 639.8972863140; 4803.325 2.44194097666 419.4846438752; 4056.325 2.92166618776 110.2063212194; 3768.630 3.64965631460 3.9321532631; 3384.684 2.41694251653 3.1813937377; 3302.200 1.26256486715 433.7117378768; 3071.382 2.32739317750 199.0720014364; 1953.036 3.56394683300 11.0457002639; 1249.348 2.62803737519 95.9792272178; 921.683 1.96089834250 227.5261894396; 705.587 4.41689249330 529.6909650946; 649.654 6.17418093659 202.2533951741; 627.603 6.11088227167 309.2783226558; 486.843 6.03998200305 853.1963817520; 468.377 4.61707843907 63.7358983034; 478.501 4.98776987984 522.5774180938; 417.010 2.11708169277 323.5054166574; 407.630 1.29949556676 209.3669421749; 343.826 3.95854178574 412.3710968744; 339.724 3.63396398752 316.3918696566; 335.936 3.77173072712 735.8765135318; 331.933 2.86077699882 210.1177017003; 352.489 2.31707079463 632.7837393132; 289.429 2.73263080235 117.3198682202; 265.801 0.54344631312 647.0108333148; 230.493 1.64428879621 216.4804891757; 280.911 5.74398845416 2.4476805548; 191.667 2.96512946582 224.3447957019; 172.891 4.07695221044 846.0828347512; 167.131 2.59745202658 21.3406410024; 136.328 2.28580246629 10.2949407385; 131.364 3.44108355646 742.9900605326; 127.838 4.09533471247 217.2312487011; 108.862 6.16141072262 415.5524906121; 93.909 3.48397279899 1052.2683831884; 92.482 3.94755499926 88.8656802170; 97.584 4.72845436677 838.9692877504; 86.600 1.21951325061 440.8252848776; 83.463 3.11269504725 625.6701923124; 77.588 6.24408938835 302.1647756550; 61.557 1.82789612597 195.1398481733; 61.900 4.29344363385 127.4717966068; 67.106 0.28961738595 4.6658664460; 56.919 5.01889578112 137.0330241624; 54.160 5.12628572382 490.3340891794; 54.585 0.28356341456 74.7815985673; 51.425 1.45766406064 536.8045120954; 65.843 5.64757042732 9.5612275556; 57.780 2.47630552035 191.9584544356; 44.444 2.70873627665 5.4166259714; 46.799 1.17721211050 149.5631971346; 40.380 3.88870105683 728.7629665310; 37.768 2.53379013859 12.5301729722; 46.649 5.14818326902 515.4638710930; 45.891 2.23198878761 956.2891559706; 40.400 0.41281520440 269.9214467406; 37.191 3.78239026411 2.9207613068; 33.778 3.21070688046 1368.6602528450; 37.969 0.64665967180 422.6660376129; 32.857 0.30063884563 351.8165923087; 33.050 5.43038091186 1066.4954771900; 30.276 2.84067004928 203.0041546995; 35.116 6.08421794089 5.6290742925; 29.667 3.39052569135 1059.3819301892; 33.217 4.64063092111 277.0349937414; 31.876 4.38622923770 1155.3611574070; 28.913 2.02614760507 330.6189636582; 28.264 2.74178953996 265.9892934775; 30.089 6.18684614308 284.1485407422; 31.329 2.43455855525 52.6901980395; 26.493 4.51214170121 340.7708920448; 21.983 5.14437352579 4.1927856940; 22.230 1.96481952451 203.7378678824; 20.824 6.16048095923 860.3099287528; 21.690 2.67578768862 942.0620619690; 22.552 5.88579123000 210.8514148832; 19.807 2.31345263487 437.6438911399; 19.447 4.76573277668 70.8494453042; 19.310 4.10209060369 18.1592472647; 22.662 4.13732273379 191.2076949102; 18.209 0.90310796389 429.7795846137; 17.667 1.84954766042 234.6397364404; 17.547 2.44735118493 423.4167971383; 15.428 4.23790088205 1162.4747044078; 14.608 3.59713247857 1045.1548361876; 14.111 2.94262468353 1685.0521225016; 16.328 4.05665272725 949.1756089698; 13.348 6.24509592240 38.1330356378; 15.918 1.06434204938 56.6223513026; 14.059 1.43503954068 408.4389436113; 13.093 5.75815864257 138.5174968707; 15.772 5.59350835225 6.1503391543; 14.962 5.77192239389 22.0914005278; 16.024 1.93900586533 1272.6810256272; 16.751 5.96673627422 628.8515860501; 12.843 4.24658666814 405.2575498736; 13.628 4.09892958087 1471.7530270636; 15.067 0.74142807591 200.7689224658; 10.961 1.55022573283 223.5940361765; 11.695 1.81237511034 124.4334152210; 10.346 3.46814088412 1375.7737998458; 12.056 1.85655834555 131.4039498699; 10.123 2.38221133049 107.0249274817; 9.855 3.95166998848 430.5303441391; 9.803 2.55389483994 99.9113804809; 10.614 5.36692189034 215.7467759928; 12.080 4.84549317054 831.8557407496; 10.210 6.07692961370 32.2433289144; 9.245 3.65417467270 142.4496501338; 8.984 1.23808405498 106.2741679563; 9.336 5.81062768434 7.1617311106; 9.717 1.38703872827 145.6310438715; 8.394 4.42341211111 703.6331846174; 8.370 5.64015188458 62.2514255951; 8.244 2.42225929772 1258.4539316256; 7.784 0.52562994711 654.1243803156; 7.626 3.75258725596 312.1990839626; 7.222 0.28429555677 0.7507595254; 8.236 6.22250515902 14.9778535270; 7.054 0.53177810740 388.4651552382; 6.567 3.48657341701 35.4247226521; 9.011 4.94919626910 208.6332289920; 8.980 0.08138173719 288.0806940053; 6.421 3.32905264657 1361.5467058442; 6.489 2.89389587598 114.1384744825; 6.244 0.54973852782 65.2203710117; 6.154 2.67885860584 2001.4439921582; 6.742 0.23586769279 8.0767548473; 7.297 4.85321224483 222.8603229936; 6.302 3.80651124694 1788.1448967202; 5.824 4.39327457448 81.7521332162; 6.102 0.88585782895 92.0470739547; 6.914 2.04631426723 99.1606209555; 5.363 5.47995103139 563.6312150384; 5.172 2.11968421583 214.2623032845; 5.117 5.76987684107 565.1156877467; 6.197 1.62553688800 1589.0728952838; 4.970 0.41949366126 76.2660712756; 6.640 5.82582210639 483.2205421786; 5.277 4.57975789757 134.5853436076; 4.974 4.20243895902 404.5067903482; 5.150 4.67582673243 212.3358875915; 4.764 4.59303997414 554.0699874828; 4.573 3.24875415786 231.4583427027; 4.811 0.46206327592 362.8622925726; 5.148 3.33570646174 1.4844727083; 4.654 5.80233066659 217.9649618840; 4.509 5.37581684215 497.4476361802; 4.443 0.11349392292 295.0512286542; 4.943 3.78020789259 1265.5674786264; 4.211 4.88306021960 98.8999885246; 4.252 5.00120115113 213.3472795478; 4.774 4.53259894142 1148.2476104062; 3.911 0.58582192963 750.1036075334; 5.069 2.20305668335 207.8824694666; 3.553 0.35374030841 333.6573450440; 3.771 0.98542435766 24.3790223882; 3.458 1.84990273999 225.8292684102; 3.401 5.31342401626 347.8844390456; 3.347 0.21414641376 635.9651330509; 3.637 1.61315058382 245.5424243524; 3.416 2.19551489078 1574.8458012822; 3.655 0.80544245690 343.2185725996; 4.260 1.80258750109 213.2509113282; 3.110 3.03815175282 1677.9385755008; 3.052 1.33858964447 543.9180590962; 3.694 0.81606028298 344.7030453079; 3.016 3.36219319026 7.8643065262; 2.937 4.86927342776 144.1465711632; 2.768 2.42707131609 2317.8358618148; 3.059 4.30820099442 6062.6632075526; 3.650 5.12802531219 218.9281697305; 2.963 3.53480751374 2104.5367663768; 3.230 2.88057019783 216.2198567448; 2.984 2.52971310583 1692.1656695024; 2.897 5.73256482240 9992.8729037722; 2.591 3.79880285744 17.2654753874; 3.495 5.29902525443 350.3321196004; 2.859 3.72804950659 6076.8903015542; 2.775 0.23549396237 357.4456666012; 2.976 2.48769315964 46.4704229160; 2.487 4.37868078530 217.4918811320; 2.711 5.15376840150 10007.0999977738; 3.127 1.92343235583 17.4084877393; 3.181 1.72419900322 1169.5882514086; 2.348 0.77373103004 414.0680179038; 2.606 3.42836913440 31.0194886370; 2.556 0.91735028377 479.2883889155; 2.399 4.82440545738 1279.7945726280; 2.245 3.76323995584 425.1137181677; 3.020 0.25310250109 120.3582496060; 2.503 2.10679832121 168.0525127994; 2.564 1.63158205055 182.2796068010 ]; l2 = [ 116441.181 1.17987850633 7.1135470008; 91920.844 0.07425261094 213.2990954380; 90592.251 0.00000000000 0.0000000000; 15276.909 4.06492007503 206.1855484372; 10631.396 0.25778277414 220.4126424388; 10604.979 5.40963595885 426.5981908760; 4265.368 1.04595556630 14.2270940016; 1215.527 2.91860042123 103.0927742186; 1164.684 4.60942128971 639.8972863140; 1081.967 5.69130351670 433.7117378768; 1020.079 0.63369182642 3.1813937377; 1044.754 4.04206453611 199.0720014364; 633.582 4.38825410036 419.4846438752; 549.329 5.57303134242 3.9321532631; 456.914 1.26840971349 110.2063212194; 425.100 0.20935499279 227.5261894396; 273.739 4.28841011784 95.9792272178; 161.571 1.38139149420 11.0457002639; 129.494 1.56586884170 309.2783226558; 117.008 3.88120915956 853.1963817520; 105.415 4.90003203599 647.0108333148; 100.967 0.89270493100 21.3406410024; 95.227 5.62561150598 412.3710968744; 81.948 1.02477558315 117.3198682202; 74.857 4.76178468163 210.1177017003; 82.727 6.05030934786 216.4804891757; 95.659 2.91093561539 316.3918696566; 63.696 0.35179804917 323.5054166574; 84.860 5.73472777961 209.3669421749; 60.647 4.87517850190 632.7837393132; 66.459 0.48297940601 10.2949407385; 67.184 0.45648612616 522.5774180938; 53.281 2.74730541387 529.6909650946; 45.827 5.69296621745 440.8252848776; 45.293 1.66856699796 202.2533951741; 42.330 5.70768187703 88.8656802170; 32.140 0.07050050346 63.7358983034; 31.573 1.67190022213 302.1647756550; 31.150 4.16379537691 191.9584544356; 24.631 5.65564728570 735.8765135318; 26.558 0.83256214407 224.3447957019; 20.108 5.94364609981 217.2312487011; 17.511 4.90014736798 625.6701923124; 17.130 1.62593421274 742.9900605326; 13.744 3.76497167300 195.1398481733; 12.236 4.71789723976 203.0041546995; 11.940 0.12620714199 234.6397364404; 16.040 0.57886320845 515.4638710930; 11.154 5.92216844780 536.8045120954; 14.068 0.20675293700 838.9692877504; 11.013 5.60207982774 728.7629665310; 11.718 3.12098483554 846.0828347512; 9.962 4.15472049127 860.3099287528; 10.601 3.20327613035 1066.4954771900; 10.072 0.25709351996 330.6189636582; 9.490 0.46379969328 956.2891559706; 10.240 4.98736656070 422.6660376129; 8.287 2.13990364272 269.9214467406; 7.238 5.39724715258 1052.2683831884; 7.730 5.24602742309 429.7795846137; 6.353 4.46211130731 284.1485407422; 5.935 5.40967847103 149.5631971346; 7.550 4.03401153929 9.5612275556; 5.779 4.29380891110 415.5524906121; 6.082 5.93416924841 405.2575498736; 5.711 0.01824076994 124.4334152210; 5.676 6.02235682150 223.5940361765; 4.757 4.92804854717 654.1243803156; 4.727 2.27461984667 18.1592472647; 4.509 4.40688707557 942.0620619690; 5.621 0.29694719379 127.4717966068; 5.453 5.53868222772 949.1756089698; 4.130 4.68673560379 74.7815985673; 4.098 5.30851262200 1045.1548361876; 4.223 2.89014939299 56.6223513026; 4.887 3.20022991216 277.0349937414; 3.905 3.30270187305 490.3340891794; 3.923 6.09732996823 81.7521332162; 3.755 4.93065184796 52.6901980395; 4.602 6.13908576681 1155.3611574070; 3.714 0.40648076787 137.0330241624; 3.407 4.28514461015 99.9113804809; 3.579 0.20402442077 1272.6810256272; 3.946 0.36500928968 12.5301729722; 3.246 1.56761884227 1059.3819301892; 4.063 0.29084229143 831.8557407496; 3.688 0.15467406177 437.6438911399; 2.895 3.13473183482 70.8494453042; 2.800 0.32727938074 191.2076949102; 2.672 1.87612402267 295.0512286542; 3.454 4.77197610696 423.4167971383; 2.623 5.15237415384 1368.6602528450; 2.457 3.89612890177 210.8514148832; 2.461 1.58522876760 32.2433289144; 2.595 3.59007068361 131.4039498699; 2.289 4.76825865118 351.8165923087; 2.357 5.83099000562 106.2741679563; 2.221 5.98277491515 6062.6632075526; 2.221 2.05930402282 6076.8903015542; 2.183 5.94985336393 145.6310438715; 2.718 3.37801252354 408.4389436113; 2.288 3.14000619320 22.0914005278; 2.090 1.12304173562 9992.8729037722; 2.089 3.48276230686 10007.0999977738; 2.570 5.12167203704 265.9892934775; 1.835 4.15379879659 1258.4539316256; 1.820 5.05340615445 1361.5467058442; 1.760 4.13532689228 107.0249274817; 1.921 4.51790997496 138.5174968707; 1.707 1.35864593280 231.4583427027; 1.956 5.87006093798 1471.7530270636; 2.133 5.23409848720 1265.5674786264; 1.595 5.61962698786 447.9388318784; 1.609 3.74893709671 628.8515860501; 1.490 0.48352404940 340.7708920448; 1.560 5.97095003614 430.5303441391; 1.352 0.71405348653 28.4541880032; 1.355 2.91219493604 215.7467759928; 1.298 5.84254169775 543.9180590962; 1.664 6.23834873469 1148.2476104062; 1.205 2.83373725021 200.7689224658; 1.192 3.52219428945 497.4476361802; 1.122 2.60571030270 1279.7945726280; 1.217 6.23528359211 1589.0728952838; 1.420 0.85079202155 6069.7767545534; 1.120 4.95656566453 1685.0521225016; 1.010 3.39689646619 1073.6090241908; 1.352 2.27575429523 9999.9864507730; 0.979 1.58571463442 1375.7737998458; 1.159 0.71823181781 508.3503240922; 1.014 2.40759054741 703.6331846174; 0.956 2.66256831556 134.5853436076; 1.110 1.19713920197 618.5566453116; 0.945 4.68155456977 362.8622925726; 0.953 4.20749172571 288.0806940053; 1.033 1.08781255146 184.8449074348; 0.942 2.43465223460 222.8603229936; 0.909 4.51769385360 38.1330356378; 1.002 1.38543153271 483.2205421786 ]; l3 = [ 16038.734 5.73945377424 7.1135470008; 4249.793 4.58539675603 213.2990954380; 1906.524 4.76082050205 220.4126424388; 1465.687 5.91326678323 206.1855484372; 1162.041 5.61973132428 14.2270940016; 1066.581 3.60816533142 426.5981908760; 239.377 3.86088273439 433.7117378768; 236.975 5.76826451465 199.0720014364; 165.641 5.11641150216 3.1813937377; 131.409 4.74327544615 227.5261894396; 151.352 2.73594641861 639.8972863140; 61.630 4.74287052463 103.0927742186; 63.365 0.22850089497 419.4846438752; 40.437 5.47298059144 21.3406410024; 40.205 5.96420266720 95.9792272178; 38.746 5.83386199529 110.2063212194; 28.025 3.01235311514 647.0108333148; 25.029 0.98808170740 3.9321532631; 18.101 1.02506397063 412.3710968744; 17.879 3.31913418974 309.2783226558; 16.208 3.89825272754 440.8252848776; 15.763 5.61667809625 117.3198682202; 19.014 1.91614237463 853.1963817520; 18.262 4.96738415934 10.2949407385; 12.947 1.18068953942 88.8656802170; 17.919 4.20376505349 216.4804891757; 11.453 5.57520615096 11.0457002639; 10.548 5.92906266269 191.9584544356; 10.389 3.94838736947 209.3669421749; 8.650 3.39335369698 302.1647756550; 7.580 4.87736913157 323.5054166574; 6.697 0.38198725552 632.7837393132; 5.864 1.05621157685 210.1177017003; 5.449 4.64268475485 234.6397364404; 6.327 2.25492722762 522.5774180938; 3.602 2.30677010956 515.4638710930; 3.229 2.20309400066 860.3099287528; 3.701 3.14159265359 0.0000000000; 2.583 4.93447677059 224.3447957019; 2.543 0.42393884183 625.6701923124; 2.213 3.19814958289 202.2533951741; 2.421 4.76621391814 330.6189636582; 2.850 0.58604395010 529.6909650946; 1.965 4.39525359412 124.4334152210; 2.154 1.35488209144 405.2575498736; 2.296 3.34809165905 429.7795846137; 2.018 3.06693569701 654.1243803156; 1.979 1.02981005658 728.7629665310; 1.868 3.09383546177 422.6660376129; 1.846 4.15225985450 536.8045120954; 2.194 1.18918501013 1066.4954771900; 2.090 4.15631351317 223.5940361765; 1.481 0.37916705169 316.3918696566; 1.720 5.82865773356 195.1398481733; 1.460 1.57663426355 81.7521332162; 1.623 6.03706764648 742.9900605326; 1.286 1.66154726117 63.7358983034; 1.304 5.02409881054 956.2891559706; 1.446 2.10575519127 838.9692877504; 1.245 3.88109752770 269.9214467406; 1.018 3.72599601656 295.0512286542; 1.323 1.38492882986 735.8765135318; 1.318 2.33460998999 217.2312487011; 0.943 2.75813531246 284.1485407422; 0.906 0.71155526266 846.0828347512; 0.886 3.83754799777 447.9388318784; 0.943 3.31480217015 18.1592472647; 0.800 4.71386673963 56.6223513026; 0.908 2.02119147951 831.8557407496; 0.787 0.80410269937 1045.1548361876; 0.709 4.27064410504 437.6438911399; 0.651 6.17565900032 942.0620619690; 0.785 2.40767785311 203.0041546995; 0.702 1.64585301418 423.4167971383; 0.543 2.86326941725 184.8449074348; 0.532 6.25762144463 1059.3819301892; 0.521 3.43013038466 149.5631971346; 0.484 4.88366060720 1272.6810256272; 0.437 5.40220619672 408.4389436113; 0.388 2.57589594168 508.3503240922; 0.421 4.05836524024 543.9180590962; 0.375 1.22747948298 2324.9494088156; 0.347 0.59237194522 22.0914005278; 0.433 1.69090148012 1155.3611574070; 0.389 1.46170367972 1073.6090241908; 0.307 1.82185086955 628.8515860501; 0.409 1.21858750514 1052.2683831884; 0.309 0.33610530663 6076.8903015542; 0.309 1.42279282226 6062.6632075526; 0.340 1.83325770310 1141.1340634054; 0.303 2.41584747330 127.4717966068; 0.305 5.34154702988 131.4039498699; 0.298 2.28594631393 635.9651330509; 0.372 1.03723911390 313.2104759189; 0.338 0.69100012338 1361.5467058442; 0.325 1.78816356937 1148.2476104062; 0.322 1.18628805010 721.6494195302 ]; l4 = [ 1661.894 3.99826248978 7.1135470008; 257.107 2.98436499013 220.4126424388; 236.344 3.90241428075 14.2270940016; 149.418 2.74110824208 213.2990954380; 109.598 1.51515739251 206.1855484372; 113.953 3.14159265359 0.0000000000; 68.390 1.72120953337 426.5981908760; 37.699 1.23795458356 199.0720014364; 40.060 2.04644897412 433.7117378768; 31.219 3.01094184090 227.5261894396; 15.111 0.82897064529 639.8972863140; 9.444 3.71485300868 21.3406410024; 5.690 2.41995290633 419.4846438752; 4.470 1.45120818748 95.9792272178; 5.608 1.15607095740 647.0108333148; 4.463 2.11783225176 440.8252848776; 3.229 4.09278077834 110.2063212194; 2.871 2.77203153866 412.3710968744; 2.796 3.00730249564 88.8656802170; 2.638 0.00255721254 853.1963817520; 2.574 0.39246854091 103.0927742186; 1.862 5.07955457727 309.2783226558; 2.225 3.77689198137 117.3198682202; 1.769 5.19176876406 302.1647756550; 1.921 2.82884328662 234.6397364404; 1.805 2.23816036743 216.4804891757; 1.211 1.54685246534 191.9584544356; 0.765 3.44501766503 323.5054166574; 0.763 4.83197222448 210.1177017003; 0.613 4.19052656353 515.4638710930; 0.648 2.28591710303 209.3669421749; 0.616 4.03194472161 522.5774180938; 0.630 2.37952532019 632.7837393132; 0.639 0.29772678242 860.3099287528; 0.559 2.17110060530 124.4334152210; 0.442 2.23500083592 447.9388318784; 0.407 5.44515970990 1066.4954771900; 0.469 1.26889429317 654.1243803156; 0.488 3.20329778617 405.2575498736; 0.415 3.12435410343 330.6189636582; 0.442 3.38933498625 81.7521332162; 0.332 4.12464206608 838.9692877504; 0.320 3.18332026736 529.6909650946; 0.312 1.40962796637 429.7795846137; 0.291 3.18885372262 1464.6394800628; 0.333 2.94355912397 728.7629665310; 0.235 3.67049647573 1148.2476104062; 0.286 2.57895004576 1045.1548361876; 0.223 3.57980034401 1155.3611574070; 0.261 2.04564143519 1677.9385755008; 0.218 2.61967125327 536.8045120954; 0.262 2.48322150677 625.6701923124; 0.191 4.39064160974 1574.8458012822; 0.176 1.26161895188 422.6660376129; 0.190 2.32693171200 223.5940361765; 0.185 1.08713469614 742.9900605326; 0.168 0.69946458053 824.7421937488 ]; l5 = [ 123.615 2.25923345732 7.1135470008; 34.190 2.16250652689 14.2270940016; 27.546 1.19868150215 220.4126424388; 5.818 1.21584270184 227.5261894396; 5.318 0.23550400093 433.7117378768; 3.677 6.22669694355 426.5981908760; 3.057 2.97372046322 199.0720014364; 2.861 4.28710932685 206.1855484372; 1.617 6.25265362286 213.2990954380; 1.279 5.27612561266 639.8972863140; 0.932 5.56741549127 647.0108333148; 0.756 6.17716234487 191.9584544356; 0.760 0.69475544472 302.1647756550; 1.038 0.23516951637 440.8252848776; 1.007 3.14159265359 0.0000000000; 0.549 4.87733288264 88.8656802170; 0.504 4.77955496203 419.4846438752; 0.346 4.31847547394 853.1963817520; 0.392 5.69922389094 654.1243803156 ]; .. latitude series b0 = [ 4330678.040 3.60284428399 213.2990954380; 240348.303 2.85238489390 426.5981908760; 84745.939 0.00000000000 0.0000000000; 30863.357 3.48441504465 220.4126424388; 34116.063 0.57297307844 206.1855484372; 14734.070 2.11846597870 639.8972863140; 9916.668 5.79003189405 419.4846438752; 6993.564 4.73604689179 7.1135470008; 4807.587 5.43305315602 316.3918696566; 4788.392 4.96512927420 110.2063212194; 3432.125 2.73255752123 433.7117378768; 1506.129 6.01304536144 103.0927742186; 1060.298 5.63099292414 529.6909650946; 969.071 5.20434966103 632.7837393132; 942.050 1.39646678088 853.1963817520; 707.645 3.80302329547 323.5054166574; 552.313 5.13149109045 202.2533951741; 399.675 3.35891413961 227.5261894396; 316.063 1.99716764199 647.0108333148; 319.380 3.62571550980 209.3669421749; 284.494 4.88648481625 224.3447957019; 314.225 0.46510272410 217.2312487011; 236.442 2.13887472281 11.0457002639; 215.354 5.94982610103 846.0828347512; 208.522 2.12003893769 415.5524906121; 178.958 2.95361514672 63.7358983034; 207.213 0.73021462851 199.0720014364; 139.140 1.99821990940 735.8765135318; 134.884 5.24500819605 742.9900605326; 140.585 0.64417620299 490.3340891794; 121.669 3.11537140876 522.5774180938; 139.240 4.59535168021 14.2270940016; 115.524 3.10891547171 216.4804891757; 114.218 0.96261442133 210.1177017003; 96.376 4.48164339766 117.3198682202; 80.593 1.31692750150 277.0349937414; 72.952 3.05988482370 536.8045120954; 69.261 4.92378633635 309.2783226558; 74.302 2.89376539620 149.5631971346; 68.040 2.18002263974 351.8165923087; 61.734 0.67728106562 1066.4954771900; 56.598 2.60963391288 440.8252848776; 48.864 5.78725874107 95.9792272178; 48.243 2.18211837430 74.7815985673; 38.304 5.29151303843 1059.3819301892; 36.323 1.63348365121 628.8515860501; 35.055 1.71279210041 1052.2683831884; 34.270 2.45740470599 422.6660376129; 34.313 5.97994514798 412.3710968744; 33.787 1.14073392951 949.1756089698; 31.633 4.14722153007 437.6438911399; 36.833 6.27769966148 1162.4747044078; 26.980 1.27154816810 860.3099287528; 23.516 2.74936525342 838.9692877504; 23.460 0.98962849901 210.8514148832; 23.600 4.11386961467 3.9321532631; 23.631 3.07427204313 215.7467759928; 20.813 3.51084686918 330.6189636582; 19.509 2.81857577372 127.4717966068; 17.103 3.89784279922 214.2623032845; 17.635 6.19715516746 703.6331846174; 17.824 2.28524493886 388.4651552382; 20.935 0.14356167048 430.5303441391; 16.551 1.66649120724 38.1330356378; 19.100 2.97699096081 137.0330241624; 15.517 4.54798410406 956.2891559706; 17.065 0.16611115812 212.3358875915; 14.169 0.48937283445 213.3472795478; 19.027 6.27326062836 423.4167971383 ]; b1 = [ 397554.998 5.33289992556 213.2990954380; 49478.641 3.14159265359 0.0000000000; 18571.607 6.09919206378 426.5981908760; 14800.587 2.30586060520 206.1855484372; 9643.981 1.69674660120 220.4126424388; 3757.161 1.25429514018 419.4846438752; 2716.647 5.91166664787 639.8972863140; 1455.309 0.85161616532 433.7117378768; 1290.595 2.91770857090 7.1135470008; 852.630 0.43572078997 316.3918696566; 284.386 1.61881754773 227.5261894396; 292.185 5.31574251270 853.1963817520; 275.090 3.88864137336 103.0927742186; 297.726 0.91909206723 632.7837393132; 172.359 0.05215146556 647.0108333148; 127.731 1.20711452525 529.6909650946; 166.237 2.44351613165 199.0720014364; 158.220 5.20850125766 110.2063212194; 109.839 2.45695551627 217.2312487011; 81.759 2.75839171353 210.1177017003; 81.010 2.86038377187 14.2270940016; 68.658 1.65537623146 202.2533951741; 59.281 1.82410768234 323.5054166574; 65.161 1.25527521313 216.4804891757; 61.024 1.25273412095 209.3669421749; 46.386 0.81534705304 440.8252848776; 36.163 1.81851057689 224.3447957019; 34.041 2.83971297997 117.3198682202; 32.164 1.18676132343 846.0828347512; 33.114 1.30557080010 412.3710968744; 27.282 4.64744847591 1066.4954771900; 22.805 4.12923703368 415.5524906121; 27.128 4.44228739187 11.0457002639; 18.100 5.56392353608 860.3099287528; 20.851 1.40999273740 309.2783226558; 14.947 1.34302610607 95.9792272178; 15.316 1.22393617996 63.7358983034; 14.601 1.00753704970 536.8045120954; 12.842 2.27059911053 742.9900605326; 12.832 4.88898877901 522.5774180938; 13.137 2.45991904379 490.3340891794; 11.883 1.87308666696 423.4167971383; 13.027 3.21731634178 277.0349937414; 9.946 3.11650057543 625.6701923124; 12.710 0.29501589197 422.6660376129; 9.644 1.74586356703 330.6189636582; 8.079 2.41931187953 430.5303441391; 8.245 4.68121931659 215.7467759928; 8.958 0.46482448501 429.7795846137; 6.547 3.01351967549 949.1756089698; 7.251 5.97098186912 149.5631971346; 6.056 1.49115011100 234.6397364404; 5.791 5.36720639912 735.8765135318; 5.994 0.02442871989 654.1243803156; 6.647 3.90879134581 351.8165923087; 6.824 1.52456408861 437.6438911399; 5.134 3.81149834833 74.7815985673; 3.959 5.63505813057 210.8514148832; 3.811 2.63992803111 3.1813937377; 3.643 1.73267151007 1059.3819301892; 3.554 4.98621474362 3.9321532631; 4.568 4.33599514584 628.8515860501; 3.145 2.51404811765 1162.4747044078; 3.522 1.16093567319 223.5940361765; 2.933 2.06057834252 956.2891559706 ]; b2 = [ 20629.977 0.50482422817 213.2990954380; 3719.555 3.99833475829 206.1855484372; 1627.158 6.18189939500 220.4126424388; 1346.067 0.00000000000 0.0000000000; 705.842 3.03914308836 419.4846438752; 365.042 5.09928680706 426.5981908760; 329.632 5.27899210039 433.7117378768; 219.335 3.82841533795 639.8972863140; 139.393 1.04272623499 7.1135470008; 103.980 6.15730992966 227.5261894396; 92.961 1.97994412845 316.3918696566; 71.242 4.14754353431 199.0720014364; 51.927 2.88364833898 632.7837393132; 48.961 4.43390206741 647.0108333148; 41.373 3.15927770079 853.1963817520; 28.602 4.52978327558 210.1177017003; 23.969 1.11595912146 14.2270940016; 20.511 4.35095844197 217.2312487011; 19.532 5.30779711223 440.8252848776; 18.263 0.85391476786 110.2063212194; 15.742 4.25767226302 103.0927742186; 16.840 5.68112084135 216.4804891757; 13.613 2.99904334066 412.3710968744; 11.567 2.52679928410 529.6909650946; 7.963 3.31512423920 202.2533951741; 6.599 0.28766025146 323.5054166574; 6.312 1.16121321336 117.3198682202; 5.891 3.58260177246 309.2783226558; 6.648 5.55714129949 209.3669421749; 5.590 2.47783944511 1066.4954771900; 6.192 3.61231886519 860.3099287528; 4.231 3.02212363572 846.0828347512; 3.612 4.79935735435 625.6701923124; 3.398 3.76732731354 423.4167971383; 3.387 6.04222745633 234.6397364404; 2.578 5.63610668746 735.8765135318; 2.831 4.81642822334 429.7795846137; 2.817 4.47516563908 654.1243803156; 2.573 0.22467245054 522.5774180938; 2.610 3.29126967191 95.9792272178; 2.419 0.02986335489 415.5524906121; 2.112 4.55964179603 422.6660376129; 2.304 6.25081073546 330.6189636582; 1.758 5.53430456858 536.8045120954; 1.814 5.05675881426 277.0349937414; 1.550 5.60375604692 223.5940361765; 1.457 4.47767649852 430.5303441391; 1.607 5.53599550100 224.3447957019; 1.172 4.71017775994 203.0041546995; 1.231 0.25115931880 3.9321532631; 1.105 1.01595427676 21.3406410024; 0.868 4.84623483952 949.1756089698; 0.939 1.35429452093 742.9900605326; 0.693 6.03599130692 124.4334152210; 0.712 4.45550701473 191.9584544356; 0.690 5.44243765037 437.6438911399; 0.810 0.46198177342 515.4638710930; 0.694 5.23748122403 447.9388318784; 0.604 2.95749705544 88.8656802170 ]; b3 = [ 666.252 1.99006340181 213.2990954380; 632.350 5.69778316807 206.1855484372; 398.051 0.00000000000 0.0000000000; 187.838 4.33779804809 220.4126424388; 91.884 4.84104208217 419.4846438752; 42.369 2.38073239056 426.5981908760; 51.548 3.42149490328 433.7117378768; 25.661 4.40167213109 227.5261894396; 20.551 5.85313509872 199.0720014364; 18.081 1.99321433229 639.8972863140; 10.874 5.37344546547 7.1135470008; 9.590 2.54901825866 647.0108333148; 7.085 3.45518372721 316.3918696566; 6.002 4.80055225135 632.7837393132; 5.778 0.01680378777 210.1177017003; 4.881 5.63719730884 14.2270940016; 4.501 1.22424419010 853.1963817520; 5.542 3.51756747774 440.8252848776; 3.548 4.71299370890 412.3710968744; 2.851 0.62679207578 103.0927742186; 2.173 3.71982274459 216.4804891757; 1.991 6.10867071657 217.2312487011; 1.435 1.69177141453 860.3099287528; 1.217 4.30778838827 234.6397364404; 1.157 5.75027789902 309.2783226558; 0.795 5.69026441157 117.3198682202; 0.733 0.59842720676 1066.4954771900; 0.713 0.21700311697 625.6701923124; 0.773 5.48361981990 202.2533951741; 0.897 2.65577866867 654.1243803156; 0.509 2.86079833766 429.7795846137; 0.462 4.17742567173 529.6909650946; 0.390 6.11288036049 191.9584544356; 0.505 4.51905764563 323.5054166574; 0.379 3.74436004151 223.5940361765; 0.332 5.49370890570 21.3406410024; 0.377 5.25624813434 95.9792272178; 0.384 4.48187414769 330.6189636582; 0.367 5.03190929680 846.0828347512; 0.281 1.14133888637 735.8765135318; 0.245 5.81618253250 423.4167971383; 0.241 1.70335120180 522.5774180938; 0.258 3.69110118716 447.9388318784 ]; b4 = [ 80.384 1.11918414679 206.1855484372; 31.660 3.12218745098 213.2990954380; 17.143 2.48073200414 220.4126424388; 11.844 3.14159265359 0.0000000000; 9.005 0.38441424927 419.4846438752; 6.164 1.56186379537 433.7117378768; 4.660 1.28235639570 199.0720014364; 4.775 2.63498295487 227.5261894396; 1.487 1.43096671616 426.5981908760; 1.424 0.66988083613 647.0108333148; 1.075 6.18092274059 639.8972863140; 1.145 1.72041928134 440.8252848776; 0.682 3.84841098180 14.2270940016; 0.655 3.49486258327 7.1135470008; 0.456 0.47338193402 632.7837393132; 0.509 0.31432285584 412.3710968744; 0.343 5.86413875355 853.1963817520; 0.270 2.50125594913 234.6397364404; 0.197 5.39156324804 316.3918696566; 0.236 2.11084590211 210.1177017003 ]; b5 = [ 7.895 2.81927558645 206.1855484372; 1.014 0.51187210270 220.4126424388; 0.772 2.99484124049 199.0720014364; 0.967 3.14159265359 0.0000000000; 0.583 5.96456944075 433.7117378768; 0.588 0.78008666397 227.5261894396 ]; .. radius vector r0 = [ 955758135.801 0.00000000000 0.0000000000; 52921382.465 2.39226219733 213.2990954380; 1873679.934 5.23549605091 206.1855484372; 1464663.959 1.64763045468 426.5981908760; 821891.059 5.93520025371 316.3918696566; 547506.899 5.01532628454 103.0927742186; 371684.449 2.27114833428 220.4126424388; 361778.433 3.13904303264 7.1135470008; 140617.548 5.70406652991 632.7837393132; 108974.737 3.29313595577 110.2063212194; 69007.015 5.94099622447 419.4846438752; 61053.350 0.94037761156 639.8972863140; 48913.044 1.55733388472 202.2533951741; 34143.794 0.19518550682 277.0349937414; 32401.718 5.47084606947 949.1756089698; 20936.573 0.46349163993 735.8765135318; 20839.118 1.52102590640 433.7117378768; 20746.678 5.33255667599 199.0720014364; 15298.457 3.05943652881 529.6909650946; 14296.479 2.60433537909 323.5054166574; 11993.314 5.98051421881 846.0828347512; 11380.261 1.73105746566 522.5774180938; 12884.128 1.64892310393 138.5174968707; 7752.769 5.85191318903 95.9792272178; 9796.061 5.20475863996 1265.5674786264; 6465.967 0.17733160145 1052.2683831884; 6770.621 3.00433479284 14.2270940016; 5850.443 1.45519636076 415.5524906121; 5307.481 0.59737534050 63.7358983034; 4695.746 2.14919036956 227.5261894396; 4043.988 1.64010323863 209.3669421749; 3688.132 0.78016133170 412.3710968744; 3376.457 3.69528478828 224.3447957019; 2885.348 1.38764077631 838.9692877504; 2976.033 5.68467931117 210.1177017003; 3419.551 4.94549148887 1581.9593482830; 3460.943 1.85088802878 175.1660598002; 3400.616 0.55386747515 350.3321196004; 2507.630 3.53851863255 742.9900605326; 2448.325 6.18412386316 1368.6602528450; 2406.138 2.96559220267 117.3198682202; 2881.181 0.17960757891 853.1963817520; 2173.959 0.01508587396 340.7708920448; 2024.483 5.05411271271 11.0457002639; 1740.254 2.34657043464 309.2783226558; 1861.397 5.93361638244 625.6701923124; 1888.436 0.02968443389 3.9321532631; 1610.859 1.17302463549 74.7815985673; 1462.631 1.92588134017 216.4804891757; 1474.547 5.67670461130 203.7378678824; 1395.109 5.93669404929 127.4717966068; 1781.165 0.76314388077 217.2312487011; 1817.186 5.77713225779 490.3340891794; 1472.392 1.40064915651 137.0330241624; 1304.089 0.77235613966 647.0108333148; 1149.773 5.74021249703 1162.4747044078; 1126.667 4.46707803791 265.9892934775; 1277.489 2.98412586423 1059.3819301892; 1207.053 0.75285933160 351.8165923087; 1071.399 1.13567265104 1155.3611574070; 1020.922 5.91233512844 1685.0521225016; 1315.042 5.11202572637 211.8146227297; 1295.553 4.69184139933 1898.3512179396; 1099.037 1.81765118601 149.5631971346; 998.462 2.63131596867 200.7689224658; 985.869 2.25992849742 956.2891559706; 932.434 3.66980793184 554.0699874828; 664.481 0.60297724821 728.7629665310; 659.850 4.66635439533 195.1398481733; 617.740 5.62092000007 942.0620619690; 626.382 5.94208232590 1478.8665740644; 482.230 1.84070179496 479.2883889155; 487.689 2.79373616806 3.1813937377; 470.086 0.83847755040 1471.7530270636; 451.817 5.64468459871 2001.4439921582; 553.128 3.41088600844 269.9214467406; 534.397 1.26443331367 275.5505210331; 472.572 1.88198584660 515.4638710930; 405.434 1.64001413521 536.8045120954; 517.196 4.44310450526 2214.7430875962; 452.848 3.00349117198 302.1647756550; 494.340 2.28626675074 278.5194664497; 489.825 5.80631420383 191.2076949102; 427.459 0.05741344372 284.1485407422; 339.763 1.40198657693 440.8252848776; 340.627 0.89091104306 628.8515860501; 385.974 1.99700402508 1272.6810256272; 288.298 1.12160250272 422.6660376129; 294.444 0.42577061903 312.1990839626 ]; r1 = [ 6182981.282 0.25843515034 213.2990954380; 506577.574 0.71114650941 206.1855484372; 341394.136 5.79635773960 426.5981908760; 188491.375 0.47215719444 220.4126424388; 186261.540 3.14159265359 0.0000000000; 143891.176 1.40744864239 7.1135470008; 49621.111 6.01744469580 103.0927742186; 20928.189 5.09245654470 639.8972863140; 19952.612 1.17560125007 419.4846438752; 18839.639 1.60819563173 110.2063212194; 12892.827 5.94330258435 433.7117378768; 13876.565 0.75886204364 199.0720014364; 5396.699 1.28852405908 14.2270940016; 4869.308 0.86793894213 323.5054166574; 4247.455 0.39299384543 227.5261894396; 3252.084 1.25853470491 95.9792272178; 2856.006 2.16731405366 735.8765135318; 2909.411 4.60679154788 202.2533951741; 3081.408 3.43662557418 522.5774180938; 1987.689 2.45054204795 412.3710968744; 1941.309 6.02393385142 209.3669421749; 1581.446 1.29191789712 210.1177017003; 1339.511 4.30801821806 853.1963817520; 1315.590 1.25296446023 117.3198682202; 1203.085 1.86654673794 316.3918696566; 1091.088 0.07527246854 216.4804891757; 954.403 5.15173410519 647.0108333148; 966.012 0.47991379141 632.7837393132; 881.827 1.88471724478 1052.2683831884; 874.215 1.40224683864 224.3447957019; 897.512 0.98343776092 529.6909650946; 784.866 3.06377517461 838.9692877504; 739.892 1.38225356694 625.6701923124; 612.961 3.03307306767 63.7358983034; 658.210 4.14362930980 309.2783226558; 649.600 1.72489486160 742.9900605326; 599.236 2.54924174765 217.2312487011; 502.886 2.12958819475 3.9321532631; 413.017 4.59334402271 415.5524906121; 356.117 2.30312127651 728.7629665310; 344.777 5.88787577835 440.8252848776; 395.004 0.53349091102 956.2891559706; 335.526 1.61614647174 1368.6602528450; 362.772 4.70691652867 302.1647756550; 321.611 0.97931764923 3.1813937377; 277.783 0.26007031431 195.1398481733; 291.173 2.83129427918 1155.3611574070; 264.971 2.42670902733 88.8656802170; 264.864 5.82860588985 149.5631971346; 316.777 3.58395655749 515.4638710930; 294.324 2.81632778983 11.0457002639; 244.864 1.04493438899 942.0620619690; 215.368 3.56535574833 490.3340891794; 264.047 1.28547685567 1059.3819301892; 246.245 0.90730313861 191.9584544356; 222.077 5.13193212050 269.9214467406; 194.973 4.56665009915 846.0828347512; 182.802 2.67913220473 127.4717966068; 181.645 4.93431600689 74.7815985673; 174.651 3.44560172182 137.0330241624; 165.515 5.99775895715 536.8045120954; 154.809 1.19720845085 265.9892934775; 169.743 4.63464467495 284.1485407422; 151.526 0.52928231044 330.6189636582; 152.461 5.43886711695 422.6660376129; 157.687 2.99559914619 340.7708920448; 140.630 2.02069760726 1045.1548361876; 139.834 1.35282959390 1685.0521225016; 140.977 1.27099900689 203.0041546995; 136.013 5.01678984678 351.8165923087; 153.391 0.26968607873 1272.6810256272; 129.476 1.14344730612 21.3406410024; 127.831 2.53876158952 1471.7530270636; 126.538 3.00310970076 277.0349937414; 100.277 3.61360169153 1066.4954771900; 103.169 0.38175114761 203.7378678824; 107.527 4.31870663477 210.8514148832 ]; r2 = [ 436902.464 4.78671673044 213.2990954380; 71922.760 2.50069994874 206.1855484372; 49766.792 4.97168150870 220.4126424388; 43220.894 3.86940443794 426.5981908760; 29645.554 5.96310264282 7.1135470008; 4141.650 4.10670940823 433.7117378768; 4720.909 2.47527992423 199.0720014364; 3789.370 3.09771025067 639.8972863140; 2963.990 1.37206248846 103.0927742186; 2556.363 2.85065721526 419.4846438752; 2208.457 6.27588858707 110.2063212194; 2187.621 5.85545832218 14.2270940016; 1956.896 4.92448618045 227.5261894396; 2326.801 0.00000000000 0.0000000000; 923.840 5.46392422737 323.5054166574; 705.936 2.97081280098 95.9792272178; 546.115 4.12854181522 412.3710968744; 373.838 5.83435991809 117.3198682202; 360.882 3.27703082368 647.0108333148; 356.350 3.19152043942 210.1177017003; 390.627 4.48106176893 216.4804891757; 431.485 5.17825414612 522.5774180938; 325.598 2.26867601656 853.1963817520; 405.018 4.17294157872 209.3669421749; 204.494 0.08774848590 202.2533951741; 206.854 4.02188336738 735.8765135318; 178.474 4.09716541453 440.8252848776; 180.143 3.59704903955 632.7837393132; 153.656 3.13470530382 625.6701923124; 147.779 0.13614300541 302.1647756550; 123.189 4.18895309647 88.8656802170; 133.076 2.59350469420 191.9584544356; 100.367 5.46056190585 3.1813937377; 131.975 5.93293968941 309.2783226558; 97.235 4.01832604356 728.7629665310; 110.709 4.77853798276 838.9692877504; 119.053 5.55385105975 224.3447957019; 93.852 4.38395529912 217.2312487011; 108.701 5.29310899841 515.4638710930; 78.609 5.72525447528 21.3406410024; 81.468 5.10897365253 956.2891559706; 96.412 6.25859229567 742.9900605326; 69.228 4.04901237761 3.9321532631; 65.168 3.77713343518 1052.2683831884; 64.088 5.81235002453 529.6909650946; 62.541 2.18445116349 195.1398481733; 56.987 3.14666549033 203.0041546995; 55.979 4.84108422860 234.6397364404; 52.940 5.07780548444 330.6189636582; 50.635 2.77318570728 942.0620619690; 41.649 4.79014211005 63.7358983034; 44.858 0.56460613593 269.9214467406; 41.357 3.73496404402 316.3918696566; 52.847 3.92623831484 949.1756089698; 38.398 3.73966157784 1045.1548361876; 37.583 4.18924633757 536.8045120954; 35.285 2.90795856092 284.1485407422; 33.576 3.80465978802 149.5631971346; 41.073 4.57870454147 1155.3611574070; 30.412 2.48140171991 860.3099287528; 31.373 4.84075951849 1272.6810256272; 30.218 4.35186294470 405.2575498736; 39.430 3.50858482049 422.6660376129; 29.658 1.58886982096 1066.4954771900; 35.202 5.94478241578 1059.3819301892 ]; r3 = [ 20315.005 3.02186626038 213.2990954380; 8923.581 3.19144205755 220.4126424388; 6908.677 4.35174889353 206.1855484372; 4087.129 4.22406927376 7.1135470008; 3879.041 2.01056445995 426.5981908760; 1070.788 4.20360341236 199.0720014364; 907.332 2.28344368029 433.7117378768; 606.121 3.17458570534 227.5261894396; 596.639 4.13455753351 14.2270940016; 483.181 1.17345973258 639.8972863140; 393.174 0.00000000000 0.0000000000; 229.472 4.69838526383 419.4846438752; 188.250 4.59003889007 110.2063212194; 149.508 3.20199444400 103.0927742186; 121.442 3.76831374104 323.5054166574; 101.215 5.81884137755 412.3710968744; 102.146 4.70974422803 95.9792272178; 93.078 1.43531270909 647.0108333148; 72.601 4.15395598507 117.3198682202; 84.347 2.63462379693 216.4804891757; 62.198 2.31239345505 440.8252848776; 45.145 4.37317047297 191.9584544356; 49.536 2.38854232908 209.3669421749; 54.829 0.30526468471 853.1963817520; 40.498 1.83836569765 302.1647756550; 38.089 5.94455115525 88.8656802170; 32.243 4.01146349387 21.3406410024; 40.671 0.68845183210 522.5774180938; 28.209 5.77193013961 210.1177017003; 24.976 3.06249709014 234.6397364404; 20.824 4.92570695678 625.6701923124; 25.070 0.73137425284 515.4638710930; 17.485 5.73135068691 728.7629665310; 18.009 1.45593152612 309.2783226558; 16.927 3.52771580455 3.1813937377; 13.437 3.36479898106 330.6189636582; 11.090 3.37212682914 224.3447957019; 11.082 3.41719974793 956.2891559706; 9.978 1.58791582772 202.2533951741; 11.551 5.99093726182 735.8765135318; 10.500 6.06911092266 405.2575498736; 9.144 2.93557421591 124.4334152210; 8.737 4.65432480769 632.7837393132; 10.023 0.58247011625 860.3099287528; 7.482 4.50669216436 942.0620619690; 10.091 0.28268774007 838.9692877504; 9.243 2.57034547708 223.5940361765; 8.652 1.75808100881 429.7795846137; 7.564 1.45635107202 654.1243803156 ]; r4 = [ 1202.050 1.41499446465 220.4126424388; 707.796 1.16153570102 213.2990954380; 516.121 6.23973568330 206.1855484372; 426.664 2.46924890293 7.1135470008; 267.736 0.18659206741 426.5981908760; 170.171 5.95926972384 199.0720014364; 145.113 1.44211060143 227.5261894396; 150.339 0.47970167140 433.7117378768; 121.033 2.40527320817 14.2270940016; 47.332 5.56857488676 639.8972863140; 15.745 2.90112466278 110.2063212194; 16.668 0.52920774279 440.8252848776; 18.954 5.85626429118 647.0108333148; 14.074 1.30343550656 412.3710968744; 12.708 2.09349305926 323.5054166574; 14.724 0.29905316786 419.4846438752; 11.133 2.46304825990 117.3198682202; 11.320 0.21785507019 95.9792272178; 9.233 2.28127318068 21.3406410024; 9.246 1.56496312830 88.8656802170; 8.970 0.68301278041 216.4804891757; 7.674 3.59367715368 302.1647756550; 7.823 4.48688804175 853.1963817520; 8.360 1.27239488455 234.6397364404; 9.552 3.14159265359 0.0000000000; 4.834 2.58836294602 515.4638710930; 6.059 5.16774448740 103.0927742186; 4.410 0.02211643085 191.9584544356; 4.364 1.59622746023 330.6189636582; 3.676 3.29899839673 210.1177017003; 4.364 5.97349927933 654.1243803156; 4.447 4.97415112184 860.3099287528; 3.220 2.72684237392 522.5774180938; 4.005 1.59858435636 405.2575498736; 3.099 0.75235436533 209.3669421749; 2.464 1.19167306488 124.4334152210; 3.088 1.32258934286 728.7629665310; 2.220 3.28087994088 203.0041546995; 2.127 6.14648095022 429.7795846137; 2.110 0.75462855247 295.0512286542; 2.020 3.89394929749 1066.4954771900; 2.248 0.49319150178 447.9388318784; 2.180 0.72761059998 625.6701923124; 1.809 0.09057839517 942.0620619690; 1.672 1.39635398184 224.3447957019; 1.641 3.02468307550 184.8449074348; 1.772 0.81879250825 223.5940361765 ]; r5 = [ 128.612 5.91282565136 220.4126424388; 32.273 0.69256228602 7.1135470008; 26.698 5.91428528629 227.5261894396; 19.923 0.67370653385 14.2270940016; 20.223 4.95136801768 433.7117378768; 13.537 1.45669521408 199.0720014364; 14.097 2.67074280191 206.1855484372; 13.364 4.58826996370 426.5981908760; 7.257 4.62966127155 213.2990954380; 4.876 3.61448275002 639.8972863140; 3.136 4.65661021909 191.9584544356; 2.917 0.48665273315 323.5054166574; 3.759 4.89624165044 440.8252848776; 3.303 4.07190859545 647.0108333148; 2.883 3.18003019204 419.4846438752; 2.338 3.69553554327 88.8656802170; 1.950 5.32729247780 302.1647756550 ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); sr5 = zzterm(r5, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4 + sr5*t5)/100000000; return [lambda, beta, delta]; endfunction function huranus(day) ## returns the heliocentric ecliptic latitude of Uranus given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## More terms than Meeus included as an experiment .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC l0 = [ ]; l1 = [ ]; l2 = [ ]; l3 = [ ]; l4 = [ ]; l5 = [ ]; .. latitude series b0 = [ ]; b1 = [ ]; b2 = [ ]; b3 = [ ]; b4 = [ ]; b5 = [ ]; .. radius vector r0 = [ ]; r1 = [ ]; r2 = [ ]; r3 = [ ]; r4 = [ ]; r5 = [ ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); sr5 = zzterm(r5, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4 + sr5*t5)/100000000; return [lambda, beta, delta]; endfunction function hneptune(day) ## returns the heliocentric ecliptic latitude of Neptune given ## the UT instant. Mean coordinates referred to equinox of date ## vector returned gives [lambda, beta, delta (au) ] ## See chapter 31 of Meeus for method - coefficients from NASA ADC ## More terms than Meeus included as an experiment .. Read in the coefficients for the series, these have been .. recopied from the VSOP82 series from NASA ADC l0 = [ ]; l1 = [ ]; l2 = [ ]; l3 = [ ]; l4 = [ ]; l5 = [ ]; .. latitude series b0 = [ ]; b1 = [ ]; b2 = [ ]; b3 = [ ]; b4 = [ ]; b5 = [ ]; .. radius vector r0 = [ ]; r1 = [ ]; r2 = [ ]; r3 = [ ]; r4 = [ ]; r5 = [ ]; t = day/365250; .. note, julian millenia not centuries t2 = t*t; t3 = t2*t; t4 = t2*t2; t5 = t4 * t; tpi = 2*pi(); sl0 = zzterm(l0, t); sl1 = zzterm(l1, t); sl2 = zzterm(l2, t); sl3 = zzterm(l3, t); sl4 = zzterm(l4, t); sl5 = zzterm(l5, t); lambda = (sl0 + sl1*t + sl2*t2 + sl3*t3 + sl4*t4 + sl5 * t5)/100000000; lambda = mod(360/tpi * lambda, 360); if lambda < 0; lambda = lambda + 360; endif; sb0 = zzterm(b0, t); sb1 = zzterm(b1, t); sb2 = zzterm(b2, t); sb3 = zzterm(b3, t); sb4 = zzterm(b4, t); sb5 = zzterm(b5, t); beta = (sb0 + sb1*t + sb2*t2 + sb3*t3 + sb4*t4 + sb5*t5)/100000000; beta = 360/tpi * beta; sr0 = zzterm(r0, t); sr1 = zzterm(r1, t); sr2 = zzterm(r2, t); sr3 = zzterm(r3, t); sr4 = zzterm(r4, t); sr5 = zzterm(r5, t); delta = (sr0 + sr1*t + sr2*t2 + sr3*t3 + sr4*t4 + sr5*t5)/100000000; return [lambda, beta, delta]; endfunction function moon(day) ## returns apparent geocentric equatorial coordinates of Moon given ## the UT instant day. .. find dphi and deps nut = nutation(day); .. find mean ecliptic coords of Moon meanmoon = gmoon(day); .. apparent longitude of Moon l = meanmoon[1] + nut[1]; b = meanmoon[2]; r = meanmoon[3]; e = obliquity(day) + nut[2]; alpha = datan2(dsin(l) *dcos(e) - dtan(b) * dsin(e), dcos(l)); if alpha < 0; alpha = alpha + 360; endif; delta = dasin(dsin(b) * dcos(e) + dcos(b) * dsin(e) * dsin(l)); return [alpha, delta, r]; endfunction function equatorial(day, ecliptic, e) ## given: the UT instant in day, a vector containing the ecliptic ## cooridnates and the obliquity of the ecliptic ## in e. The geocentric distance r is 'passed through'. ## returns: the equatorial coordinates as ra (degrees) and dec and r. ## note: obliquity argument allows choice of apparent, mean or ## J2000.0 equator. Function obiquity(day) returns the mean ## obliquity, add nutation[2] for apparent. l = ecliptic[1]; b = ecliptic[2]; r = ecliptic[3]; alpha = datan2(dsin(l) *dcos(e) - dtan(b) * dsin(e), dcos(l)); if alpha < 0; alpha = alpha + 360; endif; delta = dasin(dsin(b) * dcos(e) + dcos(b) * dsin(e) * dsin(l)); return [alpha, delta, r]; endfunction function apparent(day, ecliptic) ## given: the UT instant in day and a vector containing the mean ## ecliptic coordinates. ## returns: the apparent equatorial coordinates as ra (degrees) ## dec and r in a.u. referred to the equinox and true ecliptic ## of date nut = nutation(day); e = obliquity(day) + nut[2]; l = ecliptic[1] + nut[1]; b = ecliptic[2]; r = ecliptic[3]; alpha = datan2(dsin(l) *dcos(e) - dtan(b) * dsin(e), dcos(l)); if alpha < 0; alpha = alpha + 360; endif; delta = dasin(dsin(b) * dcos(e) + dcos(b) * dsin(e) * dsin(l)); return [alpha, delta, r]; endfunction function obliquity(day) ## returns: mean obliquity of ecliptic ## given: TD Instant as days since J2000.0 in 'day' t = day/ 36525; e = 84381.4 - 46.8150 * t - 0.00059 * t^2 + 0.001813 * t^3; e = e/3600; return e; endfunction function raltaz(day, station, equatorial) ## returns: altitude and azimuth of object ## given: 'day' - TD instant ## 'station' - vector containing geographical longitude, latitude ## height and air temperature and pressure of observer. ## (west, south negative, decimal degrees, metres asl, degrees C ## and millibars) ## note: for the Moon, use topocentric equatorial coords from ## tmoon - for other objects use 'apparent' coordinates ## refraction correction applied as in Meeus ch 15, p102 ## function altaz has no refraction correction glong = station[1]; glat = station[2]; height = station[3]; temp = station[4]; pressure = station[5]; ra = equatorial[1]; dec = equatorial[2]; lst = mod(gst(day) + glong, 360); if lst < 0; lst = lst + 360; endif; ha = mod(lst - ra, 360); if ha < 0; ha = ha + 360; endif; salt = dsin(dec)*dsin(glat) + dcos(dec)*dcos(glat)*dcos(ha); alt = dasin(salt); y = -dcos(dec)*dcos(glat)*dsin(ha); x = dsin(dec) - dsin(glat) * salt; az = datan2(y, x); .. refraction correction using Saemundsson's formula (see Meeus ch15) refrac = 1.02 / (dtan(alt + 10.3/(alt + 5.11))); refrac = refrac * (pressure /1010 * 283/(273 + temp)) / 60; alt = alt + refrac; return [az, alt, equatorial[3]]; endfunction function altaz(day, station, equatorial) ## returns: altitude and azimuth of object ## given: 'day' - TD instant ## 'station' - vector containing geographical longitude, latitude ## height and air temperature and pressure of observer. ## (west, south negative, decimal degrees, metres asl, degrees C ## and millibars) ## 'equatorial' - vector of equatorial coordinates of object ## note: for the Moon, use topocentric equatorial coords from ## tmoon - for other objects use 'apparent' coordinates ## This function does NOT correct for refraction - see raltaz glong = station[1]; glat = station[2]; height = station[3]; temp = station[4]; pressure = station[5]; ra = equatorial[1]; dec = equatorial[2]; lst = mod(gst(day) + glong, 360); if lst < 0; lst = lst + 360; endif; ha = mod(lst - ra, 360); if ha < 0; ha = ha + 360; endif; salt = dsin(dec)*dsin(glat) + dcos(dec)*dcos(glat)*dcos(ha); alt = dasin(salt); y = -dcos(dec)*dcos(glat)*dsin(ha); x = dsin(dec) - dsin(glat) * salt; az = datan2(y, x); return [az, alt, equatorial[3]]; endfunction function cart(sph) ## returns: cartesian coordinates of object in vector ## given: vector with [ra/long, dec/lat, r] of object ## note: assumes decimal degrees for angles r = sph[3]; l = sph[1]; b = sph[2]; x = r*dcos(b)*dcos(l); y = r*dcos(b)*dsin(l); z = r*dsin(b); return [x, y, z]; endfunction function sph(cart) ## returns: spherical coordinates of object in vector ## given: vector with cartesian coordinates of object ## note: returns decimal degrees for angles r = sqrt(cart.cart'); l = datan2(cart[2], cart[1]); b = datan(cart[3] / sqrt(cart[1]^2 + cart[2]^2)); return [l, b, r]; endfunction function sun(day) ## returns: apparent equatorial coordinates of the Sun (geocentric) ## given: TD instant in day ## note: coordinates returned as vector [ra (degs), dec, r (au)] ## based on Meeus truncation of VSOP87 - he claims 1" earth = hearth(day); nut = nutation(day); ob = obliquity(day) + nut[2]; b = -1 * earth[2]; r = earth[3]; abb = (-20.4898 / r)/3600; l = earth[1] + 180 + abb + nut[1]; out = equatorial(day, [l, b, r], ob); return out; endfunction function mean(day, ecliptic) ## given: the TD instant in day and a vector containing the mean ## ecliptic geocentric coordinates. ## returns: the mean equatorial coordinates as ra (degrees) ## dec and r in a.u. referred to the equinox and mean ecliptic ## of date ## note: no nutation or aberration correction at all nut = nutation(day); e = obliquity(day) ; l = ecliptic[1] ; b = ecliptic[2]; r = ecliptic[3]; alpha = datan2(dsin(l) *dcos(e) - dtan(b) * dsin(e), dcos(l)); if alpha < 0; alpha = alpha + 360; endif; delta = dasin(dsin(b) * dcos(e) + dcos(b) * dsin(e) * dsin(l)); return [alpha, delta, r]; endfunction function gmer(day) ## given: the TD instant in day ## returns: the geometric equatorial coordinates of Mercury ## note: no correction for light travel time, abberation or nutation p = cart(hmercury(day)); e = cart(hearth(day)); ec = sph(p - e); return mean(day, ec); endfunction function gjup(day) ## given: the TD instant in day ## returns: the geometric equatorial coordinates of Jupiter ## note: no correction for light travel time, abberation or nutation p = cart(hjupiter(day)); e = cart(hearth(day)); ec = sph(p - e); return mean(day, ec); endfunction function gven(day) ## given: the TD instant in day ## returns: the geometric equatorial coordinates of Venus ## note: no correction for light travel time, abberation or nutation p = cart(hvenus(day)); e = cart(hearth(day)); ec = sph(p - e); return mean(day, ec); endfunction function gmar(day) ## given: the TD instant in day ## returns: the geometric equatorial coordinates of Mars ## note: no correction for light travel time, abberation or nutation p = cart(hmars(day)); e = cart(hearth(day)); ec = sph(p - e); return mean(day, ec); endfunction function amar(day) ## given: the TD instant in day ## returns: the astrometric equatorial coordinates of Mars ## note: corrected for light travel time and aberration e = cart(hearth(day)); day1 = day; r1 = 0; lp = 0; repeat lp = lp + 1 p = cart(hmars(day1)); gr = p - e; r2 = r1; r1 = sqrt(gr.gr'); day1 = day - 0.0057755183 * r1; if abs(r1 - r2) < 10^-6 || lp == 6; break; endif; end; ec = sph(gr); return mean(day, ec); endfunction function venus(day) ## given: the TD instant in day ## returns: the apparent equatorial coordinates of Venus ## note: corrected for light travel time, aberration ## and nutation out = zzaplanet("hvenus", day); return out; endfunction function mercury(day) ## given: the TD instant in day ## returns: the apparent equatorial coordinates of Mercury ## note: corrected for light travel time, aberration ## and nutation out = zzaplanet("hmercury", day); return out; endfunction function jupiter(day) ## given: the TD instant in day ## returns: the apparent equatorial coordinates of Jupiter ## note: corrected for light travel time, aberration ## and nutation out = zzaplanet("hjupiter", day); return out; endfunction function mars(day) ## given: the TD instant in day ## returns: the apparent equatorial coordinates of Mars ## note: corrected for light travel time, aberration ## and nutation out = zzaplanet("hmars", day); return out; endfunction function saturn(day) ## given: the TD instant in day ## returns: the apparent equatorial coordinates of Saturn ## note: corrected for light travel time, aberration ## and nutation out = zzaplanet("hsaturn", day); return out; endfunction function zzaplanet(planet, day) ## given: the TD instant in day and heliocentric coordinate function ## in string 'planet' ## returns: the apparent equatorial coordinates of planet ## note: corrected for light travel time, aberration ## and nutation ## this function is not normally used interactively .. iterative light travel time calculation se = hearth(day); e = cart(se); day1 = day; r1 = 0; lp = 0; repeat lp = lp + 1; dummy = planet(day1); p = cart(dummy); gr = p - e; r2 = r1; r1 = sqrt(gr.gr'); day1 = day - 0.0057755183 * r1; if abs(r1 - r2) < 10^-6 || lp == 6; break; endif; end; ec = sph(gr); .. correction for aberration due to motion of Earth t = day/ 36525; t2 = t * t; ecc = 0.016708617 - 0.000042037*t - 0.0000001236 * t2; epi = 102.93735 + 1.71953 * t + 0.00046 *t2; sunl = se[1] + 180; k = -20.49552 /3600; dl = (k*dcos(sunl - ec[1]) + ecc*k*dcos(epi - ec[1])) / dcos(ec[2]); db = k*dsin(ec[2])* (sin(sunl - ec[1]) - ecc* dsin(epi - sunl)); ec[1] = ec[1] + dl; ec[2] = ec[2] +db; return apparent(day, ec); endfunction function topocentric(day, station, pos) ## returns: topcentric equatorial coordinates of body ## given: function 'pos' to obtain apparent coords of body ## TD instant in 'day' ## station - lat, long, height, temp, pressure of observing ## location .. get constants rho sin (phi') and rho cos(phi') .. Meeus ch 10 geo = pos(day); f1 = 0.99664719; phi = station[2]; h = station[3]; u = datan(f1 * dtan(phi)); rhosin = f1 * dsin(u) + h / 6378140 * dsin(phi); rhocos = dcos(u) + h/ 6378140 * dcos(phi); .. apply the differential formulas for topcentric ha and dec .. find if this is the Moon or not and get parallax if geo[3] > 1000; sinpar = 6378.14 / geo[3]; else; sinpar = 8.794 / (geo[3] * 3600); endif; .. find LST and hence geocentric hour angle lst = gst(day) + station[1]; ha = lst - geo[1]; gdec = geo[2]; a = dcos(gdec) * dsin(ha); b = dcos(gdec) * dcos(ha) - rhocos * sinpar; c = dsin(gdec) - rhosin * sinpar; q = sqrt(a*a + b*b + c*c); tha = datan2(a, b); tdec = dasin(c/q); ra = lst - tha; if ra < 0; ra = ra + 360; endif; r = q* geo[3]; return [ra, tdec, r ]; endfunction function brum() ## returns: a vector with the longitude, latitude, height, temperature ## and pressure typical to Birmingham UK ## note: modify temperature and pressure if accurate refraction ## adjustments needed. return [-1.91667, 52.5, 120, 10, 1012 ]; endfunction function settle() ## returns: a vector with the longitude, latitude, height, temperature ## and pressure typical to Settle, near Skipton, Yorkshire UK ## note: modify temperature and pressure if accurate refraction ## adjustments needed. return [-2 + 18/60, 54 + 5/60, 200, 10, 1012 ]; endfunction function reekie() ## returns: a vector with the longitude, latitude, height, temperature ## and pressure typical to Edinburgh, Scotland ## note: modify temperature and pressure if accurate refraction ## adjustments needed. return [-3 + 12/60, 55 + 57/60, 50, 10, 1012 ]; endfunction function tmoon(day, station) ## returns: topcentric equatorial coordinates of Moon ## given: TD instant in 'day' ## 'station' - vector lat, long, height (m), temp (C), ## and pressure (mB) of observing location return topocentric(day, station, "moon"); endfunction function librate(day) ## returns: selenographic longitude and latitude of the sub earth point ## position angle of the Moon's rotation axis ## given: TD instant in days since J2000.0 ## note: Geocentric librations as in Meeus Ch 51 .. arguments (put these in a separate function eventually) t = day/36525; t2 = t*t; t3 = t2*t; t4 = t2*t2; .. Mean longitude including light travel time and referred to equinox of date l1 = mod(218.3164591 + 481267.88134236 *t - 0.0013268 *t2 + t3/538841 - t4/65194000, 360); .. Mean elongation d = mod(297.8502042 + 445267.1115168 *t - 0.0016300 *t2 + t3/545868 - t4/113065000, 360); .. Sun's mean anomaly m = mod(357.5291092 + 35999.0502909 *t - 0.0001536 *t2 + t3/24490000, 360); .. Moon's mean anomaly m1 = mod(134.9634114 + 477198.8676313 *t + 0.0089970 *t2 + t3/69699 - t4/14712000, 360); ..Moon's argument of latitude (mean distance from ascending node) f = mod(93.2720993 + 483202.0175273 *t - 0.0034029 *t2 - t3/3526000 + t4/863310000, 360); if f < 0; f = f + 360; endif; ..Moon's longitude of mean ascending node om = mod(125.0445550 - 1934.1361849 * t + 0.0020762 * t2 + t3 / 467410 - t4 / 60616000, 360); if om < 0; om = om + 360; endif; ..eccentricity of Earth's orbit .. Eccentricity of earth's orbit round Sun (affects terms with M or 2M) e = 1 - 0.002516 * t - 0.0000074 *t2; .. Optical librations i = 1 + 32/60 + 32.7 /3600; mm = gmoon(day); am = amoon(day); lambda = mm[1]; beta = am[2]; w = lambda - om; y = dsin(w) * dcos(beta) * dcos(i) - dsin(beta) * dsin(i); x = dcos(w) * dcos(beta); a = datan2(y , x); la = a - f; ba = dasin(-dsin(w) * dcos(beta) * dsin(i) - dsin(beta)* dcos(i)); .. Physical libration calculation - see Meeus ch51 p 343 k1 = 119.75 + 131.849 * t; k2 = 72.56 + 20.186 * t; rho = -0.02752* dcos(m1); rho = rho - 0.02245 * dsin(f); rho = rho + 0.00684 * dcos(m1 - 2*f); rho = rho - 0.00293 * dcos(2*f); rho = rho - 0.00085 * dcos(2*f - 2*d); rho = rho - 0.00054 * dcos(m1 - 2 * d); rho = rho - 0.00020 * dsin(m1 + f); rho = rho - 0.00020 * dcos(m1 + 2*f); rho = rho - 0.00020 * dcos(m1 - f); rho = rho + 0.00014 * dcos(m1 + 2*f - 2*d); sigma = - 0.02816 * dsin(m1); sigma = sigma + 0.02244 * dcos(f); sigma = sigma - 0.00682 * dsin(m1 - 2*f); sigma = sigma - 0.00279 * dsin(2*f); sigma = sigma - 0.00083 * dsin(2*f - 2*d); sigma = sigma + 0.00069 * dsin(m1 - 2*d); sigma = sigma + 0.00040 * dcos(m1 + f); sigma = sigma - 0.00025 * dsin(2*m1); sigma = sigma - 0.00023 * dsin(m1 + 2*f); sigma = sigma + 0.00020 * dcos(m1 - f); sigma = sigma + 0.00019 * dsin(m1 - f); sigma = sigma + 0.00013 * dsin(m1 + 2*f - 2*d); sigma = sigma - 0.00010 * dcos(m1 - 3*f); tau = + 0.02520 * e * dsin(m); tau = tau + 0.00473 * dsin(2*m1 - 2*f); tau = tau - 0.00467 * dsin(m1); tau = tau + 0.00396 * dsin(k1); tau = tau + 0.00276 * dsin(2*m1 - 2*d); tau = tau + 0.00196 * dsin(om); tau = tau - 0.00183 * dcos(m1 - f); tau = tau + 0.00115 * dsin(m1 - 2*d); tau = tau - 0.00096 * dsin(m1 - d); tau = tau + 0.00046 * dsin(2*f - 2*d); tau = tau - 0.00039 * dsin(m1 - f); tau = tau - 0.00032 * dsin(m1 - m - d); tau = tau + 0.00027 * dsin(2*m1 - m - 2*d); tau = tau + 0.00023 * dsin(k2); tau = tau - 0.00014 * dsin(2*d); tau = tau + 0.00014 * dcos(2*m1 - 2*f); tau = tau - 0.00012 * dsin(m1 - 2*f); tau = tau - 0.00012 * dsin(2*m1); tau = tau + 0.00011 * dsin(2*m1 - 2*m - 2*d); laa = -tau + ( rho * dcos(a) + sigma * dsin(a)) * dtan(ba); baa = sigma * dcos(a) - rho * dsin(a); ..total libration l0 = la + laa; b0 = ba + baa; ..position angle of Moon's rotation axis measured from ..celestial north pole nut = nutation(day); epsilon = obliquity(day) + nut[2]; apparent = moon(day); alpha = apparent[1]; v = om + nut[1] + sigma / dsin(i); x = dsin(i + rho) * dsin(v); y = dsin(i + rho) * dcos(v) * dcos(epsilon) - dcos(i + rho) * dsin(epsilon); w2 = datan2(x, y); p = dasin(sqrt(x*x + y*y) * dcos(alpha - w2) / dcos(b0)); return [l0, b0, p]; endfunction