home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / sun / volume1 / calentool / part07 / riseset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-05-27  |  12.9 KB  |  444 lines
  1. /*
  2.  * $Header: riseset.c,v 2.1 89/05/09 14:19:49 billr Exp $
  3.  */
  4. /*
  5.  * riseset.c
  6.  *
  7.  * Compute various useful times
  8.  *
  9.  * Written by Marc T. Kaufman
  10.  *            14100 Donelson Place
  11.  *            Los Altos Hills, CA 94022
  12.  *            (415) 948-3777
  13.  *
  14.  * Based on : "Explanatory Supplement to the Astronomical Ephemeris
  15.  *             and the American Ephemeris and Nautical Almanac",
  16.  *             H.M. Nautical Almanac Office, London.  Updated from
  17.  *             equations in the 1985 Astronomical Almanac.
  18.  *
  19.  * Copyright 1986 by Marc Kaufman
  20.  *
  21.  * Permission to use this program is granted, provided it is not sold.
  22.  *
  23.  * This program was originally written on a VAX, under 4.2bsd.
  24.  *  it was then ported to a 68000 system under REGULUS (Alcyon's version
  25.  *  of UNIX system III).  Major differences included: no 'double' and
  26.  *  a default integer length of 'short'.  Having been through all that,
  27.  *  porting to your machine should be easy.  Watch out for 'time' related
  28.  *  functions and make sure your 'atan2' program works right.
  29.  *
  30.  *    850210    revised to 1985 Ephemeris - mtk
  31.  *
  32.  * Modified 1/26/89 by Bill Randle, Tektronix, Inc., to interface
  33.  * to the calentool program. Also removed some of the REGULUS stuff.
  34.  * Modifications Copyright (C) 1989, Tektronix, Inc.  All Rights Reserved
  35.  *
  36.  * Permission is hereby granted to use and modify this code in source
  37.  * or binary form as long as it is not sold for profit and this copyright
  38.  * notice remains intact.
  39.  */
  40. #include "ct.h"        /* for the NO_SUN_MOON #define */
  41. #ifndef NO_SUN_MOON
  42.  
  43. #include <stdio.h>
  44. #include <sys/time.h>
  45. #include <math.h>
  46. #include "riseset.h"
  47.  
  48. extern char *daynames[], *monthnames[];
  49.  
  50. long    UTC, TDT, tim, tim2;
  51. double    Julian_Day, MJD, Tu, Ru, T70, Local, GMST, LST;
  52. double    Eqt, Tua, L, G, e, eps, g, alpha, delta, sd, cd, lha, lhr, sh, ch;
  53. double    la, lf, S, C, sp, cp, tp, Az, alt;
  54. double    Lm, lm, px, SD, am, dm;
  55. double    zs, x;
  56. double    fabs(), fmod(), asin(), acos();
  57. struct    tm *ltm, *t, *gmtime(), *localtime();
  58. char    *tdate, *gmctime(), *asctime(), *timezone();
  59. struct    timeval tvp;
  60. struct    timezone tzp;
  61. char    riseset_buf[B_SIZE][64];
  62. char    rbuf[24];
  63.  
  64. #define Pi            3.1415926535
  65. #define Degree_to_Radian    ((2.0 * Pi)/ 360.)
  66. #define Asec_Radian        ((2.0 * Pi)/(360. * 60. * 60.))
  67. #define Tsec_to_Radian        ((2.0 * Pi)/( 24. * 60.* 60.))
  68. #define Asec_to_Tsec        (24./360.)
  69. #define Sec_per_day        (24 * 60 * 60)
  70. #define Round            0.5        /* for rounding to integer */
  71.  
  72. #define J1900    /* 24 */15020.0    /* Julian Day number at Epoch 1900.0 */
  73. #define    J1970    /* 24 */40587.5    /* VAX clock Epoch 1970 Jan 1 (0h UT) */
  74. #define    J1985    /* 24 */46065.5    /* Epoch 1985 Jan 1 (0h UT) */
  75. #define    J2000    /* 24 */51545.0    /* Epoch 2000 Jan 1 (12h UT) */
  76. #define Delta_T        (54.6 + 0.9*(Julian_Day - J1985)/365.)    /* TDT - UT */
  77. /* (This is the position of my house ) */
  78. #ifndef LONG_DEGREES
  79. #    define Longitude    0.
  80. #else
  81. #    define Longitude    (((LONG_DEGREES)*60. +  LONG_MINUTES)*60. +  LONG_SECONDS)    /* Arc-seconds West */
  82. #endif
  83. #ifndef LAT_DEGREES
  84. #    define Latitude        37.*360.
  85. #else
  86. #    define Latitude    (((LAT_DEGREES)*60. + LAT_MINUTES)*60. + LAT_SECONDS)    /* Arc-seconds North */
  87. #endif
  88. #define f1            (1. - (1./298.25))    /* 1 - flattening of Earth */
  89. /* the following alternate values are useful when debugging */
  90. /*#define Longitude    (((000.)*60. +  0.)*60. +  0.)    /* Arc-seconds West */
  91. /*#define Latitude    ((( 35.)*60. +  0.)*60. +  0.)    /* Arc-seconds North */
  92. /*#define f1        1.            /* 1 - flattening of Earth */
  93.  
  94. riseset(epoch_seconds)
  95. long epoch_seconds;
  96. {
  97.     char *moonrise(), *sunrise();
  98.  
  99. /* at this point we digress to discuss UNIX differences.
  100.  * In UCB UNIX we dont have ctime(), but do instead have asctime(),
  101.  *  which works from the structures created by gmtime() and localtime().
  102.  * However, system time is kept in UTC (Greenwich), and the localtime
  103.  *  routine correctly handles daylight savings time.
  104.  * Since the Regulus system only knows local time, a few direct
  105.  *  fiddles are needed.
  106.  */
  107.  
  108. /* correct apparent latitude for shape of Earth */
  109.  
  110.     lf= atan(f1*f1 * tan(Latitude * Asec_Radian));
  111.     sp= sin(lf);
  112.     cp= cos(lf);
  113.     tp= sp/cp;
  114.  
  115.     UTC = epoch_seconds;
  116.     gettimeofday(&tvp, &tzp);    /* for timezone info */
  117. #ifdef LONG_DEGREES
  118.     Local = - Longitude/15.;    /* Local apparent time correction */
  119. #else
  120.     Local = (double)(-tzp.minuteswest * 60.);
  121. #endif
  122.  
  123.     sprintf(riseset_buf[B_GMT], "%.24s GMT", gmctime(&UTC));
  124.  
  125.     stuff(UTC);            /* start with local time info */
  126.  
  127. /* Compute Terrestrial Dynamical Time (this used to be called Ephemeris Time) */
  128.  
  129.     TDT = UTC + (long)(Delta_T + Round);
  130.     tdate= gmctime(&TDT);
  131.     *(tdate+19) = '\0';
  132.     sprintf(riseset_buf[B_TDT], "           %.8s      Terrestrial Dynamical Time", tdate+11);
  133.  
  134.     ltm = localtime(&UTC);
  135.     sprintf(riseset_buf[B_LCT], "%.24s Local Civil Time", asctime(ltm));
  136.  
  137.     sprintf(riseset_buf[B_DMY], "%s %s %d, %d", daynames[ltm->tm_wday],
  138.         monthnames[ltm->tm_mon], ltm->tm_mday, ltm->tm_year+1900);
  139.     tim2 = UTC + (long)(Local + Round);    /* Compute Local Solar Time */
  140.     tdate= gmctime(&tim2);
  141.     *(tdate+19) = '\0';
  142.     sprintf(riseset_buf[B_LMT], "           %.8s      Local Mean Time", tdate+11);
  143.  
  144. /* compute phase of moon */
  145.  
  146.     moondata(UTC);
  147.     Lm = fmod(Lm-L, 360.);    /* phase is Lm - L (longitude of Sun) */
  148.     lm = fmod(Lm, 90.);    /* excess over phase boundary */
  149.     sprintf(riseset_buf[B_POM], "The Moon is %3.1f days past ", lm*36525./481267.883);
  150.     if    (Lm <  90.)    strcat(riseset_buf[B_POM], "New");
  151.     else if (Lm < 180.)    strcat(riseset_buf[B_POM], "First Quarter");
  152.     else if (Lm < 270.)    strcat(riseset_buf[B_POM], "Full");
  153.     else            strcat(riseset_buf[B_POM], "Last Quarter");
  154.  
  155.     sprintf(riseset_buf[B_JLD], "Julian Day  24%10.4f", Julian_Day);
  156.  
  157.     tim2 = GMST + Round;
  158.     tdate= gmctime(&tim2);
  159.     *(tdate+19) = '\0';
  160.     sprintf(riseset_buf[B_GST], "           %.8s      Greenwich Mean Sidereal Time", tdate+11);
  161.  
  162.     tim2 = LST + Round;
  163.     tdate= gmctime(&tim2);
  164.     *(tdate+19) = '\0';
  165.     sprintf(riseset_buf[B_LST], "           %.8s      Local Sidereal Time", tdate+11);
  166.  
  167.     tim2= lha + Round;
  168.     tdate= gmctime(&tim2);
  169.     *(tdate+19) = '\0';
  170.     sprintf(riseset_buf[B_LHA], "L.H.A. of Sun  %.8s", tdate+11);
  171.     sprintf(riseset_buf[B_SDE], "   Declination  %6.3f  Degrees",delta/3600.);
  172.     sprintf(riseset_buf[B_SAZ], "   Azimuth      %6.3f  Degrees",Az/3600.);
  173.     sprintf(riseset_buf[B_SEL], "   Elevation     %6.3f  Degrees",alt/3600.);
  174.  
  175.     /* compute sunrise and sunset */
  176.     t= localtime(&UTC);        /* compute start of day */
  177.     tim = UTC - (3600L * t->tm_hour + 60L * t->tm_min + t->tm_sec)
  178.             + Sec_per_day/2;    /* about noon */
  179.  
  180.     zs = 90. + 50./60.;        /* zenith angle of rise/set */
  181.     strcpy(riseset_buf[B_SRD], sunrise(tim, -1.0, zs));
  182.     strcpy(riseset_buf[B_SSD], sunrise(tim, 1.0, zs));
  183.     strcpy(riseset_buf[B_SRT], sunrise((long)(tim+Sec_per_day), -1.0, zs));
  184.     strcpy(riseset_buf[B_SST], sunrise((long)(tim+Sec_per_day), 1.0, zs));
  185.  
  186.     /* compute moonrise and moonset */
  187.     tim = tim - Sec_per_day/2 - 31;    /* about start of day */
  188.  
  189.     zs = 90. + 34./60.;        /* zenith angle of rise/set */
  190.     strcpy(riseset_buf[B_MRD], moonrise(tim, -1.0, zs));
  191.     strcpy(riseset_buf[B_MSD], moonrise(tim, 1.0, zs));
  192.     strcpy(riseset_buf[B_MRT], moonrise((long)(tim+Sec_per_day), -1.0, zs));
  193.     strcpy(riseset_buf[B_MST], moonrise((long)(tim+Sec_per_day), 1.0, zs));
  194. }
  195.  
  196. char *
  197. sunrise(t0, rs, z)
  198.     long t0;
  199.     double rs, z;
  200. {
  201.     double cz, dh;
  202.     long dt;
  203.     struct tm *ltm;
  204.  
  205.     cz = cos(z * Degree_to_Radian);    /* zenith distance of phenomonon */
  206.  
  207.     do {    /* iterate */
  208.         stuff(t0);    /* compute declination and current hour angle */
  209.         dh= -tp*sd/cd + cz/(cp*cd);
  210.         if ((dh < -1.0) || (dh > 1.0))
  211.             return("none");
  212.         dh=acos(dh)*rs;
  213.         dt= (dh - lhr) / Tsec_to_Radian;
  214.         t0 += dt;
  215.     } while (dt);
  216.  
  217.     t0 += 30 /* seconds, rounding to nearest minute */;
  218.     ltm = localtime(&t0);
  219.     sprintf(rbuf, "%02d:%02d %s", ltm->tm_hour, ltm->tm_min, timezone(tzp.tz_minuteswest, ltm->tm_isdst));
  220.     return(rbuf);
  221. }
  222.  
  223. char *
  224. moonrise(t0, rs, z)
  225.     long t0;
  226.     double rs, z;
  227. {
  228. #define SRATE    1.033863192    /* ratio of Moon's motion to Sun's motion */
  229.     double    cz, dh, sd, cd;
  230.     long    t1, dt;
  231.     struct tm *ltm;
  232.  
  233.     moondata(t0);    /* get starting declination of Moon */
  234.  
  235.     /* compute zenith distance of phenomonon */
  236.     cz = cos(z * Degree_to_Radian + SD /* -px */);
  237.  
  238.     /* first iteraton is forward only (to approx. phenom time) */
  239.     sd = sin(dm);
  240.     cd = cos(dm);
  241.     dh= -tp*sd/cd + cz/(cp*cd);
  242.     if ((dh < -1.0) || (dh > 1.0))
  243.         return("none");
  244.     dh= acos(dh)*rs;
  245.     dt= fmod((dh - am), 2.0*Pi) * SRATE / Tsec_to_Radian;
  246.     t1 = t0 + dt;
  247.  
  248.     do {    /* iterate */
  249.         moondata(t1);    /* compute declination and current hour angle */
  250.         cz = cos(z * Degree_to_Radian + SD /* -px */);
  251.         sd = sin(dm);
  252.         cd = cos(dm);
  253.  
  254.         dh= -tp*sd/cd + cz/(cp*cd);
  255.         if ((dh < -1.0) || (dh > 1.0))
  256.             return("none");
  257.         dh= acos(dh)*rs;
  258.         dt= (dh - am) * SRATE / Tsec_to_Radian;
  259.         t1 += dt;
  260.     } while (dt);
  261.  
  262.     if ((t1 - t0) >= Sec_per_day)
  263.         return("none");
  264.     t1 += 30 /* seconds, rounding to nearest minute */;
  265.     ltm = localtime(&t1);
  266.     sprintf(rbuf, "%02d:%02d %s", ltm->tm_hour, ltm->tm_min, timezone(tzp.tz_minuteswest, ltm->tm_isdst));
  267.     return(rbuf);
  268. }
  269.  
  270. stuff(tim)
  271. long tim;
  272. {        /* main computation loop */
  273.  
  274.     timedata(tim);
  275.  
  276. /* where is the Sun (angles are in seconds of arc) */
  277. /*    Low precision elements from 1985 Almanac   */
  278.  
  279.     L= 280.460 + 0.9856474 * MJD;        /* Mean Longitde */
  280.     L = fmod(L, 360.);        /* corrected for aberration */
  281.  
  282.     g= 357.528 + 0.9856003 * MJD;        /* Mean Anomaly */
  283.     g = fmod(g, 360.);
  284.  
  285.     eps= 23.439 - 0.0000004 * MJD;        /* Mean Obiquity of Ecliptic */
  286.  
  287.     {    /* convert to R.A. and DEC */
  288.         double Lr, gr, epsr, lr, ca, sa, R;
  289.         double sA, cA, gphi;
  290.  
  291.         Lr = L * Degree_to_Radian;
  292.         gr = g * Degree_to_Radian;
  293.         epsr = eps * Degree_to_Radian;
  294.  
  295.         lr = (L + 1.915*sin(gr) + 0.020*sin(2.0*gr)) * Degree_to_Radian;
  296.  
  297.         sd = sin(lr) * sin(epsr);
  298.         cd = sqrt(1.0 - sd*sd);
  299.         sa = sin(lr) * cos(epsr);
  300.         ca = cos(lr);
  301.  
  302.         delta = asin(sd);
  303.         alpha = atan2(sa, ca);
  304.  
  305.     /* equation of time */
  306.         Eqt= (Lr - alpha) / Tsec_to_Radian;
  307.  
  308.         delta = delta / Asec_Radian;
  309.         alpha = alpha / Tsec_to_Radian;
  310.  
  311.         lhr = (LST - alpha) * Tsec_to_Radian;
  312.         sh =  sin(lhr);
  313.         ch =  cos(lhr);
  314.         lhr= atan2(sh, ch);    /* normalized -pi to pi */
  315.         lha= lhr / Tsec_to_Radian + Sec_per_day/2;
  316.  
  317.     /* convert to Azimuth and altitude */
  318.  
  319.         alt = asin(sd*sp + cd*ch*cp);
  320.         ca =  cos(alt);
  321.         sA =  -cd * sh / ca;
  322.         cA =  (sd*cp - cd*ch*sp) / ca;
  323.         Az = atan2(sA, cA) / Asec_Radian;
  324.         Az = fmod(Az, 1296000. /* 360.*3600. */);
  325.         alt = alt / Asec_Radian;
  326.     }
  327. }
  328.  
  329. moondata(tim)
  330. long    tim;
  331. {
  332.     double    lst, beta, rm, sa, ca, sl, cl, sb, cb, x, y, z, l, m, n;
  333.  
  334. /* compute location of the moon */
  335. /* Ephemeris elements from 1985 Almanac */
  336.  
  337.     timedata(tim);
  338.  
  339.     Lm= 218.32 + 481267.883*Tu
  340.         + 6.29 * sin((134.9 + 477198.85*Tu)*Degree_to_Radian)
  341.         - 1.27 * sin((259.2 - 413335.38*Tu)*Degree_to_Radian)
  342.         + 0.66 * sin((235.7 + 890534.23*Tu)*Degree_to_Radian)
  343.         + 0.21 * sin((269.9 + 954397.70*Tu)*Degree_to_Radian)
  344.         - 0.19 * sin((357.5 +  35999.05*Tu)*Degree_to_Radian)
  345.         - 0.11 * sin((186.6 + 966404.05*Tu)*Degree_to_Radian);
  346.  
  347.     beta=      5.13 * sin(( 93.3 + 483202.03*Tu)*Degree_to_Radian)
  348.         + 0.28 * sin((228.2 + 960400.87*Tu)*Degree_to_Radian)
  349.         - 0.28 * sin((318.3 +   6003.18*Tu)*Degree_to_Radian)
  350.         - 0.17 * sin((217.6 - 407332.20*Tu)*Degree_to_Radian);
  351.  
  352.     px= 0.9508
  353.         + 0.0518 * cos((134.9 + 477198.85*Tu)*Degree_to_Radian)
  354.         + 0.0095 * cos((259.2 - 413335.38*Tu)*Degree_to_Radian)
  355.         + 0.0078 * cos((235.7 + 890534.23*Tu)*Degree_to_Radian)
  356.         + 0.0028 * cos((269.9 + 954397.70*Tu)*Degree_to_Radian);
  357.  
  358. /*    SD= 0.2725 * px;    */
  359.  
  360.     rm= 1.0 / sin(px * Degree_to_Radian);
  361.  
  362.     lst= (100.46 + 36000.77*Tu) * Degree_to_Radian
  363.         + ((tim % Sec_per_day) + Local) * Tsec_to_Radian;
  364.  
  365. /* form geocentric direction cosines */
  366.  
  367.     sl= sin(Lm * Degree_to_Radian);
  368.     cl= cos(Lm * Degree_to_Radian);
  369.     sb= sin(beta* Degree_to_Radian);
  370.     cb= cos(beta * Degree_to_Radian);
  371.  
  372.     l= cb * cl;
  373.     m= 0.9175 * cb * sl - 0.3978 * sb;
  374.     n= 0.3978 * cb * sl + 0.9175 * sb;
  375.  
  376. /* R.A. and Dec of Moon, geocentric*/
  377.  
  378.     am= atan2(m, l);
  379.     dm= asin(n);
  380.  
  381. /* topocentric rectangular coordinates */
  382.  
  383.     cd= cos(dm);
  384.     sd= n;
  385.     ca= cos(am);
  386.     sa= sin(am);
  387.     sl= sin(lst);
  388.     cl= cos(lst);
  389.  
  390.     x= rm * cd *ca - cp * cl;
  391.     y= rm * cd * sa - cp * sl;
  392.     z= rm * sd - sp;
  393.  
  394. /* finally, topocentric Hour-Angle and Dec */
  395.  
  396.     am = lst - atan2(y, x);
  397.     ca = cos(am);
  398.     sa = sin(am);
  399.     am = atan2(sa,ca);
  400.     rm = sqrt(x*x + y*y + z*z);
  401.     dm = asin(z/rm);
  402.     px = asin(1.0 / rm);
  403.     SD = 0.2725 * px;
  404. }
  405.  
  406. timedata(tim)
  407. long    tim;
  408. {
  409.  
  410. /* compute seconds from 2000 Jan 1.5 UT (Ephemeris Epoch) */
  411. /* the VAX Epoch is     1970 Jan 1.0 UT (Midnight on Jan 1) */
  412.  
  413.     Julian_Day = (tim/Sec_per_day) +
  414.                 (double)(tim % Sec_per_day)/Sec_per_day + J1970;
  415.     MJD= Julian_Day -J2000;    /* Julian Days past Epoch */
  416.     Tu = MJD/36525.;        /* Julian Centuries past Epoch */
  417.  
  418. /* compute Sidereal time */
  419.  
  420.     Ru= 24110.54841 + Tu * (8640184.812866
  421.         + Tu * (0.09304 - Tu * 6.2e-6));    /* seconds */
  422.     GMST = (tim % Sec_per_day) + Sec_per_day + fmod(Ru, (double)Sec_per_day);
  423.     LST  = GMST + Local;
  424. }
  425.  
  426. /* time functions */
  427. char *gmctime(t)
  428. long *t;
  429. {
  430.     return(asctime(gmtime(t)));
  431. }
  432.  
  433. /* double precision modulus, put in range 0 <= result < m */
  434. double fmod(x, m)
  435. double    x, m;
  436. {
  437.     long i;
  438.  
  439.     i = fabs(x)/m;        /* compute integer part of x/m */
  440.     if (x < 0)    return( x + (i+1)*m);
  441.     else        return( x - i*m);
  442. }
  443. #endif    /* NO_SUN_MOON */
  444.