/*
** Astrolog (Version 4.40) File: placalc2.c
**
** IMPORTANT NOTICE: The graphics database and chart display routines
** used in this program are Copyright (C) 1991-1995 by Walter D. Pullen
** (astara@u.washington.edu). Permission is granted to freely use and
** distribute these routines provided one doesn't sell, restrict, or
** profit from them in any way. Modification is allowed provided these
** notices remain with any altered or edited versions of the program.
**
** The main planetary calculation routines used in this program have
** been Copyrighted and the core of this program is basically a
** conversion to C of the routines created by James Neely as listed in
** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
** available from Matrix Software. The copyright gives us permission to
** use the routines for personal use but not to sell them or profit from
** them in any way.
**
** The PostScript code within the core graphics routines are programmed
** and Copyright (C) 1992-1993 by Brian D. Willoughby
** (brianw@sounds.wa.com). Conditions are identical to those above.
**
** The extended accurate ephemeris databases and formulas are from the
** calculation routines in the program "Placalc" and are programmed and
** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
** (alois@azur.ch). The use of that source code is subject to
** regulations made by Astrodienst Zurich, and the code is not in the
** public domain. This copyright notice must not be changed or removed
** by any user of this program.
**
** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
** X Window graphics initially programmed 10/23-29/1991.
** PostScript graphics initially programmed 11/29-30/1992.
** Last code change made 1/29/1995.
*/

#include "placalc.h"


#ifdef PLACALC
/*
** ---------------------------------------------------------------
** | Copyright Astrodienst AG and Alois Treindl, 1989,1991,1993  |
** | The use of this source code is subject to regulations made  |
** | by Astrodienst Zurich. The code is NOT in the public domain.|
** |                                                             |
** | This copyright notice must not be changed or removed        |
** | by any user of this program.                                |
** ---------------------------------------------------------------
**
** Important changes:
** 11-jun-93 revision 1.12: fixed error which affected Mercury between -2100
** and -3100 (it jumped wildly).
*/

/***********************************************************
** $Header$
**
** definition module for planetary elements
** and disturbation coefficients
** version HP-UX C  for new version with stored outer planets
** 31-jul-88
** by Alois Treindl 
**
** ---------------------------------------------------------------
** | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
** | The use of this source code is subject to regulations made |
** | by Astrodienst Zurich. The code is NOT in the public domain.|
** |                                                             |
** | This copyright notice must not be changed or removed        |
** | by any user of this program.                                |
** ---------------------------------------------------------------
**
***********************************************************/

/*
** In the elements degrees were kept as the units for the constants. This
** requires conversion to radians, when the actual calculations are performed.
** This approach is not the most efficient, but safer for development.
** Constant conversion could be done by writing all degree constants with
** value * DEGTORAD
*/

#define TIDAL_26  TRUE  /* decide wheter to use new or old lunar tidal
term; a consistent system of delta t must be 
used */
#define MOON_TEST_CORR FALSE  /* to include more lunar terms in longitude */

REAL8 ekld[4] = {23.452294, -46.845, -.0059, 0.00181};
/* ecliptic with epoch1900, Ekd(0..3) in basic */

struct eledata pd[MARS + 1] = {
{/*earth*/  1.00000023, 365.25636042, EPOCH1900,
99.696678,  .9856473354,  1.089,    0,  
101.220833, 6189.03,  1.63,   0.012,
0.01675104, -0.00004180,  -0.000000126,
0,    0,    0,    0,
0,    0,    0},
/*
** note 29 June 88 by Alois: G.M.Clemence, Astronomical Journal
** vol.53,p. 178 (1948) gives a correction to the perihel motion
** of -4.78" T, giving 6184.25 for the linear Term above. We have
** not yet applied this correction. It has been used in APAE 22,4
** on the motion of mars and does make an official impression.
*/
{/*moon*/ 0.0025955307, 27.321661,  EPOCH1900,
# if ! TIDAL_26 
/*
** values from Improved Lunar Ephemeris, corresponding to tidal
** term -22.44"/cy and  consistent with delta t ~ 29.949 T*T
*/
270.4341638,  13.176396526808121, -4.08,  0.0068,
# endif
# if TIDAL_26 
/*
** new values from Morrison 1979, with tidal term -26"/cy as
** stated in A.E. 1986 onwards, consistent with delta t ~ 44.3 T*T 
** correction: -1.54" + 2.33" T - 1.78" T*T
*/
270.4337361,  13.176396544528099, -5.86,  0.0068,
# endif
334.329556, 14648522.52,  -37.17,   -0.045,
0.054900489,  0,    0,
259.183275, -6962911.23,  7.48 ,    0.008,
5.145388889,  0,    0},
{/*mercury*/  .3870986, 87.969252,  EPOCH1900,
178.179078, 4.0923770233, 1.084,    0,
75.89969722,  5599.76,  1.061,    0,
0.20561421, .00002046,   -.000000030,
47.145944444, 4266.75,  .626,   0,
7.0028805555, 6.699,    -.066},
{/*venus*/  .72333162,  224.700726, EPOCH1900,
342.767053, 1.6021687039, 1.1148,   0,
130.16383333, 5068.93,  -3.515,   0,
0.00682069, -.00004774, .000000091,
75.7796472223,3239.46,  1.476,    0,
3.3936305555, 3.621,    .0035},
{/*mars*/ 1.5236914620, 686.9296097,  EPOCH1900,
/* These are the corrected elements by Ross */
293.74762778, .524071163814,  1.1184,   0,
334.21820278, 6626.73,  .4675,    -0.0043,
0.09331290, .000092064, -.000000077,
48.786441667, 2775.57,  -.005,    -0.0192,
1.85033333, -2.430,   .0454}
};

/*
** mimimum and maximum distances computed over 1000 years with plamimax,
** required for relative distances rgeo, where the distance is given
** as 100 when a planet is closest and as 0 when farthest from earth.
*/
REAL8 rmima[CALC_N][2] = {  
  { 0.98296342,  1.01704665},
  { 0.00238267,  0.00271861},
  { 0.54900496,  1.45169607},
  { 0.26411287,  1.73597885},
  { 0.37289847,  2.67626927},
  { 3.94877993,  6.45627627},
  { 7.99362824, 11.09276636},
  {17.28622633, 21.10714104},
  {28.81374786, 31.33507284},
  {28.67716748, 50.29208774},
  { 0.00,  0.00259553}, /* nodes don't get a real value */
  { 0.00,  0.00259553},
  { 7.36277475, 19.86585062}};


struct sdat _sd [SDNUM] = {
  114.50, 585.17493,
  109.856,  191.39977,
  148.031,  30.34583,
  284.716,  12.21794,
  114.508,  585.17656,
  -0.56,    359.99213,
  148.03, 30.34743,
  284.72, 12.2196,
  248.07, 1494.726615,
  359.44, 359.993595,
  109.86, 191.402867,
  148.02, 30.348930,
  114.503,  585.173715,
  359.444,  359.989285,
  148.021,  30.344620,
  284.716,  12.21669,
  148.0315, 30.34906264,
  284.7158, 12.22117085,
  220.1695,   4.284931111,
  291.8024,   2.184704167
};

REAL8 sa[SDNUM];

/*
** delta long = lampl * COS (lphase - arg) in seconds of arc
** delta rad  = rampl * COS (rphase - arg) in ninth place of log
** arg = j * sa (k) + i * ma (this planet)
** ma = mean anomaly
** sa = mean anomaly of disturbing planet, where this
** is taken from the aproximate value in sa[]
** For the COS (phase - arg) it is good enough to compute
** with 32 bit reals, because ampl and phase have only 
** four to five significant digits. 
** While saving constant space, it is costing execution time due
** to float/double conversions.
**
** In basic, all correction terms for sun, mercury, venus and mars
** were contained in one array K(0..142,0..6); Nk(N,0) contained
** the index of the first term of planet N and Nk(N,1) the number
** of terms for this planet. Here, we use a  0 in the first column
** kor.j to indicate the end of the table for a planet.
** K(*) was a basic INTEGER array, therefore the amplitudes and phases
** had to be expressed as
** K(i,2) = ampl. of longitude in 0.001 seconds of arc
** K(i,3) = phase of longitude in 0.01 degrees
** K(i,4) = ampl. of radius in 9th place of log
** K(i,5) = phase of radius in 0.01 degrees.
** Here we have converted the amplitude of long. to seconds of arc
** and the phases to degrees.
*/

struct kor ARR earthkor[] = {  /* 11-jul-88 all terms to 0.020" longitude */
  /*  j     i       lampl   lphase  rampl   rphase  k */
  -1, 1,  0.013,  243,  28, 335,  8,  /* mercury */
  -1, 3,  0.015,  357,  18, 267,  8,
  -1, 4,  0.023,  326,  5,  239,  8,
  -1, 0,  0.075,  296.6,  94, 205.0,  0,  /* Venus */
  -1, 1,  4.838,  299.10, 2359, 209.08, 0,  
  -1, 2,  0.074,  207.9,  69, 348.5,  0,
  -1, 3,  0.009,  249,  16, 330,  0,
  -2, 1,  .116, 148.90, 160,  58.40,  0,  
  -2, 2,  5.526,  148.31, 6842, 58.32,  0,  
  -2, 3,  2.497,  315.94, 869,  226.70, 0,  
  -2, 4,  0.044,  311.4,  52, 38.8, 0,
  -3, 2,  0.013,  176,  21, 90, 0,
  -3, 3,  .666, 177.71, 1045, 87.57,  0,  
  -3, 4,  1.559,  -14.75, 1497, 255.25, 0,  
  -3, 5,  1.024,  318.15, 194,  49.50,  0,  
  -3, 6,  0.017,  315,  19, 43, 0,
  -4, 4,  .210, 206.20, 376,  116.28, 0,  
  -4, 5,  .144, 195.40, 196,  105.20, 0,  
  -4, 6,  .152, -16.20, 94, 254.80, 0,  
  -5, 5,  0.084,  235.6,  163,  145.4,  0,
  -5, 6,  0.037,  221.8,  59, 132.2,  0,
  -5, 7,  .123, 195.30, 141,  105.40, 0,  
  -5, 8,  .154, -.40, 26, 270.00, 0,  
  -6, 6,  0.038,  264.1,  80, 174.3,  0,
  -6, 7,  0.014,  253,  25, 164,  0,
  -6, 8,  0.01, 230,  14, 135,  0,
  -6, 9,  0.014,  12, 12, 284,  0,
  -7, 7,  0.020,  294,  42, 203.5,  0,
  -7, 8,  0.006,  279,  12, 194,  0,
  -8, 8,  0.011,  322,  24, 234,  0,
  -8, 12, 0.042,  259.2,  44, 169.7,  0,
  -8, 14, 0.032,  48.8, 33, 138.7,  0,
  -9, 9,  0.006,  351,  13, 261,  0,
  1,  -1, .273, 217.70, 150,  127.70, 1,  /* mars */
  1,  0,  0.048,  260.3,  28, 347,  1,
  2,  -3, 0.041,  346,  52, 255.4,  1,
  2,  -2, 2.043,  343.89, 2057, 253.83, 1,
  2,  -1, 1.770,  200.40, 151,  295.00, 1,
  2,  0,  0.028,  148,  31, 234.3,  1,  
  3,  -3, .129, 294.20, 168,  203.50, 1,
  3,  -2, .425, -21.12, 215,  249.00, 1,
  4,  -4, 0.034,  71, 49, 339.7,  1,
  4,  -3, .500, 105.18, 478,  15.17,  1,
  4,  -2, .585, -25.94, 105,  65.90,  1,
  5,  -4, 0.085,  54.6, 107,  324.6,  1,
  5,  -3, .204, 100.80, 89, 11.00,  1,
  6,  -5, 0.020,  186,  30, 95.7, 1,
  6,  -4, .154, 227.40, 139,  137.30, 1,
  6,  -3, .101, 96.30,  27, 188.00, 1,
  7,  -5, 0.049,  176.5,  60, 86.2, 1,
  7,  -4, .106, 222.70, 38, 132.90, 1,
  8,  -5, 0.052,  348.9,  45, 259.7,  1,
  8,  -4, 0.021,  215.2,  8,  310,  1,
  8,  -6, 0.010,  307,  15, 217,  1,
  9,  -6, 0.028,  298,  34, 208.1,  1,
  9,  -5, 0.062,  346,  17, 257,  1,
  10, -6, 0.019,  111,  15, 23, 1,
  11, -7, 0.017,  59, 20, 330,  1,
  11, -6, 0.044,  105.9,  9,  21, 1,
  13, -8, 0.013,  184,  15, 94, 1,
  13, -7, 0.045,  227.8,  5,  143,  1,
  15, -9, 0.021,  309,  22, 220,  1,
  17, -9, 0.026,  113,  0,  0,  1,
  1,  -2, .163, 198.60, 208,  112.00, 2,  /* jupiter */
  1,  -1, 7.208,  179.53, 7067, 89.55,  2,
  1,  0,  2.600,  263.22, 244,  -21.40, 2,
  1,  1,  0.073,  276.3,  80, 6.5,  2,
  2,  -3, 0.069,  80.8, 103,  350.5,  2,
  2,  -2, 2.731,  87.15,  4026, -2.89,  2,
  2,  -1, 1.610,  109.49, 1459, 19.47,  2,
  2,  0,  0.073,  252.6,  8,  263,  2,
  3,  -3, .164, 170.50, 281,  81.20,  2,
  3,  -2, .556, 82.65,  803,  -7.44,  2,
  3,  -1, .210, 98.50,  174,  8.60, 2,
  4,  -4, 0.016,  259,  29, 170,  2,
  4,  -3, 0.044,  168.2,  74, 79.9, 2,
  4,  -2, 0.080,  77.7, 113,  347.7,  2,
  4,  -1, 0.023,  93, 17, 3,  2,
  5,  -2, 0.009,  71, 14, 343,  2,
  1,  -2, 0.011,  105,  15, 11, 3,  /* saturn */
  1,  -1, .419, 100.58, 429,  10.60,  3,
  1,  0,  .320, 269.46, 8,  -7.00,  3,
  2,  -2, .108, 290.60, 162,  200.60, 3,
  2,  -1, .112, 293.60, 112,  203.10, 3,
  3,  -2, 0.021,  289,  32, 200.1,  3,
  3,  -1, 0.017,  291,  17, 201,  3,
  ENDMARK
};

struct kor ARR mercurykor[] = {
  1,  -1, .711, 35.47,  491,  305.28, 4,
  2,  -3, .552, 161.15, 712,  71.12,  4,
  2,  -2, 2.100,  161.15, 2370, 71.19,  4,
  2,  -1, 3.724,  160.69, 899,  70.49,  4,
  2,  0,  .729, 159.76, 763,  250.00, 4,
  3,  -3, .431, 105.37, 541,  15.53,  4,
  3,  -2, 1.329,  104.78, 1157, 14.84,  4,
  3,  -1, .539, 278.95, 14, 282.00, 4,
  4,  -2, .484, 226.40, 234,  136.02, 4,
  5,  -4, .685, -10.43, 849,  259.51, 4,
  5,  -3, 2.810,  -10.14, 2954, 259.92, 4,
  5,  -2, 7.356,  -12.22, 282,  255.43, 4,
  5,  -1, 1.471,  -12.30, 1550, 77.75,  4,
  5,  0,  .375, -12.29, 472,  77.70,  4,
  2,  -1, .443, 218.48, 256,  128.36, 5,
  4,  -2, .374, 151.81, 397,  61.63,  5,
  4,  -1, .808, 145.93, 13, 35.00,  5,
  1,  -1, .697, 181.07, 708,  91.38,  6,
  1,  0,  .574, 236.72, 75, 265.40, 6,
  2,  -2, .938, 36.98,  1185, 306.97, 6,
  2,  -1, 3.275,  37.00,  3268, 306.99, 6,
  2,  0,  .499, 31.91,  371,  126.90, 6,
  3,  -1, .353, 25.84,  347,  295.76, 6,
  2,  -1, .380, 239.87, 0,  0,  7,
  ENDMARK
};

struct kor ARR venuskor[] = {
  -1, 2,  .264, -19.20, 175,  251.10, 8,
  -2, 5,  .361, 167.68, 55, 77.20,  8,
  1,  -1, 4.889,  119.11, 2246, 29.11,  9,
  2,  -2, 11.261, 148.23, 9772, 58.21,  9,
  3,  -3, 7.128,  -2.57,  8271, 267.42, 9,
  3,  -2, 3.446,  135.91, 737,  47.37,  9,
  4,  -4, 1.034,  26.54,  1426, 296.49, 9,
  4,  -3, .677, 165.32, 445,  75.70,  9,
  5,  -5, .330, 56.88,  510,  -33.36, 9,
  5,  -4, 1.575,  193.93, 1572, 104.21, 9,
  5,  -3, 1.439,  138.08, 162,  229.90, 9,
  6,  -6, .143, 84.40,  236,  -5.80,  9,
  6,  -5, .205, 44.20,  256,  314.20, 9,
  6,  -4, .176, 164.30, 70, 75.70,  9,
  8,  -5, .231, 180.00, 25, 75.00,  9,
  3,  -2, .673, 221.62, 717,  131.60, 10,
  3,  -1, 1.208,  237.57, 29, 149.00, 10,
  1,  -1, 2.966,  208.09, 2991, 118.09, 11,
  1,  0,  1.563,  268.31, 91, -7.60,  11,
  2,  -2, .889, 145.16, 1335, 55.17,  11,
  2,  -1, .480, 171.01, 464,  80.95,  11,
  3,  -2, .169, 144.20, 250,  54.00,  11,
  ENDMARK
};

struct kor ARR marskor[] = {
  -1, 1,  .115, 65.84,  684,  156.14, 12,
  -1, 2,  .623, 246.03, 812,  155.77, 12,
  -1, 3,  6.368,  57.60,  556,  -32.06, 12,
  -1, 4,  .588, 57.24,  616,  147.28, 12,
  -2, 5,  .138, 39.18,  157,  309.39, 12,
  -2, 6,  .459, 217.58, 82, 128.10, 12,
  -1, -1, .106, 33.60,  141,  303.45, 13,
  -1, 0,  .873, 34.34,  1112, 304.05, 13,
  -1, 1,  8.559,  35.10,  6947, 304.45, 13,
  -1, 2,  13.966, 20.50,  2875, 113.20, 13,
  -1, 3,  1.487,  22.18,  1619, 112.38, 13,
  -1, 4,  .175, 22.46,  225,  112.15, 13,
  -2, 2,  .150, 18.96,  484,  266.42, 13,
  -2, 3,  7.355,  158.64, 6412, 68.62,  13,
  -2, 4,  4.905,  154.09, 1985, 244.70, 13,
  -2, 5,  .489, 154.39, 543,  244.50, 13,
  -3, 3,  .216, 111.06, 389,  21.10,  13,
  -3, 4,  .355, 110.64, 587,  19.17,  13,
  -3, 5,  2.641,  280.58, 2038, 190.60, 13,
  -3, 6,  .970, 276.06, 587,  6.75, 13,
  -3, 7,  .100, 276.20, 116,  6.40, 13,
  -4, 5,  .152, 232.48, 259,  142.60, 13,
  -4, 6,  .264, 230.47, 387,  139.75, 13,
  -4, 7,  1.156,  41.64,  749,  312.67, 13,
  -4, 8,  .259, 37.92,  205,  128.80, 13,
  -5, 8,  .172, -8.99,  234,  260.70, 13,
  -5, 9,  .575, 164.48, 308,  74.60,  13,
  -6, 10, .115, 113.70, 145,  23.53,  13,
  -6, 11, .363, 285.69, 144,  196.00, 13,
  -7, 13, .353, 48.83,  85, 319.10, 13,
  -8, 15, 1.553,  170.14, 110,  81.00,  13,
  -8, 16, .148, 170.74, 154,  259.94, 13,
  -9, 17, .193, 293.70, 23, 22.80,  13,
  1,  -3, .382, 46.48,  521,  316.25, 14,
  1,  -2, 3.144,  46.78,  3894, 316.39, 14,
  1,  -1, 25.384, 48.96,  23116,  318.87, 14,
  1,  0,  3.732,  -17.62, 1525, 117.81, 14,
  1,  1,  .474, -34.60, 531,  59.67,  14,
  2,  -4, .265, 192.88, 396,  103.12, 14,
  2,  -3, 2.108,  192.72, 3042, 102.89, 14,
  2,  -2, 16.035, 191.90, 22144,  101.99, 14,
  2,  -1, 21.869, 188.35, 16624,  98.33,  14,
  2,  0,  1.461,  189.66, 1478, 279.04, 14,
  2,  1,  .167, 191.04, 224,  280.81, 14,
  3,  -4, .206, 167.11, 338,  76.13,  14,
  3,  -3, 1.309,  168.27, 2141, 76.24,  14,
  3,  -2, 2.607,  228.41, 3437, 139.74, 14,
  3,  -1, 3.174,  207.20, 1915, 115.83, 14,
  3,  0,  .232, 207.78, 240,  298.06, 14,
  4,  -4, .178, 127.25, 322,  36.16,  14,
  4,  -3, .241, 200.69, 389,  110.02, 14,
  4,  -2, .330, 267.57, 413,  179.86, 14,
  4,  -1, .416, 221.88, 184,  128.17, 14,
  1,  -2, .155, -38.20, 191,  231.58, 15,
  1,  -1, 1.351,  -34.10, 1345, 235.85, 15,
  1,  0,  .884, 288.05, 111,  39.90,  15,
  1,  1,  .132, 284.88, 144,  15.67,  15,
  2,  -2, .620, 35.15,  869,  305.30, 15,
  2,  -1, 1.768,  32.50,  1661, 302.51, 15,
  2,  0,  .125, 18.73,  103,  119.90, 15,
  3,  -2, .141, 47.59,  199,  318.06, 15,
  3,  -1, .281, 40.95,  248,  310.75, 15,
  ENDMARK
};

struct m45dat m45[NUM_MOON_CORR] = {
  /*    l,     l',  F,  D,      Long,    Lat,         Par),*/
  { 0,   0, 0,   4,    13.902,     14.06,  0.2607},
  { 0,   0, 0,   2,    2369.912,   2373.36,     28.2333},
  { 1,   0, 0,   4,     1.979,      6.98,  0.0433},
  { 1,   0, 0,   2,     191.953,    192.72,  3.0861},
  { 1,   0, 0,   0,   22639.500,  22609.1,     186.5398},
  { 1,   0, 0,  -2,   -4586.465,  -4578.13,     34.3117},
  { 1,   0, 0,  -4, -38.428,    -38.64,  0.6008},
  { 1,   0, 0,  -6,  -0.393,     -1.43,  0.0086},
  { 0,   1, 0,   4,  -0.289,     -1.59, -0.0053},
  { 0,   1, 0,   2, -24.420,    -25.10, -0.3000},
  { 0,   1, 0,   0,    -668.146,   -126.98,     -0.3997},
  { 0,   1, 0,  -2,    -165.145,   -165.06,      1.9178},
  { 0,   1, 0,  -4,  -1.877,     -6.46,  0.0339},
  { 0,   0, 0,   3,   0.403,     -4.01,  0.0023},
  { 0,   0, 0,   1,    -125.154,   -112.79,     -0.9781},
  { 2,   0, 0,   4,   0.213,      1.02,  0.0054},
  { 2,   0, 0,   2,  14.387,     14.78,  0.2833},
  { 2,   0, 0,   0, 769.016,    767.96, 10.1657},
  { 2,   0, 0,  -2,    -211.656,   -152.53,     -0.3039},
  { 2,   0, 0,  -4, -30.773,    -34.07,  0.3722},
  { 2,   0, 0,  -6,  -0.570,     -1.40,  0.0109},
  { 1,   1, 0,   2,  -2.921,    -11.75, -0.0484},
  { 1,   1, 0,   0,    -109.673,   -115.18,     -0.9490},
  { 1,   1, 0,  -2,    -205.962,   -182.36,      1.4437},
  { 1,  1,  0,  -4,  -4.391,     -9.66,  0.0673},
  { 1,  -1, 0,  4,    0.283,      1.53,  0.0060},
  { 1,  -1, 0,  2,   14.577,     31.70,  0.2302},
  { 1,  -1, 0,  0,  147.687,    138.76,  1.1528},
  { 1,  -1, 0,  -2,  28.475,     23.59, -0.2257},
  { 1,  -1, 0,  -4,   0.636,      2.27, -0.0102},
  { 0,  2,  0,  2,   -0.189,     -1.68, -0.0028},
  { 0,  2,  0,  0,   -7.486,     -0.66, -0.0086},
  { 0,  2,  0,  -2,  -8.096,    -16.35,  0.0918},
  { 0,  0,  2,  2,   -5.741,     -0.04, -0.0009},
  { 0,  0,  2,  0,  -411.608,    -0.2,      -0.0124},
  { 0,  0,  2,  -2, -55.173,    -52.14, -0.1052},
  { 0,  0,  2,  -4,   0.025,     -1.67,  0.0031},
  { 1,  0,  0,  1,   -8.466,    -13.51, -0.1093},
  { 1,  0,  0,  -1,  18.609,      3.59,  0.0118},
  { 1,  0,  0,  -3,   3.215,      5.44, -0.0386},
  { 0,  1,  0,  1,   18.023,     17.93,  0.1494},
  { 0,  1,  0,  -1,   0.560,      0.32, -0.0037},
  { 3,  0,  0,  2,    1.060,      2.96,  0.0243},
  { 3,  0,  0,  0,   36.124,     50.64,  0.6215},
  { 3,  0,  0,  -2, -13.193,    -16.40, -0.1187},
  { 3,  0,  0,  -4,  -1.187,     -0.74,  0.0074},
  { 3,  0,  0,  -6,  -0.293,     -0.31,  0.0046},
  { 2,  1,  0,  2,   -0.290,     -1.45, -0.0051},
  { 2,  1,  0,  0,   -7.649,    -10.56, -0.1038},
  { 2,  1,  0,  -2,  -8.627,     -7.59, -0.0192},
  { 2,  1,  0,  -4,  -2.740,     -2.54,  0.0324},
  { 2,  -1, 0,  2,    1.181,      3.32,  0.0213},
  { 2,  -1, 0,  0,    9.703,     11.67,  0.1268},
  { 2,  -1, 0,  -2,  -2.494,     -1.17, -0.0017},
  { 2,  -1, 0,  -4,   0.360,      0.20, -0.0043},
  { 1,  2,  0,  0,   -1.167,     -1.25, -0.0106},
  { 1,  2,  0,  -2,  -7.412,     -6.12,  0.0484},
  { 1,  2,  0,  -4,  -0.311,     -0.65,  0.0044},
  { 1,  -2, 0,  2,    0.757,      1.82,  0.0112},
  { 1,  -2, 0,  0,    2.580,      2.32,  0.0196},
  { 1,  -2, 0,  -2,   2.533,      2.40, -0.0212},
  { 0,  3,  0,  -2,  -0.344,     -0.57,  0.0036},
  { 1,  0,  2,  2,   -0.992,     -0.02,  0},
  { 1,  0,  2,  0,  -45.099,     -0.02, -0.0010},
  { 1,  0,  2,  -2,  -0.179,     -9.52, -0.0833},
  { 1,  0,  -2, 2,   -6.382,     -3.37, -0.0481},
  { 1,  0,  -2, 0,   39.528,     85.13, -0.7136},
  { 1,  0,  -2, -2,   9.366,      0.71,     -0.0112},
  { 0,  1,  2,  0,    0.415,      0.10,    0.0013},
  { 0,  1,  2,  -2,  -2.152,     -2.26, -0.0066},
  { 0,  1,  -2, 2,   -1.440,     -1.30,  0.0014},
  { 0,  1,  -2, -2,   0.384,      0.0,   0.0},
  { 2,  0,  0,  1,   -0.586,     -1.20, -0.0100},
  { 2,  0,  0,  -1,   1.750,      2.01,  0.0155},
  { 2,  0,  0,  -3,   1.225,      0.91,     -0.0088},
  { 1,  1,  0,  1,    1.267,      1.52,  0.0164},
  { 1,  -1, 0,  -1,  -1.089,      0.55,      0},
  { 0,  0,  2,  -1,   0.584,      8.84,  0.0071},
  { 4,  0,  0,  0,    1.938,      3.60,  0.0401},
  { 4,  0,  0,  -2,  -0.952,     -1.58, -0.0130},
  { 3,  1,  0,  0,   -0.551,      0.94, -0.0097},
  { 3,  1,  0,  -2,  -0.482,     -0.57, -0.0045},
  { 3,  -1, 0,  0,    0.681,      0.96,      0.0115},
  { 2,  0,  2,  0,   -3.996,      0,   0.0004},
  { 2,  0,  2,  -2,   0.557,     -0.75,     -0.0090},
  { 2,  0,  -2, 2,   -0.459,     -0.38, -0.0053},
  { 2,  0,  -2, 0,   -1.298,      0.74,      0.0004},
  { 2,  0,  -2, -2,   0.538,      1.14, -0.0141},
  { 1,  1,  -2, -2,   0.426,      0.07,     -0.0006},
  { 1,  -1, 2,  0,   -0.304,      0.03,  0.0003},
  { 1,  -1, -2, 2,   -0.372,     -0.19, -0.0027},
  { 0,  0,  4,  0,    0.418,      0,   0},
  { 2,  -1, 0,  -1,  -0.352,     -0.37,     -0.0028}
};

# if MOON_TEST_CORR
/* moon additional correction terms */
struct m5dat {    
  REAL8 lng;
  int  i0,i1,i2,i3;
} m5[] = {
  /*    lng,  l,     l',  F,  D,      */
  0.127,  0,  0,  0,  6,
  -0.151, 0,  2,  0,  -4,
  -0.085, 0,  0,  2,  4,
  0.150,  0,  1,  0,  3,
  -0.091, 2,  1,  0,  -6,
  -0.103, 0,  3,  0,  0,
  -0.301, 1,  0,  2,  -4,
  0.202,  1,  0,  -2, -4,
  0.137,  1,  1,  0,  -1,
  0.233,  1,  1,  0,  -3,
  -0.122, 1,  -1, 0,  1,
  -0.276, 1,  -1, 0,  -3,
  0.255,  0,  0,  2,  1,
  0.254,  0,  0,  2,  -3,
  -0.100, 3,  1,  0,  -4,
  -0.183, 3,  -1, 0,  -2,
  -0.297, 2,  2,  0,  -2,
  -0.161, 2,  2,  0,  -4,
  0.197,  2,  -2, 0,  0,
  0.254,  2,  -2, 0,  -2,
  -0.250, 1,  3,  0,  -2,
  -0.123, 2,  0,  2,  2,
  0.173,  2,  0,  -2, -4,
  0.263,  1,  1,  2,  0,
  0.130,  3,  0,  0,  -1,
  0.113,  5,  0,  0,  0,
  0.092,  3,  0,  2,  -2,
  0,  99, 0,  0,  0 /* end mark */
};
# endif /* MOON_TEST_CORR */


/*****************************************************
$Header: deltat.c,v 1.10 93/01/27 14:37:06 alois Exp $
deltat.c
deltat(t): returns delta t (in julian days) from universal time t
is included by users
ET = UT +  deltat

---------------------------------------------------------------
| Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
| The use of this source code is subject to regulations made  |
| by Astrodienst Zurich. The code is NOT in the public domain.|
|                                                             |
| This copyright notice must not be changed or removed        |
| by any user of this program.                                |
---------------------------------------------------------------

******************************************************/

double deltat(jd_ad)
double jd_ad;
{ 
  static short int ARR dt[] = { /* in centiseconds */
  /*
  ** dt from 1637 to 2000, as tabulated in A.E. 
  ** the values 1620 - 1636 are not taken, as they fit
  ** badly the parabola 25.5 t*t for the next range. The
  ** best crossing point to switch to the parabola is
  ** 1637, where we have fitted the value for continuity
  */
  6780, 6500, 6300, 
  6200, 6000, 5800, 5700, 5500,
  5400, 5300, 5100, 5000, 4900,
  4800, 4700, 4600, 4500, 4400,
  4300, 4200, 4100, 4000, 3800, /* 1655 - 59 */
  3700, 3600, 3500, 3400, 3300,
  3200, 3100, 3000, 2800, 2700,
  2600, 2500, 2400, 2300, 2200,
  2100, 2000, 1900, 1800, 1700,
  1600, 1500, 1400, 1400, 1300,
  1200, 1200, 1100, 1100, 1000,
  1000, 1000, 900,  900,  900,
  900,  900,  900,  900,  900,
  900,  900,  900,  900,  900,  /* 1700 - 1704 */
  900,  900,  900,  1000, 1000,
  1000, 1000, 1000, 1000, 1000,
  1000, 1000, 1100, 1100, 1100,
  1100, 1100, 1100, 1100, 1100,
  1100, 1100, 1100, 1100, 1100,
  1100, 1100, 1100, 1100, 1200, /* 1730 - 1734 */
  1200, 1200, 1200, 1200, 1200,
  1200, 1200, 1200, 1200, 1300,
  1300, 1300, 1300, 1300, 1300,
  1300, 1400, 1400, 1400, 1400,
  1400, 1400, 1400, 1500, 1500,
  1500, 1500, 1500, 1500, 1500, /* 1760 - 1764 */
  1600, 1600, 1600, 1600, 1600,
  1600, 1600, 1600, 1600, 1600,
  1700, 1700, 1700, 1700, 1700,
  1700, 1700, 1700, 1700, 1700,
  1700, 1700, 1700, 1700, 1700,
  1700, 1700, 1600, 1600, 1600, /* 1790 - 1794 */
  1600, 1500, 1500, 1400, 1400,
  1370, 1340, 1310, 1290, 1270, /* 1800 - 1804 */
  1260, 1250, 1250, 1250, 1250,
  1250, 1250, 1250, 1250, 1250,
  1250, 1250, 1240, 1230, 1220,
  1200, 1170, 1140, 1110, 1060,
  1020, 960,  910,  860,  800,
  750,  700,  660,  630,  600,  /* 1830 - 1834 */
  580,  570,  560,  560,  560,
  570,  580,  590,  610,  620,
  630,  650,  660,  680,  690,
  710,  720,  730,  740,  750,
  760,  770,  770,  780,  780,
  788,  782,  754,  697,  640,  /* 1860 - 1864 */
  602,  541,  410,  292,  182,
  161,  10, -102, -128, -269,
  -324, -364, -454, -471, -511,
  -540, -542, -520, -546, -546,
  -579, -563, -564, -580, -566,
  -587, -601, -619, -664, -644, /* 1890 - 1894 */
  -647, -609, -576, -466, -374,
  -272, -154, -2, 124,  264,
  386,  537,  614,  775,  913,
  1046, 1153, 1336, 1465, 1601,
  1720, 1824, 1906, 2025, 2095,
  2116, 2225, 2241, 2303, 2349, /* 1920 - 1924 */
  2362, 2386, 2449, 2434, 2408,
  2402, 2400, 2387, 2395, 2386,
  2393, 2373, 2392, 2396, 2402,
  2433, 2483, 2530, 2570, 2624,
  2677, 2728, 2778, 2825, 2871,
  2915, 2957, 2997, 3036, 3072, /* 1950 - 1954 */
  3107, 3135, 3168, 3218, 3268,
  3315, 3359, 3400, 3447, 3503,
  3573, 3654, 3743, 3829, 3920,
  4018, 4117, 4223, 4337, 4449,
  4548, 4646, 4752, 4853, 4959,
  5054, 5138, 5217, 5296, 5379, /* 1980 - 1984 */
  5434, 5487, 5532, 5582, 5630, /* 1985 - 89 from AE 1991 */
  5686, 5757, 5900, 5900, 6000, /* AE 1993 and extrapol */
  6050, 6100, 6150, 6200, 6250, /* 1995 - 1999 */
  6300};          /* 2000 */
  double yr, cy, delta;
  long iyr, i;
  yr = (jd_ad + 18262) / 365.25 + 100.0;    /*  year  relative 1800 */
  cy = yr / 100;
  iyr =  (long) (RFloor(yr) + 1800);   /* truncated to integer, rel 0 */
#if TIDAL_26    /* Stephenson formula only when 26" tidal
  term in lunar motion */
  if (iyr >= 1637  && iyr < 2000) {
  i = iyr - 1637;
  delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
  } else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
  delta = 25.5 * cy * cy  - 25.5 * 4 + 63.00; 
  } else if (iyr >= 948) {  /* from 948 - 1637 use parabola */
  delta = 25.5 * cy * cy; 
  } else {  /* before 984 use other parabola */ 
  delta = 1361.7  + 320 * cy + 44.3 * cy * cy;  /* fits at 948 */
  }
#else    /* use Clemence value + 5 sec before 1690, new dt afterwards */
  cy -= 1;  /* epoch 1900 */
  if (iyr >= 1690  && iyr < 2000) {
  i = iyr - 1637;
  delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
  } else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
  delta = 29.949 * cy * cy  - 29.949 * 4 + 63.0;  
  } else {  
  delta = 5 + 24.349 + 72.3165 * cy + 29.949 * cy * cy; /* fits at 1690 */
  }
#endif
  return delta / 86400.0;
}


/*******************************************
$Header: d2l.c,v 1.9 91/11/16 16:24:20 alois Exp $

double to long with rounding, no overflow check
*************************************/ 
long d2l(x)
double x;
{
  if (x >=0)
    return ((long) (x + 0.5));
  else
    return (-(long) (0.5 - x));
}


/*
* $Header$
*
* A collection of useful functions for centisec calculations.

---------------------------------------------------------------
| Copyright Astrodienst Zurich AG and Alois Treindl, 1991.    |
| The use of this source code is subject to regulations made  |
| by Astrodienst Zurich. The code is NOT in the public domain.|
|                                                             |
| This copyright notice must not be changed or removed        |
| by any user of this program.                                |
---------------------------------------------------------------
*******************************************************/

double degnorm(p)
double p;
{
  if (p < 0) 
    do {
      p += 360.0;
    } while (p < 0);
  else if (p >= 360.0)
    do {
      p -= 360.0;
    } while (p >= 360.0);
  return (p);
}


/*********************************************************
$Header: julday.c,v 1.9 91/11/16 16:25:06 alois Exp $
*********************************************************/

/*
** This function returns the absolute Julian day number (JD)
** for a given calendar date.
** The aruguments are a calendar date: day, month, year as integers,
** hour as double with decimal fraction.
** If gregflag = 1, Gregorian calendar is assumed, gregflag = 0
** Julian calendar is assumed.
**
** The Julian day number is system of numbering all days continously
** within the time range of known human history. It should be familiar
** for every astrological or astronomical programmer. The time variable
** in astronomical theories is usually expressed in Julian days or
** Julian centuries (36525 days per century) relative to some start day;
** the start day is called 'the epoch'.
** The Julian day number is a double representing the number of
** days since JD = 0.0 on 1 Jan -4712, 12:00 noon.
** Midnight has always a JD with fraction .5, because traditionally
** the astronomical day started at noon.
**
** NOTE: The Julian day number is named after the monk Julianus. It must
** not be confused with the Julian calendar system, which is named after
** Julius Cesar, the Roman politician who introduced this calendar.
**
** Original author: Marc Pottenger, Los Angeles.
** with bug fix for year < -4711   15-aug-88 by Alois Treindl
**
** References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
**             Verlag Sterne und Weltraum (1987), p.49 ff
**
** related functions: revjul() reverse Julian day number: compute the
**             calendar date from a given JD
*/

double julday(month, day, year, hour, gregflag)
int month;
int day;
int year;
double hour;
int gregflag;
{
  double jd, u, u0, u1, u2;

  u = year;
  if (month < 3)
    u -=1;
  u0 = u + 4712.0;
  u1 = month + 1.0;
  if (u1 < 4)
    u1 += 12.0;
  jd = RFloor(u0*365.25)
    + RFloor(30.6*u1+0.000001)
    + day + hour/24.0 - 63.5;
  if (gregflag) {
    u2 = RFloor(ABS8(u) / 100) - RFloor(ABS8(u) / 400);
    if (u < 0.0)
      u2 = -u2;
    jd = jd - u2 + 2;            
    if ((u < 0.0) && (u/100 == RFloor(u/100)) && (u/400 != RFloor(u/400)))
      jd -= 1;
  }
  return jd;
}


/*********************************************************
$Header: revjul.c,v 1.9 91/11/16 16:25:37 alois Exp $
*********************************************************/

/*
** revjul() is the inverse function to julday(), see the description there.
** Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
** return values are the calendar day, month, year and the hour of
** the day with decimal fraction (0 .. 23.999999).
**
** Original author Mark Pottenger, Los Angeles.
** with bug fix for year < -4711 16-aug-88 Alois Treindl
*/

void revjul(jd, gregflag, jmon, jday, jyear, jut)
double jd;
int gregflag;
int *jmon;
int *jday;
int *jyear;
double *jut;
{
  double u0, u1, u2, u3, u4;

  u0 = jd + 32082.5;
  if (gregflag) {
    u1 = u0 + RFloor(u0/36525.0) - RFloor(u0/146100.0) - 38.0;
    if (jd >= 1830691.5) u1 +=1;
      u0 = u0 + RFloor(u1/36525.0) - RFloor(u1/146100.0) - 38.0;
  }
  u2 = RFloor(u0 + 123.0);
  u3 = RFloor((u2 - 122.2) / 365.25);
  u4 = RFloor((u2 - RFloor(365.25 * u3)) / 30.6001);
  *jmon = (int)(u4-1.0);
  if (*jmon > 12)
    *jmon -= 12;
  *jday = (int)(u2 - RFloor(365.25 * u3) - RFloor(30.6001 * u4));
  *jyear = (int)(u3 + RFloor((u4 - 2.0) / 12.0) - 4800);
  *jut = (jd - RFloor(jd + 0.5) + 0.5) * 24.0;
}
#endif /* PLACALC */

/* placalc2.c */