home *** CD-ROM | disk | FTP | other *** search
/ Nebula 2 / Nebula Two.iso / Apps / Astro / astrolog / Source / placalc2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-12  |  35.6 KB  |  1,059 lines

  1. /*
  2. ** Astrolog (Version 4.10) File: placalc2.c
  3. **
  4. ** IMPORTANT NOTICE: the graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1994 by Walter D. Pullen
  6. ** (cruiser1@stein.u.washington.edu). Permission is granted to freely
  7. ** use and distribute these routines provided one doesn't sell,
  8. ** restrict, or profit from them in any way. Modification is allowed
  9. ** provided these notices remain with any altered or edited versions of
  10. ** the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 3/19/1994.
  36. */
  37.  
  38. #ifdef PLACALC
  39. #ifdef ASTROLOG
  40. /* Begin contents of placalc2.c */
  41. #endif
  42. /*****************************************************
  43. $Header: placalc.c,v 1.14 93/07/19 22:13:07 alois Exp $
  44.  
  45.  
  46.   ---------------------------------------------------------------
  47.   | Copyright Astrodienst AG and Alois Treindl, 1989,1991,1993  |
  48.   | The use of this source code is subject to regulations made    |
  49.   | by Astrodienst Zurich. The code is NOT in the public domain.|
  50.   |                                |
  51.   | This copyright notice must not be changed or removed    |
  52.   | by any user of this program.                |
  53.   ---------------------------------------------------------------
  54.  
  55. Important changes:
  56. 11-jun-93 revision 1.12: fixed error which affected Mercury between -2100 and
  57.             -3100 (it jumped wildly).
  58.  
  59. *******************************************************/
  60.  
  61. #ifndef ASTROLOG
  62. #include "placalc.h"    /* includes ourdef.h */
  63. #include "astrolib.h"    /* includes ourdef.h */
  64. #include "helconst.c"    /* orbital elements and disturbations
  65.             for inner planets and moon */
  66. #include "deltat.c"
  67. #endif
  68. #ifdef ASTROLOG
  69. /* Begin contents of helconst.c */
  70. #endif
  71. /***********************************************************
  72.  * $Header$
  73.  
  74.  definition module for planetary elements
  75.  and disturbation coefficients
  76.  version HP-UX C  for new version with stored outer planets
  77.  31-jul-88
  78.  by Alois Treindl 
  79.  
  80.   ---------------------------------------------------------------
  81.   | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
  82.   | The use of this source code is subject to regulations made    |
  83.   | by Astrodienst Zurich. The code is NOT in the public domain.|
  84.   |                                |
  85.   | This copyright notice must not be changed or removed    |
  86.   | by any user of this program.                |
  87.   ---------------------------------------------------------------
  88.  
  89. ***********************************************************/
  90. /* In the elements degrees were kept as the units for the constants. This
  91. requires conversion to radians, when the actual calculations are performed.
  92. This approach is not the most efficient, but safer for development.
  93. Constant conversion could be done by writing all degree constants with
  94.  value * DEGTORAD */
  95.  
  96. # define TIDAL_26    TRUE    /* decide wheter to use new or old lunar tidal
  97.                 term; a consistent system of delta t must be 
  98.                 used */
  99. # define MOON_TEST_CORR FALSE    /* to include more lunar terms in longitude */
  100.  
  101. REAL8  ekld [4] = { 23.452294, -46.845, -.0059, 0.00181 };
  102.     /* ecliptic with epoch1900, Ekd(0..3) in basic */
  103. struct eledata {
  104.   REAL8 axis,        /* mean distance, in a.u., A(N) in basic */
  105.     period,        /* days for one revolution, P(N) in basic */
  106.     epoch,        /* relative juldate of epoch, Ep(N) in basic */
  107.             /* T = distance to epoch in jul.centuries 36525 day*/
  108.     lg0,lg1,lg2,lg3,/* deg(epoch), degree/day, seconds/T^2,seconds/T^3 */
  109.             /* Pd(N,0..2) in basic, lg3 was not present */
  110.     pe0,pe1,pe2,pe3,/* deg(epoch), seconds/T,  seconds/T^2,seconds/T^3 */
  111.             /* Pd(N,3..5) in basic, pe3 was not present */
  112.     ex0,ex1,ex2,    /* ecl(epoch), 1/T, 1/T^2    */
  113.             /* Pd(N,6..8) in basic */
  114.     kn0,kn1,kn2,kn3,/* node(epoch),seconds/T,  seconds/T^2,seconds/T^3 */
  115.             /* Pd(N,9..11) in basic, kn3 was not present */
  116.     in0,in1,in2;    /* incl(epoch),1/T, 1/T^2    */
  117.             /* Pd(N,12..14) in basic */
  118.     } pd [MARS + 1] = {
  119.     {/*earth*/    1.00000023,    365.25636042,    EPOCH1900,
  120.       99.696678,    .9856473354,    1.089,        0,    
  121.       101.220833,    6189.03,    1.63,        0.012,
  122.       0.01675104,    -0.00004180,    -0.000000126,
  123.       0,        0,        0,        0,
  124.       0,        0,        0},
  125.       /* note 29 June 88 by Alois: G.M.Clemence, Astronomical Journal
  126.       vol.53,p. 178 (1948) gives a correction to the perihel motion
  127.       of -4.78" T, giving 6184.25 for the linear Term above. We have
  128.       not yet applied this correction. It has been used in APAE 22,4
  129.       on the motion of mars and does make an official impression. */
  130.     {/*moon*/    0.0025955307,    27.321661,    EPOCH1900,
  131. # if ! TIDAL_26 
  132.     /* values from Improved Lunar Ephemeris, corresponding to tidal
  133.        term -22.44"/cy and  consistent with delta t ~ 29.949 T*T
  134.     */
  135.       270.4341638,    13.176396526808121, -4.08,    0.0068,
  136. # endif
  137. # if TIDAL_26 
  138.     /* new values from Morrison 1979, with tidal term -26"/cy as
  139.        stated in A.E. 1986 onwards, consistent with delta t ~ 44.3 T*T 
  140.        correction: -1.54" + 2.33" T - 1.78" T*T
  141.     */
  142.       270.4337361,    13.176396544528099, -5.86,    0.0068,
  143. # endif
  144.       334.329556,    14648522.52,    -37.17,        -0.045,
  145.       0.054900489,    0,        0,
  146.       259.183275,    -6962911.23,    7.48 ,        0.008,
  147.       5.145388889,    0,        0},
  148.     {/*mercury*/    .3870986,    87.969252,    EPOCH1900,
  149.       178.179078,    4.0923770233,    1.084,        0,
  150.       75.89969722,    5599.76,    1.061,        0,
  151.       0.20561421,    .00002046,     -.000000030,
  152.       47.145944444,    4266.75,    .626,        0,
  153.       7.0028805555,    6.699,        -.066},
  154.     {/*venus*/    .72333162,    224.700726,    EPOCH1900,
  155.       342.767053,    1.6021687039,    1.1148,        0,
  156.       130.16383333,    5068.93,    -3.515,        0,
  157.       0.00682069,    -.00004774,    .000000091,
  158.       75.7796472223,3239.46,    1.476,        0,
  159.       3.3936305555,    3.621,        .0035},
  160.     {/*mars*/    1.5236914620,    686.9296097,    EPOCH1900,
  161.       /* These are the corrected elements by Ross */
  162.       293.74762778,    .524071163814,    1.1184,        0,
  163.       334.21820278,    6626.73,    .4675,        -0.0043,
  164.       0.09331290,    .000092064,    -.000000077,
  165.       48.786441667,    2775.57,    -.005,        -0.0192,
  166.       1.85033333,    -2.430,        .0454}
  167. };
  168.  
  169. /*
  170.  * mimimum and maximum distances computed over 1000 years with plamimax,
  171.  * required for relative distances rgeo, where the distance is given
  172.  * as 100 when a planet is closest and as 0 when farthest from earth.
  173.  */
  174. REAL8 rmima[CALC_N][2] = {    
  175.     { 0.98296342,  1.01704665},
  176.     { 0.00238267,  0.00271861},
  177.     { 0.54900496,  1.45169607},
  178.     { 0.26411287,  1.73597885},
  179.     { 0.37289847,  2.67626927},
  180.     { 3.94877993,  6.45627627},
  181.     { 7.99362824, 11.09276636},
  182.     {17.28622633, 21.10714104},
  183.     {28.81374786, 31.33507284},
  184.     {28.67716748, 50.29208774},
  185.     { 0.00,     0.00259553},    /* nodes don't get a real value*/
  186.     { 0.00,  0.00259553},
  187.     { 7.36277475, 19.86585062}};
  188.       
  189.     
  190.  
  191. #define SDNUM 20
  192. struct sdat {        /* 0..19 mean anomalies of disturbing planets 
  193.             Sd(0..19,0..1) in basic */
  194.     REAL8 sd0,    /* mean anomaly at epoch 1850 */
  195.           sd1;    /* degrees/year */
  196.     } _sd [SDNUM] = {
  197.      114.50,    585.17493,
  198.      109.856,    191.39977,
  199.      148.031,    30.34583,
  200.      284.716,    12.21794,
  201.      114.508,    585.17656,
  202.      -0.56,        359.99213,
  203.      148.03,    30.34743,
  204.      284.72,    12.2196,
  205.      248.07,    1494.726615,
  206.      359.44,    359.993595,
  207.      109.86,    191.402867,
  208.      148.02,    30.348930,
  209.      114.503,    585.173715,
  210.      359.444,    359.989285,
  211.      148.021,    30.344620,
  212.      284.716,    12.21669,
  213.      148.0315,    30.34906264,
  214.      284.7158,    12.22117085,
  215.      220.1695,     4.284931111,
  216.      291.8024,     2.184704167
  217. };
  218.  
  219. REAL8  sa [SDNUM];
  220.  
  221. struct kor {
  222.     int     j, i;
  223.     REAL8   lampl;    /* amplitude of disturbation in long, seconds of arc */
  224.     REAL8    lphase; /* phase of disturbation in long, degrees */
  225.     INT4    rampl;  /* ampl. of disturbation in radius, 9th place of log */
  226.     REAL8   rphase; /* phase of disturbation in radius, degrees */
  227.     int     k;    /* index into disturbing planet anomaly table sa[] */
  228. };
  229.        /* delta long = lampl * COS (lphase - arg) in seconds of arc
  230.       delta rad  = rampl * COS (rphase - arg) in ninth place of log
  231.       arg = j * sa (k) + i * ma (this planet)
  232.       ma = mean anomaly
  233.       sa = mean anomaly of disturbing planet, where this
  234.            is taken from the aproximate value in sa[]
  235.       For the COS (phase - arg) it is good enough to compute
  236.       with 32 bit reals, because ampl and phase have only 
  237.       four to five significant digits. 
  238.       While saving constant space, it is costing execution time due
  239.       to float/double conversions.
  240.     */
  241.         /* In basic, all correction terms for sun, mercury, venus and mars
  242.        were contained in one array K(0..142,0..6); Nk(N,0) contained
  243.        the index of the first term of planet N and Nk(N,1) the number
  244.        of terms for this planet. Here, we use a  0 in the first column
  245.        kor.j to indicate the end of the table for a planet.
  246.        K(*) was a basic INTEGER array, therefore the amplitudes and phases
  247.        had to be expressed as
  248.        K(i,2) = ampl. of longitude in 0.001 seconds of arc
  249.        K(i,3) = phase of longitude in 0.01 degrees
  250.        K(i,4) = ampl. of radius in 9th place of log
  251.        K(i,5) = phase of radius in 0.01 degrees.
  252.        Here we have converted the amplitude of long. to seconds of arc
  253.        and the phases to degrees.
  254.      */
  255.  
  256. struct kor earthkor[] = {    /* 11-jul-88 all terms to 0.020" longitude */
  257.       /*  j     i       lampl   lphase  rampl   rphase  k */
  258.      -1,    1,    0.013,    243,    28,    335,    8,    /* mercury */
  259.      -1,    3,    0.015,    357,    18,    267,    8,
  260.      -1,    4,    0.023,    326,    5,    239,    8,
  261.      -1,    0,    0.075,    296.6,    94,    205.0,    0,    /* Venus */
  262.      -1,    1,    4.838,    299.10,    2359,    209.08,    0,    
  263.      -1,    2,    0.074,    207.9,    69,    348.5,    0,
  264.      -1,    3,    0.009,    249,    16,    330,    0,
  265.         -2,    1,    .116,    148.90,    160,    58.40,    0,    
  266.      -2,    2,    5.526,    148.31,    6842,    58.32,    0,    
  267.      -2,    3,    2.497,    315.94,    869,    226.70,    0,    
  268.      -2,    4,    0.044,    311.4,    52,    38.8,    0,
  269.      -3,    2,    0.013,    176,    21,    90,    0,
  270.      -3,    3,    .666,    177.71,    1045,    87.57,    0,    
  271.      -3,    4,    1.559,    -14.75,    1497,    255.25,    0,    
  272.      -3,    5,    1.024,    318.15,    194,    49.50,    0,    
  273.      -3,    6,    0.017,    315,    19,    43,    0,
  274.      -4,    4,    .210,    206.20,    376,    116.28,    0,    
  275.      -4,    5,    .144,    195.40,    196,    105.20,    0,    
  276.      -4,    6,    .152,    -16.20,    94,    254.80,    0,    
  277.      -5,    5,    0.084,    235.6,    163,    145.4,    0,
  278.      -5,    6,    0.037,    221.8,    59,    132.2,    0,
  279.      -5,    7,    .123,    195.30,    141,    105.40,    0,    
  280.      -5,    8,    .154,    -.40,    26,    270.00,    0,    
  281.      -6,    6,    0.038,    264.1,    80,    174.3,    0,
  282.      -6,    7,    0.014,    253,    25,    164,    0,
  283.      -6,    8,    0.01,    230,    14,    135,    0,
  284.      -6,    9,    0.014,    12,    12,    284,    0,
  285.      -7,    7,    0.020,    294,    42,    203.5,    0,
  286.      -7,    8,    0.006,    279,    12,    194,    0,
  287.      -8,    8,    0.011,    322,    24,    234,    0,
  288.      -8,    12,    0.042,    259.2,    44,    169.7,    0,
  289.      -8,    14,    0.032,    48.8,    33,    138.7,    0,
  290.      -9,    9,    0.006,    351,    13,    261,    0,
  291.      1,    -1,    .273,    217.70,    150,    127.70,    1,    /* mars */
  292.      1,    0,    0.048,    260.3,    28,    347,    1,
  293.      2,    -3,    0.041,    346,    52,    255.4,    1,
  294.      2,    -2,    2.043,    343.89,    2057,    253.83,    1,
  295.      2,    -1,    1.770,    200.40,    151,    295.00,    1,
  296.      2,    0,    0.028,    148,    31,    234.3,    1,  
  297.      3,    -3,    .129,    294.20,    168,    203.50,    1,
  298.      3,    -2,    .425,    -21.12,    215,    249.00,    1,
  299.      4,    -4,    0.034,    71,    49,    339.7,    1,
  300.      4,    -3,    .500,    105.18,    478,    15.17,    1,
  301.      4,    -2,    .585,    -25.94,    105,    65.90,    1,
  302.      5,    -4,    0.085,    54.6,    107,    324.6,    1,
  303.      5,    -3,    .204,    100.80,    89,    11.00,    1,
  304.      6,    -5,    0.020,    186,    30,    95.7,    1,
  305.      6,    -4,    .154,    227.40,    139,    137.30,    1,
  306.      6,    -3,    .101,    96.30,    27,    188.00,    1,
  307.      7,    -5,    0.049,    176.5,    60,    86.2,    1,
  308.      7,    -4,    .106,    222.70,    38,    132.90,    1,
  309.      8,    -5,    0.052,    348.9,    45,    259.7,    1,
  310.      8,    -4,    0.021,    215.2,    8,    310,    1,
  311.      8,    -6,    0.010,    307,    15,    217,    1,
  312.      9,    -6,    0.028,    298,    34,    208.1,    1,
  313.      9,    -5,    0.062,    346,    17,    257,    1,
  314.      10,    -6,    0.019,    111,    15,    23,    1,
  315.      11,    -7,    0.017,    59,    20,    330,    1,
  316.      11,    -6,    0.044,    105.9,    9,    21,    1,
  317.      13,    -8,    0.013,    184,    15,    94,    1,
  318.      13,    -7,    0.045,    227.8,    5,    143,    1,
  319.      15,    -9,    0.021,    309,    22,    220,    1,
  320.      17,    -9,    0.026,    113,    0,    0,    1,
  321.      1,    -2,    .163,    198.60,    208,    112.00,    2,    /* jupiter */
  322.      1,    -1,    7.208,    179.53,    7067,    89.55,    2,
  323.      1,    0,    2.600,    263.22,    244,    -21.40,    2,
  324.      1,    1,    0.073,    276.3,    80,    6.5,    2,
  325.      2,    -3,    0.069,    80.8,    103,    350.5,    2,
  326.      2,    -2,    2.731,    87.15,    4026,    -2.89,    2,
  327.      2,    -1,    1.610,    109.49,    1459,    19.47,    2,
  328.      2,    0,    0.073,    252.6,    8,    263,    2,
  329.      3,    -3,    .164,    170.50,    281,    81.20,    2,
  330.      3,    -2,    .556,    82.65,    803,    -7.44,    2,
  331.      3,    -1,    .210,    98.50,    174,    8.60,    2,
  332.      4,    -4,    0.016,    259,    29,    170,    2,
  333.      4,    -3,    0.044,    168.2,    74,    79.9,    2,
  334.      4,    -2,    0.080,    77.7,    113,    347.7,    2,
  335.      4,    -1,    0.023,    93,    17,    3,    2,
  336.      5,    -2,    0.009,    71,    14,    343,    2,
  337.      1,    -2,    0.011,    105,    15,    11,    3,    /* saturn */
  338.      1,    -1,    .419,    100.58,    429,    10.60,    3,
  339.      1,    0,    .320,    269.46,    8,    -7.00,    3,
  340.      2,    -2,    .108,    290.60,    162,    200.60,    3,
  341.      2,    -1,    .112,    293.60,    112,    203.10,    3,
  342.      3,    -2,    0.021,    289,    32,    200.1,    3,
  343.      3,    -1,    0.017,    291,    17,    201,    3,
  344.      ENDMARK
  345. };
  346.  
  347. struct kor mercurykor[] = {
  348.      1,    -1,    .711,    35.47,    491,    305.28,    4,
  349.      2,    -3,    .552,    161.15,    712,    71.12,    4,
  350.      2,    -2,    2.100,    161.15,    2370,    71.19,    4,
  351.      2,    -1,    3.724,    160.69,    899,    70.49,    4,
  352.      2,    0,    .729,    159.76,    763,    250.00,    4,
  353.      3,    -3,    .431,    105.37,    541,    15.53,    4,
  354.      3,    -2,    1.329,    104.78,    1157,    14.84,    4,
  355.      3,    -1,    .539,    278.95,    14,    282.00,    4,
  356.      4,    -2,    .484,    226.40,    234,    136.02,    4,
  357.      5,    -4,    .685,    -10.43,    849,    259.51,    4,
  358.      5,    -3,    2.810,    -10.14,    2954,    259.92,    4,
  359.      5,    -2,    7.356,    -12.22,    282,    255.43,    4,
  360.      5,    -1,    1.471,    -12.30,    1550,    77.75,    4,
  361.      5,    0,    .375,    -12.29,    472,    77.70,    4,
  362.      2,    -1,    .443,    218.48,    256,    128.36,    5,
  363.      4,    -2,    .374,    151.81,    397,    61.63,    5,
  364.      4,    -1,    .808,    145.93,    13,    35.00,    5,
  365.      1,    -1,    .697,    181.07,    708,    91.38,    6,
  366.      1,    0,    .574,    236.72,    75,    265.40,    6,
  367.      2,    -2,    .938,    36.98,    1185,    306.97,    6,
  368.      2,    -1,    3.275,    37.00,    3268,    306.99,    6,
  369.      2,    0,    .499,    31.91,    371,    126.90,    6,
  370.      3,    -1,    .353,    25.84,    347,    295.76,    6,
  371.      2,    -1,    .380,    239.87,    0,    0,    7,
  372.      ENDMARK
  373. };
  374.  
  375. struct kor venuskor[] = {
  376.      -1,    2,    .264,    -19.20,    175,    251.10,    8,
  377.      -2,    5,    .361,    167.68,    55,    77.20,    8,
  378.      1,    -1,    4.889,    119.11,    2246,    29.11,    9,
  379.      2,    -2,    11.261,    148.23,    9772,    58.21,    9,
  380.      3,    -3,    7.128,    -2.57,    8271,    267.42,    9,
  381.      3,    -2,    3.446,    135.91,    737,    47.37,    9,
  382.      4,    -4,    1.034,    26.54,    1426,    296.49,    9,
  383.      4,    -3,    .677,    165.32,    445,    75.70,    9,
  384.      5,    -5,    .330,    56.88,    510,    -33.36,    9,
  385.      5,    -4,    1.575,    193.93,    1572,    104.21,    9,
  386.      5,    -3,    1.439,    138.08,    162,    229.90,    9,
  387.      6,    -6,    .143,    84.40,    236,    -5.80,    9,
  388.      6,    -5,    .205,    44.20,    256,    314.20,    9,
  389.      6,    -4,    .176,    164.30,    70,    75.70,    9,
  390.      8,    -5,    .231,    180.00,    25,    75.00,    9,
  391.      3,    -2,    .673,    221.62,    717,    131.60,    10,
  392.      3,    -1,    1.208,    237.57,    29,    149.00,    10,
  393.      1,    -1,    2.966,    208.09,    2991,    118.09,    11,
  394.      1,    0,    1.563,    268.31,    91,    -7.60,    11,
  395.      2,    -2,    .889,    145.16,    1335,    55.17,    11,
  396.      2,    -1,    .480,    171.01,    464,    80.95,    11,
  397.      3,    -2,    .169,    144.20,    250,    54.00,    11,
  398.      ENDMARK
  399. };
  400.  
  401. struct kor marskor[] = {
  402.      -1,    1,    .115,    65.84,    684,    156.14,    12,
  403.      -1,    2,    .623,    246.03,    812,    155.77,    12,
  404.      -1,    3,    6.368,    57.60,    556,    -32.06,    12,
  405.      -1,    4,    .588,    57.24,    616,    147.28,    12,
  406.      -2,    5,    .138,    39.18,    157,    309.39,    12,
  407.      -2,    6,    .459,    217.58,    82,    128.10,    12,
  408.      -1,    -1,    .106,    33.60,    141,    303.45,    13,
  409.      -1,    0,    .873,    34.34,    1112,    304.05,    13,
  410.      -1,    1,    8.559,    35.10,    6947,    304.45,    13,
  411.      -1,    2,    13.966,    20.50,    2875,    113.20,    13,
  412.      -1,    3,    1.487,    22.18,    1619,    112.38,    13,
  413.      -1,    4,    .175,    22.46,    225,    112.15,    13,
  414.      -2,    2,    .150,    18.96,    484,    266.42,    13,
  415.      -2,    3,    7.355,    158.64,    6412,    68.62,    13,
  416.      -2,    4,    4.905,    154.09,    1985,    244.70,    13,
  417.      -2,    5,    .489,    154.39,    543,    244.50,    13,
  418.      -3,    3,    .216,    111.06,    389,    21.10,    13,
  419.      -3,    4,    .355,    110.64,    587,    19.17,    13,
  420.      -3,    5,    2.641,    280.58,    2038,    190.60,    13,
  421.      -3,    6,    .970,    276.06,    587,    6.75,    13,
  422.      -3,    7,    .100,    276.20,    116,    6.40,    13,
  423.      -4,    5,    .152,    232.48,    259,    142.60,    13,
  424.      -4,    6,    .264,    230.47,    387,    139.75,    13,
  425.      -4,    7,    1.156,    41.64,    749,    312.67,    13,
  426.      -4,    8,    .259,    37.92,    205,    128.80,    13,
  427.      -5,    8,    .172,    -8.99,    234,    260.70,    13,
  428.      -5,    9,    .575,    164.48,    308,    74.60,    13,
  429.      -6,    10,    .115,    113.70,    145,    23.53,    13,
  430.      -6,    11,    .363,    285.69,    144,    196.00,    13,
  431.      -7,    13,    .353,    48.83,    85,    319.10,    13,
  432.      -8,    15,    1.553,    170.14,    110,    81.00,    13,
  433.      -8,    16,    .148,    170.74,    154,    259.94,    13,
  434.      -9,    17,    .193,    293.70,    23,    22.80,    13,
  435.      1,    -3,    .382,    46.48,    521,    316.25,    14,
  436.      1,    -2,    3.144,    46.78,    3894,    316.39,    14,
  437.      1,    -1,    25.384,    48.96,    23116,    318.87,    14,
  438.      1,    0,    3.732,    -17.62,    1525,    117.81,    14,
  439.      1,    1,    .474,    -34.60,    531,    59.67,    14,
  440.      2,    -4,    .265,    192.88,    396,    103.12,    14,
  441.      2,    -3,    2.108,    192.72,    3042,    102.89,    14,
  442.      2,    -2,    16.035,    191.90,    22144,    101.99,    14,
  443.      2,    -1,    21.869,    188.35,    16624,    98.33,    14,
  444.      2,    0,    1.461,    189.66,    1478,    279.04,    14,
  445.      2,    1,    .167,    191.04,    224,    280.81,    14,
  446.      3,    -4,    .206,    167.11,    338,    76.13,    14,
  447.      3,    -3,    1.309,    168.27,    2141,    76.24,    14,
  448.      3,    -2,    2.607,    228.41,    3437,    139.74,    14,
  449.      3,    -1,    3.174,    207.20,    1915,    115.83,    14,
  450.      3,    0,    .232,    207.78,    240,    298.06,    14,
  451.      4,    -4,    .178,    127.25,    322,    36.16,    14,
  452.      4,    -3,    .241,    200.69,    389,    110.02,    14,
  453.      4,    -2,    .330,    267.57,    413,    179.86,    14,
  454.      4,    -1,    .416,    221.88,    184,    128.17,    14,
  455.      1,    -2,    .155,    -38.20,    191,    231.58,    15,
  456.      1,    -1,    1.351,    -34.10,    1345,    235.85,    15,
  457.      1,    0,    .884,    288.05,    111,    39.90,    15,
  458.      1,    1,    .132,    284.88,    144,    15.67,    15,
  459.      2,    -2,    .620,    35.15,    869,    305.30,    15,
  460.      2,    -1,    1.768,    32.50,    1661,    302.51,    15,
  461.      2,    0,    .125,    18.73,    103,    119.90,    15,
  462.      3,    -2,    .141,    47.59,    199,    318.06,    15,
  463.      3,    -1,    .281,    40.95,    248,    310.75,    15,
  464.      ENDMARK
  465. };
  466.  
  467. #define NUM_MOON_CORR 93
  468. /* moon correction data; revised 30-jul-88: all long. to 0.3" */
  469. struct m45dat {        
  470.     int  i0,i1,i2,i3;
  471.     REAL8 lng,lat,par;
  472.     } m45 [NUM_MOON_CORR] = {  
  473.     /*      l,     l',    F,    D,      Long,       Lat,            Par),*/
  474.     { 0,     0,    0,     4,      13.902,     14.06,     0.2607},
  475.     { 0,     0,    0,     2,    2369.912,   2373.36,     28.2333},
  476.     { 1,     0,    0,     4,       1.979,      6.98,     0.0433},
  477.     { 1,     0,    0,     2,     191.953,    192.72,     3.0861},
  478.     { 1,     0,    0,     0,   22639.500,  22609.1,     186.5398},
  479.     { 1,     0,    0,    -2,   -4586.465,  -4578.13,     34.3117},
  480.     { 1,     0,    0,    -4,    -38.428,    -38.64,     0.6008},
  481.     { 1,     0,    0,    -6,     -0.393,     -1.43,     0.0086},
  482.     { 0,     1,    0,     4,     -0.289,     -1.59,    -0.0053},
  483.     { 0,     1,    0,     2,    -24.420,    -25.10,    -0.3000},
  484.     { 0,     1,    0,     0,    -668.146,   -126.98,     -0.3997},
  485.     { 0,     1,    0,    -2,    -165.145,   -165.06,      1.9178},
  486.     { 0,     1,    0,    -4,     -1.877,     -6.46,     0.0339},
  487.     { 0,     0,    0,     3,      0.403,     -4.01,     0.0023},
  488.     { 0,     0,    0,     1,    -125.154,   -112.79,     -0.9781},
  489.     { 2,     0,    0,     4,      0.213,      1.02,     0.0054},
  490.     { 2,     0,    0,     2,     14.387,     14.78,     0.2833},
  491.     { 2,     0,    0,     0,    769.016,    767.96,    10.1657},
  492.     { 2,     0,    0,    -2,    -211.656,   -152.53,     -0.3039},
  493.     { 2,     0,    0,    -4,    -30.773,    -34.07,     0.3722},
  494.     { 2,     0,    0,    -6,     -0.570,     -1.40,     0.0109},
  495.     { 1,     1,    0,     2,     -2.921,    -11.75,    -0.0484},
  496.     { 1,     1,    0,     0,    -109.673,   -115.18,     -0.9490},
  497.     { 1,     1,    0,    -2,    -205.962,   -182.36,      1.4437},
  498.     { 1,    1,    0,    -4,     -4.391,     -9.66,     0.0673},
  499.     { 1,    -1,    0,    4,      0.283,      1.53,     0.0060},
  500.     { 1,    -1,    0,    2,     14.577,     31.70,     0.2302},
  501.     { 1,    -1,    0,    0,    147.687,    138.76,     1.1528},
  502.     { 1,    -1,    0,    -2,     28.475,     23.59,    -0.2257},
  503.     { 1,    -1,    0,    -4,      0.636,      2.27,    -0.0102},
  504.     { 0,    2,    0,    2,     -0.189,     -1.68,    -0.0028},
  505.     { 0,    2,    0,    0,     -7.486,     -0.66,    -0.0086},
  506.     { 0,    2,    0,    -2,     -8.096,    -16.35,     0.0918},
  507.     { 0,    0,    2,    2,     -5.741,     -0.04,    -0.0009},
  508.     { 0,    0,    2,    0,    -411.608,    -0.2,      -0.0124},
  509.     { 0,    0,    2,    -2,    -55.173,    -52.14,    -0.1052},
  510.     { 0,    0,    2,    -4,      0.025,     -1.67,     0.0031},
  511.     { 1,    0,    0,    1,     -8.466,    -13.51,    -0.1093},
  512.     { 1,    0,    0,    -1,     18.609,      3.59,     0.0118},
  513.     { 1,    0,    0,    -3,      3.215,      5.44,    -0.0386},
  514.     { 0,    1,    0,    1,     18.023,     17.93,     0.1494},
  515.     { 0,    1,    0,    -1,      0.560,      0.32,    -0.0037},
  516.     { 3,    0,    0,    2,      1.060,      2.96,     0.0243},
  517.     { 3,    0,    0,    0,     36.124,     50.64,     0.6215},
  518.     { 3,    0,    0,    -2,    -13.193,    -16.40,    -0.1187},
  519.     { 3,    0,    0,    -4,     -1.187,     -0.74,     0.0074},
  520.     { 3,    0,    0,    -6,     -0.293,     -0.31,     0.0046},
  521.     { 2,    1,    0,    2,     -0.290,     -1.45,    -0.0051},
  522.     { 2,    1,    0,    0,     -7.649,    -10.56,    -0.1038},
  523.     { 2,    1,    0,    -2,     -8.627,     -7.59,    -0.0192},
  524.     { 2,    1,    0,    -4,     -2.740,     -2.54,     0.0324},
  525.     { 2,    -1,    0,    2,      1.181,      3.32,     0.0213},
  526.     { 2,    -1,    0,    0,      9.703,     11.67,     0.1268},
  527.     { 2,    -1,    0,    -2,     -2.494,     -1.17,    -0.0017},
  528.     { 2,    -1,    0,    -4,      0.360,      0.20,    -0.0043},
  529.     { 1,    2,    0,    0,     -1.167,     -1.25,    -0.0106},
  530.     { 1,    2,    0,    -2,     -7.412,     -6.12,     0.0484},
  531.     { 1,    2,    0,    -4,     -0.311,     -0.65,     0.0044},
  532.     { 1,    -2,    0,    2,      0.757,      1.82,     0.0112},
  533.     { 1,    -2,    0,    0,      2.580,      2.32,     0.0196},
  534.     { 1,    -2,    0,    -2,      2.533,      2.40,    -0.0212},
  535.     { 0,    3,    0,    -2,     -0.344,     -0.57,     0.0036},
  536.     { 1,    0,    2,    2,     -0.992,     -0.02,     0},
  537.     { 1,    0,    2,    0,    -45.099,     -0.02,    -0.0010},
  538.     { 1,    0,    2,    -2,     -0.179,     -9.52,    -0.0833},
  539.     { 1,    0,    -2,    2,     -6.382,     -3.37,    -0.0481},
  540.     { 1,    0,    -2,    0,     39.528,     85.13,    -0.7136},
  541.     { 1,    0,    -2,    -2,      9.366,      0.71,     -0.0112},
  542.     { 0,    1,    2,    0,      0.415,      0.10,      0.0013},
  543.     { 0,    1,    2,    -2,     -2.152,     -2.26,    -0.0066},
  544.     { 0,    1,    -2,    2,     -1.440,     -1.30,     0.0014},
  545.     { 0,    1,    -2,    -2,      0.384,      0.0,     0.0},
  546.     { 2,    0,    0,    1,     -0.586,     -1.20,    -0.0100},
  547.     { 2,    0,    0,    -1,      1.750,      2.01,     0.0155},
  548.     { 2,    0,    0,    -3,      1.225,      0.91,     -0.0088},
  549.     { 1,    1,    0,    1,      1.267,      1.52,     0.0164},
  550.     { 1,    -1,    0,    -1,     -1.089,      0.55,      0},
  551.     { 0,    0,    2,    -1,      0.584,      8.84,     0.0071},
  552.     { 4,    0,    0,    0,      1.938,      3.60,     0.0401},
  553.     { 4,    0,    0,    -2,     -0.952,     -1.58,    -0.0130},
  554.     { 3,    1,    0,    0,     -0.551,      0.94,    -0.0097},
  555.     { 3,    1,    0,    -2,     -0.482,     -0.57,    -0.0045},
  556.     { 3,    -1,    0,    0,      0.681,      0.96,      0.0115},
  557.     { 2,    0,    2,    0,     -3.996,      0,     0.0004},
  558.     { 2,    0,    2,    -2,      0.557,     -0.75,     -0.0090},
  559.     { 2,    0,    -2,    2,     -0.459,     -0.38,    -0.0053},
  560.     { 2,    0,    -2,    0,     -1.298,      0.74,      0.0004},
  561.     { 2,    0,    -2,    -2,      0.538,      1.14,    -0.0141},
  562.     { 1,    1,    -2,    -2,      0.426,      0.07,        -0.0006},
  563.     { 1,    -1,    2,    0,     -0.304,      0.03,     0.0003},
  564.     { 1,    -1,    -2,    2,     -0.372,     -0.19,    -0.0027},
  565.     { 0,    0,    4,    0,      0.418,      0,     0},
  566.     { 2,    -1,    0,    -1,     -0.352,     -0.37,     -0.0028}
  567. };
  568.  
  569. # if MOON_TEST_CORR
  570. /* moon additional correction terms */
  571. struct m5dat {        
  572.     REAL8 lng;
  573.     int  i0,i1,i2,i3;
  574.     } m5 [] = {  
  575.     /*      lng,    l,     l',    F,    D,      */
  576.     0.127,    0,    0,    0,    6,
  577.     -0.151,    0,    2,    0,    -4,
  578.     -0.085,    0,    0,    2,    4,
  579.     0.150,    0,    1,    0,    3,
  580.     -0.091,    2,    1,    0,    -6,
  581.     -0.103,    0,    3,    0,    0,
  582.     -0.301,    1,    0,    2,    -4,
  583.     0.202,    1,    0,    -2,    -4,
  584.     0.137,    1,    1,    0,    -1,
  585.     0.233,    1,    1,    0,    -3,
  586.     -0.122,    1,    -1,    0,    1,
  587.     -0.276,    1,    -1,    0,    -3,
  588.     0.255,    0,    0,    2,    1,
  589.     0.254,    0,    0,    2,    -3,
  590.     -0.100,    3,    1,    0,    -4,
  591.     -0.183,    3,    -1,    0,    -2,
  592.     -0.297,    2,    2,    0,    -2,
  593.     -0.161,    2,    2,    0,    -4,
  594.     0.197,    2,    -2,    0,    0,
  595.     0.254,    2,    -2,    0,    -2,
  596.     -0.250,    1,    3,    0,    -2,
  597.     -0.123,    2,    0,    2,    2,
  598.     0.173,    2,    0,    -2,    -4,
  599.     0.263,    1,    1,    2,    0,
  600.     0.130,    3,    0,    0,    -1,
  601.     0.113,    5,    0,    0,    0,
  602.     0.092,    3,    0,    2,    -2,
  603.     0,    99,    0,    0,    0    /* end mark */
  604.     };
  605. # endif    /* MOON_TEST_CORR */
  606. #ifdef ASTROLOG
  607. /* End contents of helconst.c */
  608. #endif
  609. #ifdef ASTROLOG
  610. /* Begin contents of deltat.c */
  611. #endif
  612. /*****************************************************
  613. $Header: deltat.c,v 1.10 93/01/27 14:37:06 alois Exp $
  614. deltat.c
  615. deltat(t): returns delta t (in julian days) from universal time t
  616. is included by users
  617. ET = UT +  deltat
  618.  
  619.   ---------------------------------------------------------------
  620.   | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
  621.   | The use of this source code is subject to regulations made    |
  622.   | by Astrodienst Zurich. The code is NOT in the public domain.|
  623.   |                                |
  624.   | This copyright notice must not be changed or removed    |
  625.   | by any user of this program.                |
  626.   ---------------------------------------------------------------
  627.  
  628. ******************************************************/
  629.  
  630. #ifndef ASTROLOG
  631. double deltat (double jd_ad) /* Astrodienst relative julian date */
  632. #else
  633. double deltat(jd_ad)
  634. double jd_ad;
  635. #endif
  636. #ifndef ASTROLOG
  637.   double floor();
  638. #endif
  639.   static short int dt[] = {    /* in centiseconds */
  640.     /* dt from 1637 to 2000, as tabulated in A.E. 
  641.     the values 1620 - 1636 are not taken, as they fit
  642.     badly the parabola 25.5 t*t for the next range. The
  643.     best crossing point to switch to the parabola is
  644.     1637, where we have fitted the value for continuity */
  645.     6780,    6500,    6300,    
  646.     6200,    6000,    5800,    5700,    5500,
  647.     5400,    5300,    5100,    5000,    4900,
  648.     4800,    4700,    4600,    4500,    4400,
  649.     4300,    4200,    4100,    4000,    3800,    /* 1655 - 59 */
  650.     3700,    3600,    3500,    3400,    3300,
  651.     3200,    3100,    3000,    2800,    2700,
  652.     2600,    2500,    2400,    2300,    2200,
  653.     2100,    2000,    1900,    1800,    1700,
  654.     1600,    1500,    1400,    1400,    1300,
  655.     1200,    1200,    1100,    1100,    1000,
  656.     1000,    1000,    900,    900,    900,
  657.     900,    900,    900,    900,    900,
  658.     900,    900,    900,    900,    900,    /* 1700 - 1704 */
  659.     900,    900,    900,    1000,    1000,
  660.     1000,    1000,    1000,    1000,    1000,
  661.     1000,    1000,    1100,    1100,    1100,
  662.     1100,    1100,    1100,    1100,    1100,
  663.     1100,    1100,    1100,    1100,    1100,
  664.     1100,    1100,    1100,    1100,    1200,    /* 1730    - 1734 */
  665.     1200,    1200,    1200,    1200,    1200,
  666.     1200,    1200,    1200,    1200,    1300,
  667.     1300,    1300,    1300,    1300,    1300,
  668.     1300,    1400,    1400,    1400,    1400,
  669.     1400,    1400,    1400,    1500,    1500,
  670.     1500,    1500,    1500,    1500,    1500,    /* 1760 - 1764 */
  671.     1600,    1600,    1600,    1600,    1600,
  672.     1600,    1600,    1600,    1600,    1600,
  673.     1700,    1700,    1700,    1700,    1700,
  674.     1700,    1700,    1700,    1700,    1700,
  675.     1700,    1700,    1700,    1700,    1700,
  676.     1700,    1700,    1600,    1600,    1600,    /* 1790 - 1794 */
  677.     1600,    1500,    1500,    1400,    1400,
  678.     1370,    1340,    1310,    1290,    1270,    /* 1800 - 1804 */
  679.     1260,    1250,    1250,    1250,    1250,
  680.     1250,    1250,    1250,    1250,    1250,
  681.     1250,    1250,    1240,    1230,    1220,
  682.     1200,    1170,    1140,    1110,    1060,
  683.     1020,    960,    910,    860,    800,
  684.     750,    700,    660,    630,    600,    /* 1830 - 1834 */
  685.     580,    570,    560,    560,    560,
  686.     570,    580,    590,    610,    620,
  687.     630,    650,    660,    680,    690,
  688.     710,    720,    730,    740,    750,
  689.     760,    770,    770,    780,    780,
  690.     788,    782,    754,    697,    640,    /* 1860 - 1864 */
  691.     602,    541,    410,    292,    182,
  692.     161,    10,    -102,    -128,    -269,
  693.     -324,    -364,    -454,    -471,    -511,
  694.     -540,    -542,    -520,    -546,    -546,
  695.     -579,    -563,    -564,    -580,    -566,
  696.     -587,    -601,    -619,    -664,    -644,    /* 1890 - 1894 */
  697.     -647,    -609,    -576,    -466,    -374,
  698.     -272,    -154,    -2,    124,    264,
  699.     386,    537,    614,    775,    913,
  700.     1046,    1153,    1336,    1465,    1601,
  701.     1720,    1824,    1906,    2025,    2095,
  702.     2116,    2225,    2241,    2303,    2349,    /* 1920 - 1924 */
  703.     2362,    2386,    2449,    2434,    2408,
  704.     2402,    2400,    2387,    2395,    2386,
  705.     2393,    2373,    2392,    2396,    2402,
  706.     2433,    2483,    2530,    2570,    2624,
  707.     2677,    2728,    2778,    2825,    2871,
  708.     2915,    2957,    2997,    3036,    3072,    /* 1950 - 1954 */
  709.     3107,    3135,    3168,    3218,    3268,
  710.     3315,    3359,    3400,    3447,    3503,
  711.     3573,    3654,    3743,    3829,    3920,
  712.     4018,    4117,    4223,    4337,    4449,
  713.     4548,    4646,    4752,    4853,    4959,
  714.     5054,    5138,    5217,    5296,    5379,    /* 1980 - 1984 */
  715.     5434,    5487,    5532,    5582,    5630,    /* 1985 - 89 from AE 1991 */
  716.     5686,    5757,    5900,    5900,    6000,    /* AE 1993 and extrapol */
  717.     6050,    6100,    6150,    6200,    6250,    /* 1995 - 1999 */
  718.     6300};                    /* 2000 */
  719.   double yr, cy, delta;
  720.   long iyr, i;
  721.   yr = (jd_ad + 18262) / 365.25 + 100.0;    /*  year  relative 1800 */
  722.   cy = yr / 100;
  723.   iyr =  (long) (floor (yr) + 1800);     /* truncated to integer, rel 0 */
  724. # if TIDAL_26            /* Stephenson formula only when 26" tidal
  725.                    term in lunar motion */
  726.   if ( iyr >= 1637  && iyr < 2000 ) {
  727.     i = iyr - 1637;
  728.     delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - floor (yr)) * 0.01;
  729.   } else if (iyr >= 2000 ) {    /* parabola, fitted at value[2000] */
  730.     delta = 25.5 * cy * cy  - 25.5 * 4 + 63.00;    
  731.   } else if (iyr >= 948) {    /* from 948 - 1637 use parabola */
  732.     delta = 25.5 * cy * cy;    
  733.   } else {    /* before 984 use other parabola */    
  734.     delta = 1361.7    + 320 * cy + 44.3 * cy * cy;    /* fits at 948 */
  735.   }
  736. # else    /* use Clemence value + 5 sec before 1690, new dt afterwards */
  737.   cy -= 1;    /* epoch 1900 */
  738.   if ( iyr >= 1690  && iyr < 2000 ) {
  739.     i = iyr - 1637;
  740.     delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - floor (yr)) * 0.01;
  741.   } else if (iyr >= 2000 ) {    /* parabola, fitted at value[2000] */
  742.     delta = 29.949 * cy * cy  - 29.949 * 4 + 63.0;    
  743.   } else {    
  744.     delta = 5 + 24.349 + 72.3165 * cy + 29.949 * cy * cy; /* fits at 1690 */
  745.   }
  746. # endif
  747.   return (delta / 86400.0);
  748. }
  749. #ifdef ASTROLOG
  750. /* End contents of deltat.c */
  751. #endif
  752. #ifdef ASTROLOG
  753. /* Begin contents of d2l.c */
  754. #endif
  755. /*******************************************
  756. $Header: d2l.c,v 1.9 91/11/16 16:24:20 alois Exp $
  757. ********************************************/
  758.  
  759. /*************************************
  760. double to long with rounding, no overflow check
  761. *************************************/ 
  762. #ifndef ASTROLOG
  763. long d2l (double x)        
  764. #else
  765. long d2l(x)
  766. double x;
  767. #endif
  768. {
  769.   if (x >=0)
  770.     return ((long) (x + 0.5));
  771.   else
  772.     return (- (long) (0.5 - x));
  773. }
  774. #ifdef ASTROLOG
  775. /* End contents of d2l.c */
  776. #endif
  777. #ifdef ASTROLOG
  778. /* Begin contents of csec.c */
  779. #endif
  780. /*
  781.  * $Header$
  782.  *
  783.  * A collection of useful functions for centisec calculations.
  784.  
  785.   ---------------------------------------------------------------
  786.   | Copyright Astrodienst Zurich AG and Alois Treindl, 1991.    |
  787.   | The use of this source code is subject to regulations made    |
  788.   | by Astrodienst Zurich. The code is NOT in the public domain.|
  789.   |                                |
  790.   | This copyright notice must not be changed or removed    |
  791.   | by any user of this program.                |
  792.   ---------------------------------------------------------------
  793. *******************************************************/
  794. #ifndef ASTROLOG
  795. #include "ourdef.h"
  796. #include "astrolib.h"
  797. #include "housasp.h"
  798. #endif
  799.  
  800. /************************************
  801. normalize argument into interval [0..DEG360]
  802. *************************************/
  803. #ifndef ASTROLOG
  804. centisec csnorm(centisec p)
  805. {
  806.   if (p < 0) 
  807.     do { p += DEG360; } while (p < 0);
  808.   else if (p >= DEG360)
  809.     do { p -= DEG360; } while (p >= DEG360);
  810.   return (p);
  811. }
  812. #endif
  813.  
  814. #ifndef ASTROLOG
  815. double degnorm(double p)
  816. #else
  817. double degnorm(p)
  818. double p;
  819. #endif
  820. {
  821.   if (p < 0) 
  822.     do { p += 360.0; } while (p < 0);
  823.   else if (p >= 360.0)
  824.     do { p -= 360.0; } while (p >= 360.0);
  825.   return (p);
  826. }
  827.  
  828. #ifndef ASTROLOG
  829. /************************************
  830. distance in centisecs p1 - p2
  831. normalized to [0..360[
  832. **************************************/
  833. centisec difcsn (centisec p1, centisec p2)
  834.   return (csnorm(p1 - p2));
  835. }
  836.  
  837. double difdegn (double p1, double p2)
  838.   return (degnorm(p1 - p2));
  839. }
  840.  
  841. /************************************
  842. distance in centisecs p1 - p2
  843. normalized to [-180..180[
  844. **************************************/
  845. centisec difcs2n (centisec p1, centisec p2)
  846. { centisec dif;
  847.   dif = csnorm(p1 - p2);
  848.   if (dif  >= DEG180) return (dif - DEG360);
  849.   return (dif);
  850. }
  851.  
  852. double difdeg2n (double p1, double p2)
  853. { double dif;
  854.   dif = degnorm(p1 - p2);
  855.   if (dif  >= 180.0) return (dif - 360.0);
  856.   return (dif);
  857. }
  858.  
  859. /*************************************
  860. round second, but at 29.5959 always down
  861. *************************************/ 
  862. centisec roundsec(centisec x)    
  863. {
  864.   centisec t;
  865.   t = (x + 50) / 100 *100L;    /* round to seconds */
  866.   if (t > x && t % DEG30 == 0)  /* was rounded up to next sign */
  867.     t = x / 100 * 100L;        /* round last second of sign downwards */
  868.   return (t);
  869. }
  870.  
  871.  
  872. /******************************/
  873. double dcos(centisec x)
  874. {
  875.   return (COS8 (CSTORAD * x));
  876. }
  877.  
  878. /******************************/
  879. double dsin(centisec x)
  880. {
  881.   return (SIN8 (CSTORAD * x));
  882. }
  883.  
  884. /******************************/
  885. double dtan(centisec x)
  886. {
  887.   return (TAN8 (CSTORAD * x));
  888. }
  889.  
  890. /******************************/
  891. centisec datan(double x)
  892. {
  893.   return (d2l (RADTOCS * ATAN8 (x)) );    
  894. }
  895.  
  896. /******************************/
  897. centisec dasin(double x)
  898. {
  899.   return (d2l (RADTOCS * ASIN8 (x)));    
  900. }
  901. #endif /* ASTROLOG */
  902. #ifdef ASTROLOG
  903. /* End contents of csec.c */
  904. #endif
  905. #ifdef ASTROLOG
  906. /* Begin contents of julday.c */
  907. #endif
  908. /*********************************************************
  909.   $Header: julday.c,v 1.9 91/11/16 16:25:06 alois Exp $
  910. *********************************************************/
  911.  
  912. #ifndef ASTROLOG
  913. # include "ourdef.h"
  914. # include "astrolib.h"
  915. #endif
  916.  
  917. /*************** julday ********************************************
  918.  * This function returns the absolute Julian day number (JD)
  919.  * for a given calendar date.
  920.  * The aruguments are a calendar date: day, month, year as integers,
  921.  * hour as double with decimal fraction.
  922.  * If gregflag = 1, Gregorian calendar is assumed, gregflag = 0
  923.  * Julian calendar is assumed.
  924.  *
  925.  * The Julian day number is system of numbering all days continously
  926.  * within the time range of known human history. It should be familiar
  927.  * for every astrological or astronomical programmer. The time variable
  928.  * in astronomical theories is usually expressed in Julian days or
  929.  * Julian centuries (36525 days per century) relative to some start day;
  930.  * the start day is called 'the epoch'.
  931.  * The Julian day number is a double representing the number of
  932.  * days since JD = 0.0 on 1 Jan -4712, 12:00 noon.
  933.  * Midnight has always a JD with fraction .5, because traditionally
  934.  * the astronomical day started at noon.
  935.  *
  936.  * NOTE: The Julian day number is named after the monk Julianus. It must
  937.  * not be confused with the Julian calendar system, which is named after
  938.  * Julius Cesar, the Roman politician who introduced this calendar.
  939.  
  940.  * Original author: Marc Pottenger, Los Angeles.
  941.  * with bug fix for year < -4711   15-aug-88 by Alois Treindl
  942.  *
  943.  * References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
  944.  *             Verlag Sterne und Weltraum (1987), p.49 ff
  945.  *
  946.  * related functions: revjul() reverse Julian day number: compute the
  947.  *                   calendar date from a given JD
  948.  ****************************************************************/
  949. #ifndef ASTROLOG
  950. double julday(int month , int day , int year  , double hour  , int gregflag ) 
  951. #else
  952. double julday(month, day, year, hour, gregflag)
  953. int month;
  954. int day;
  955. int year;
  956. double hour;
  957. int gregflag;
  958. #endif
  959. {
  960.   double jd,u,u0,u1,u2,floor();
  961.   u = year;
  962.   if (month < 3) u -=1;
  963.   u0 = u + 4712.0;
  964.   u1 = month + 1.0;
  965.   if (u1 < 4) u1 += 12.0;
  966.   jd = floor(u0*365.25)
  967.      + floor(30.6*u1+0.000001)
  968.      + day + hour/24.0 - 63.5;
  969.   if (gregflag) {
  970.     u2 = floor(ABS8(u) / 100) - floor(ABS8(u) / 400);
  971.     if (u < 0.0) u2 = -u2;
  972.     jd = jd - u2 + 2;            
  973.     if ((u < 0.0) && (u/100 == floor(u/100)) && (u/400 != floor(u/400)))
  974.       jd -=1;
  975.   }
  976.   return(jd);
  977. }
  978.  
  979. #ifndef ASTROLOG
  980. /*
  981.  * monday = 0, ... sunday = 6
  982.  */
  983. int day_of_week(double jd)
  984. {
  985.   return (((int) floor (jd - 2433282 - 1.5) %7) + 7) % 7;
  986. }
  987. #endif
  988. #ifdef ASTROLOG
  989. /* End contents of julday.c */
  990. #endif
  991. #ifdef ASTROLOG
  992. /* Begin contents of revjul.c */
  993. #endif
  994. /*********************************************************
  995.   $Header: revjul.c,v 1.9 91/11/16 16:25:37 alois Exp $
  996. *********************************************************/
  997.  
  998. #ifndef ASTROLOG
  999. # include "ourdef.h"
  1000. # include "astrolib.h"
  1001. #endif
  1002.  
  1003. /*** revjul ******************************************************
  1004.   revjul() is the inverse function to julday(), see the description
  1005.   there.
  1006.   Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
  1007.   return values are the calendar day, month, year and the hour of
  1008.   the day with decimal fraction (0 .. 23.999999).
  1009.  
  1010.   Original author Mark Pottenger, Los Angeles.
  1011.   with bug fix for year < -4711 16-aug-88 Alois Treindl
  1012. *************************************************************************/
  1013. #ifndef ASTROLOG
  1014. void revjul (double jd, int gregflag,
  1015.          int *jmon, int *jday, int *jyear, double *jut)
  1016. #else
  1017. void revjul(jd, gregflag, jmon, jday, jyear, jut)
  1018. double jd;
  1019. int gregflag;
  1020. int *jmon;
  1021. int *jday;
  1022. int *jyear;
  1023. double *jut;
  1024. #endif
  1025. {
  1026.   double u0,u1,u2,u3,u4, floor();
  1027.   u0 = jd + 32082.5;
  1028.   if (gregflag) {
  1029.     u1 = u0 + floor (u0/36525.0) - floor (u0/146100.0) - 38.0;
  1030.     if (jd >= 1830691.5) u1 +=1;
  1031.     u0 = u0 + floor (u1/36525.0) - floor (u1/146100.0) - 38.0;
  1032.   }
  1033.   u2 = floor (u0 + 123.0);
  1034.   u3 = floor ( (u2 - 122.2) / 365.25);
  1035.   u4 = floor ( (u2 - floor (365.25 * u3) ) / 30.6001);
  1036. #ifndef ASTROLOG
  1037.   *jmon = u4-1.0;
  1038.   if (*jmon > 12) *jmon -= 12;
  1039.   *jday = u2 - floor (365.25 * u3) - floor (30.6001 * u4);
  1040.   *jyear = u3 + floor ( (u4 - 2.0) / 12.0) - 4800;
  1041. #else
  1042.   *jmon = (int)(u4-1.0);
  1043.   if (*jmon > 12) *jmon -= 12;
  1044.   *jday = (int)(u2 - floor (365.25 * u3) - floor (30.6001 * u4));
  1045.   *jyear = (int)(u3 + floor ( (u4 - 2.0) / 12.0) - 4800);
  1046. #endif
  1047.   *jut = (jd - floor (jd + 0.5) + 0.5) * 24.0;
  1048. }
  1049. #ifdef ASTROLOG
  1050. /* End contents of revjul.c */
  1051. #endif
  1052. #endif /* PLACALC */
  1053. #ifdef ASTROLOG
  1054. /* End contents of placalc2.c */
  1055. #endif
  1056.