File: astro.e

package info (click to toggle)
euler 1.61.0-12
links: PTS, VCS
area: main
in suites: bookworm, bullseye, forky, sid, trixie
size: 5,164 kB
sloc: ansic: 24,761; sh: 8,314; makefile: 141; cpp: 47; php: 1
file content (4674 lines) | stat: -rw-r--r-- 151,431 bytes
parent folder | download | duplicates (8)
.. 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