home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume5 / rise_set < prev    next >
Text File  |  1986-11-30  |  14KB  |  502 lines
  1. Subject: Sun and Moon rise/set program
  2. Newsgroups: mod.sources
  3. Approved: jpn@panda.UUCP
  4.  
  5. Mod.sources:  Volume 5, Issue 10
  6. Submitted by: Marc Kaufman <talcott!su-shasta.arpa:kaufman>
  7.  
  8. Enclosed is a -C- program to compute sun and moon rise and set times.
  9. Have fun. The version below is the running source for a VAX.
  10. Send bugs (if any) to <kaufman@SU-Shasta> at Stanford.
  11.  
  12. ---------------------CUT HERE---------------------CUT HERE------------
  13. /* <sdate.c>
  14.  * Compute various useful times
  15.  *
  16.  * Written by Marc T. Kaufman
  17.  *            14100 Donelson Place
  18.  *            Los Altos Hills, CA 94022
  19.  *            (415) 948-3777
  20.  *
  21.  * Based on : "Explanatory Supplement to the Astronomical Ephemeris
  22.  *             and the American Ephemeris and Nautical Almanac",
  23.  *             H.M. Nautical Almanac Office, London.  Updated from
  24.  *             equations in the 1985 Astronomical Almanac.
  25.  *
  26.  * Copyright 1986 by Marc Kaufman
  27.  *
  28.  * Permission to use this program is granted, provided it is not sold.
  29.  *
  30.  * This program was originally written on a VAX, under 4.2bsd.
  31.  *  it was then ported to a 68000 system under REGULUS (Alcyon's version
  32.  *  of UNIX system III).  Major differences included: no 'double' and
  33.  *  a default integer length of 'short'.  Having been through all that,
  34.  *  porting to your machine should be easy.  Watch out for 'time' related
  35.  *  functions and make sure your 'atan2' program works right.
  36.  *
  37.  *    850210    revised to 1985 Ephemeris - mtk
  38.  */
  39.  
  40. #include <time.h>
  41. #include <sys/types.h>
  42. #include <sys/timeb.h>
  43. #include <stdio.h>
  44. #include <math.h>
  45.  
  46. long    UTC, TDT, tim, tim2, localdst;
  47. double    Julian_Day, MJD, Tu, Ru, T70, Local, GMST, LST;
  48. double    Eqt, Tua, L, G, e, eps, g, alpha, delta, sd, cd, lha, lhr, sh, ch;
  49. double    la, lf, S, C, sp, cp, tp, Az, alt;
  50. double    Lm, lm, px, SD, am, dm;
  51. double    zs, x;
  52. double    fabs(), fmod(), asin(), acos();
  53. struct    tm *t, *Rlocaltime(), *gmtime();
  54. char    *tdate, *gmctime(), *localctime();
  55. int    ftime();
  56. struct    timeb tb;
  57.  
  58. #define Pi            3.1415926535
  59. #define Degree_to_Radian    ((2.0 * Pi)/ 360.)
  60. #define Asec_Radian        ((2.0 * Pi)/(360. * 60. * 60.))
  61. #define Tsec_to_Radian        ((2.0 * Pi)/( 24. * 60.* 60.))
  62. #define Asec_to_Tsec        (24./360.)
  63. #define Sec_per_day        (24 * 60 * 60)
  64. #define Round            0.5        /* for rounding to integer */
  65.  
  66. #define J1900    /* 24 */15020.0    /* Julian Day number at Epoch 1900.0 */
  67. #define    J1970    /* 24 */40587.5    /* VAX clock Epoch 1970 Jan 1 (0h UT) */
  68. #define    J1985    /* 24 */46065.5    /* Epoch 1985 Jan 1 (0h UT) */
  69. #define    J2000    /* 24 */51545.0    /* Epoch 2000 Jan 1 (12h UT) */
  70. #define Delta_T        (54.6 + 0.9*(Julian_Day - J1985)/365.)    /* TDT - UT */
  71. /* (This is the position of my house ) */
  72. #define Longitude    (((122.)*60. +  8.)*60. +  3.)    /* Arc-seconds West */
  73. #define Latitude    ((( 37.)*60. + 22.)*60. + 58.)    /* Arc-seconds North */
  74. #define f1            (1. - (1./298.25))    /* 1 - flattening of Earth */
  75. /* the following alternate values are useful when debugging */
  76. /*#define Longitude    (((000.)*60. +  0.)*60. +  0.)    /* Arc-seconds West */
  77. /*#define Latitude    ((( 35.)*60. +  0.)*60. +  0.)    /* Arc-seconds North */
  78. /*#define f1        1.            /* 1 - flattening of Earth */
  79.  
  80. main() {
  81.  
  82. /* at this point we digress to discuss UNIX differences.
  83.  * In UCB UNIX we dont have ctime(), but do instead have asctime(),
  84.  *  which works from the structures created by gmtime() and localtime().
  85.  * However, system time is kept in UTC (Greenwich), and the localtime
  86.  *  routine correctly handles daylight savings time.
  87.  * Since the Regulus system only knows local time, a few direct
  88.  *  fiddles are needed.
  89.  */
  90.  
  91. /* correct apparent latitude for shape of Earth */
  92.  
  93.     lf= atan(f1*f1 * tan(Latitude * Asec_Radian));
  94.     sp= sin(lf);
  95.     cp= cos(lf);
  96.     tp= sp/cp;
  97.  
  98.     time(&UTC);            /* get time */
  99.     Local= - Longitude/15.;        /* Local apparent time correction */
  100.  
  101.     { int h, m, s;        /* manual entry mode */
  102. /*    time(&tim);
  103.     t= gmtime(&tim);
  104.     tim= tim - (60 * (60 * t->tm_hour + t->tm_min) + t->tm_se);
  105.     scanf("%d %d %d", &h, &m, &s);
  106.     {UTC = tim + 60 * (60 * h + m) + s;
  107. */    }
  108.  
  109. /* !    t= gmtime(&UTC);    /* this is Regulus time */
  110.     t= localtime(&UTC);    /* VAX version */
  111.  
  112. /* Compute delta to real UTC from time zone time */
  113.  
  114.     /* do this by hand since Regulus wont */
  115.  
  116.     switch (t->tm_mon + 1)    /* months are numbered from 0 */
  117.     {
  118.         case 1:
  119.         case 2:
  120.         case 3:
  121.         case 11:
  122.         case 12:
  123.             t->tm_isdst = 0;
  124.             break;
  125.  
  126.         case 5:
  127.         case 6:
  128.         case 7:
  129.         case 8:
  130.         case 9:
  131.             t->tm_isdst = 1;
  132.             break;
  133.  
  134.         case 4:
  135.             if ((t->tm_mday < 24) || (t->tm_mday - t->tm_wday <= 24))
  136.                 t->tm_isdst = 0;
  137.             else
  138.                 t->tm_isdst = 1;
  139.             break;
  140.  
  141.         case 10:
  142.             if ((t->tm_mday < 25) || (t->tm_mday - t->tm_wday <= 25))
  143.                 t->tm_isdst = 1;
  144.             else
  145.                 t->tm_isdst = 0;
  146.             break;
  147.     }
  148.     ftime(&tb);            /* gets time-zone information */
  149.     if (tb.dstflag == 0)
  150.         t->tm_isdst = 0;    /* dst never used here */
  151.     localdst = (-tb.timezone + t->tm_isdst*60) * 60L;    /* local time correction */
  152.  
  153. /* !    UTC -= localdst;    /* this is real UTC, not what the OS gave us! */
  154.     printf("%.24s GMT\n", gmctime(&UTC));
  155.  
  156.     stuff(UTC);            /* start with local time info */
  157.  
  158. /* Compute Terrestrial Dynamical Time (this used to be called Ephemeris Time) */
  159.  
  160.     TDT = UTC + (long)(Delta_T + Round);
  161.     tdate= gmctime(&TDT);
  162.     printf("           %.8s      Terrestrial Dynamical Time\n", tdate+11);
  163.  
  164.     printf("%.24s Local Civil Time\n", localctime(&UTC));
  165.  
  166.     tim2 = UTC + (long)(Local + Round);    /* Compute Local Solar Time */
  167.     tdate= gmctime(&tim2);
  168.     printf("           %.8s      Local Mean Time\n", tdate+11);
  169.  
  170. /* compute phase of moon */
  171.  
  172.     moondata(UTC);
  173.     Lm = fmod(Lm-L, 360.);    /* phase is Lm - L (longitude of Sun) */
  174.     lm = fmod(Lm, 90.);    /* excess over phase boundary */
  175.     printf("The Moon is%4.1f days past ", lm*36525./481267.883);
  176.     if    (Lm <  90.)    printf("New\n");
  177.     else if (Lm < 180.)    printf("First Quarter\n");
  178.     else if (Lm < 270.)    printf("Full\n");
  179.     else            printf("Last Quarter\n");
  180.  
  181.     printf("Julian Day  24%9.3f\n", Julian_Day);
  182.  
  183.     tim2 = GMST + Round;
  184.     tdate= gmctime(&tim2);
  185.         printf("           %.8s      Greenwich Mean Sidereal Time\n", tdate+11);
  186.  
  187.     tim2 = LST + Round;
  188.     tdate= gmctime(&tim2);
  189.         printf("           %.8s      Local Sidereal Time\n", tdate+11);
  190.  
  191.     tim2= lha + Round;
  192.     tdate= gmctime(&tim2);
  193.         printf("           %.8s      L.H.A. of Sun\n", tdate+11);
  194.         printf("            %11.3f  Degrees Declnation\n",delta/3600.);
  195.         printf("Azimuth     %11.3f  Degrees\n",Az/3600.);
  196.         printf("Elevation   %11.3f  Degrees\n",alt/3600.);
  197.  
  198.     /* compute sunrise and sunset */
  199.     t= Rlocaltime(&UTC);        /* compute start of day */
  200.     tim = UTC - (3600L * t->tm_hour + 60L * t->tm_min + t->tm_sec)
  201.             + Sec_per_day/2;    /* about noon */
  202.  
  203.     zs = 90. + 50./60.;            /* zenith angle of rise/set */
  204.     sunrise(tim, -1.0, zs, "Sunrise ");
  205.         printf("       ");
  206.         sunrise((long)(tim+Sec_per_day), -1.0, zs, "Tomorrow");
  207.         printf("\n");
  208.     sunrise(tim, 1.0, zs, "Sunset  ");
  209.         printf("       ");
  210.         sunrise((long)(tim+Sec_per_day), 1.0, zs, "Tomorrow");
  211.         printf("\n");
  212.  
  213.     /* compute moonrise and moonset */
  214.     tim = tim - Sec_per_day/2 - 31;    /* about start of day */
  215.  
  216.     zs = 90. + 34./60.;        /* zenith angle of rise/set */
  217.     moonrise(tim, -1.0, zs, "Moonrise");
  218.         printf("       ");
  219.         moonrise((long)(tim+Sec_per_day), -1.0, zs, "Tomorrow");
  220.         printf("\n");
  221.     moonrise(tim, 1.0, zs, "Moonset ");
  222.         printf("       ");
  223.         moonrise((long)(tim+Sec_per_day), 1.0, zs, "Tomorrow");
  224.         printf("\n");
  225. }
  226.  
  227. sunrise(t0, rs, z, s)
  228.     long t0;
  229.     double rs, z;
  230.     char *s;
  231. {
  232.     double cz, dh;
  233.     long dt;
  234.  
  235.     cz = cos(z * Degree_to_Radian);    /* zenith distance of phenomonon */
  236.  
  237.     do {    /* iterate */
  238.         stuff(t0);    /* compute declination and current hour angle */
  239.         dh= -tp*sd/cd + cz/(cp*cd);
  240.         if ((dh < -1.0) || (dh > 1.0)) {
  241.             printf("%.8s   none   ", s);
  242.             return;
  243.         }
  244.         dh=acos(dh)*rs;
  245.         dt= (dh - lhr) / Tsec_to_Radian;
  246.         t0 += dt;
  247.     } while (dt);
  248.  
  249.     t0 += 30 /* seconds, rounding to nearest minute */;
  250.     tdate= localctime(&t0);
  251.         printf("%.8s   %.5s  ", s, tdate+11);
  252. }
  253.  
  254. moonrise(t0, rs, z, s)
  255.     long t0;
  256.     double rs, z;
  257.     char *s;
  258. {
  259. #define SRATE    1.033863192    /* ratio of Moon's motion to Sun's motion */
  260.     double    cz, dh, sd, cd;
  261.     long    t1, dt;
  262.  
  263.     moondata(t0);    /* get starting declination of Moon */
  264.  
  265.     /* compute zenith distance of phenomonon */
  266.     cz = cos(z * Degree_to_Radian + SD /* -px */);
  267.  
  268.     /* first iteraton is forward only (to approx. phenom time) */
  269.     sd = sin(dm);
  270.     cd = cos(dm);
  271.     dh= -tp*sd/cd + cz/(cp*cd);
  272.     if ((dh < -1.0) || (dh > 1.0)) {
  273.         printf("%.8s   none   ", s);
  274.         return;
  275.     }
  276.     dh= acos(dh)*rs;
  277.     dt= fmod((dh - am), 2.0*Pi) * SRATE / Tsec_to_Radian;
  278.     t1 = t0 + dt;
  279.  
  280.     do {    /* iterate */
  281.         moondata(t1);    /* compute declination and current hour angle */
  282.         cz = cos(z * Degree_to_Radian + SD /* -px */);
  283.         sd = sin(dm);
  284.         cd = cos(dm);
  285.  
  286.         dh= -tp*sd/cd + cz/(cp*cd);
  287.         if ((dh < -1.0) || (dh > 1.0)) {
  288.             printf("%.8s   none  ", s);
  289.             return;
  290.         }
  291.         dh= acos(dh)*rs;
  292.         dt= (dh - am) * SRATE / Tsec_to_Radian;
  293.         t1 += dt;
  294.     } while (dt);
  295.  
  296.     if ((t1 - t0) >= Sec_per_day) {
  297.             printf("%.8s   none   ", s);
  298.         return;
  299.     }
  300.     t1 += 30 /* seconds, rounding to nearest minute */;
  301.     tdate= localctime(&t1);
  302.         printf("%.8s   %.5s  ", s, tdate+11);
  303. }
  304.  
  305. stuff(tim)
  306. long tim;
  307. {        /* main computation loop */
  308.  
  309.     timedata(tim);
  310.  
  311. /* where is the Sun (angles are in seconds of arc) */
  312. /*    Low precision elements from 1985 Almanac   */
  313.  
  314.     L= 280.460 + 0.9856474 * MJD;        /* Mean Longitde */
  315.     L = fmod(L, 360.);        /* corrected for aberration */
  316.  
  317.     g= 357.528 + 0.9856003 * MJD;        /* Mean Anomaly */
  318.     g = fmod(g, 360.);
  319.  
  320.     eps= 23.439 - 0.0000004 * MJD;        /* Mean Obiquity of Ecliptic */
  321.  
  322.     {    /* convert to R.A. and DEC */
  323.         double Lr, gr, epsr, lr, ca, sa, R;
  324.         double sA, cA, gphi;
  325.  
  326.         Lr = L * Degree_to_Radian;
  327.         gr = g * Degree_to_Radian;
  328.         epsr = eps * Degree_to_Radian;
  329.  
  330.         lr = (L + 1.915*sin(gr) + 0.020*sin(2.0*gr)) * Degree_to_Radian;
  331.  
  332.         sd = sin(lr) * sin(epsr);
  333.         cd = sqrt(1.0 - sd*sd);
  334.         sa = sin(lr) * cos(epsr);
  335.         ca = cos(lr);
  336.  
  337.         delta = asin(sd);
  338.         alpha = atan2(sa, ca);
  339.  
  340.     /* equation of time */
  341.         Eqt= (Lr - alpha) / Tsec_to_Radian;
  342.  
  343.         delta = delta / Asec_Radian;
  344.         alpha = alpha / Tsec_to_Radian;
  345.  
  346.         lhr = (LST - alpha) * Tsec_to_Radian;
  347.         sh =  sin(lhr);
  348.         ch =  cos(lhr);
  349.         lhr= atan2(sh, ch);    /* normalized -pi to pi */
  350.         lha= lhr / Tsec_to_Radian + Sec_per_day/2;
  351.  
  352.     /* convert to Azimuth and altitude */
  353.  
  354.         alt = asin(sd*sp + cd*ch*cp);
  355.         ca =  cos(alt);
  356.         sA =  -cd * sh / ca;
  357.         cA =  (sd*cp - cd*ch*sp) / ca;
  358.         Az = atan2(sA, cA) / Asec_Radian;
  359.         Az = fmod(Az, 1296000. /* 360.*3600. */);
  360.         alt = alt / Asec_Radian;
  361.     }
  362. }
  363.  
  364. moondata(tim)
  365. long    tim;
  366. {
  367.     double    lst, beta, rm, sa, ca, sl, cl, sb, cb, x, y, z, l, m, n;
  368.  
  369. /* compute location of the moon */
  370. /* Ephemeris elements from 1985 Almanac */
  371.  
  372.     timedata(tim);
  373.  
  374.     Lm= 218.32 + 481267.883*Tu
  375.         + 6.29 * sin((134.9 + 477198.85*Tu)*Degree_to_Radian)
  376.         - 1.27 * sin((259.2 - 413335.38*Tu)*Degree_to_Radian)
  377.         + 0.66 * sin((235.7 + 890534.23*Tu)*Degree_to_Radian)
  378.         + 0.21 * sin((269.9 + 954397.70*Tu)*Degree_to_Radian)
  379.         - 0.19 * sin((357.5 +  35999.05*Tu)*Degree_to_Radian)
  380.         - 0.11 * sin((186.6 + 966404.05*Tu)*Degree_to_Radian);
  381.  
  382.     beta=      5.13 * sin(( 93.3 + 483202.03*Tu)*Degree_to_Radian)
  383.         + 0.28 * sin((228.2 + 960400.87*Tu)*Degree_to_Radian)
  384.         - 0.28 * sin((318.3 +   6003.18*Tu)*Degree_to_Radian)
  385.         - 0.17 * sin((217.6 - 407332.20*Tu)*Degree_to_Radian);
  386.  
  387.     px= 0.9508
  388.         + 0.0518 * cos((134.9 + 477198.85*Tu)*Degree_to_Radian)
  389.         + 0.0095 * cos((259.2 - 413335.38*Tu)*Degree_to_Radian)
  390.         + 0.0078 * cos((235.7 + 890534.23*Tu)*Degree_to_Radian)
  391.         + 0.0028 * cos((269.9 + 954397.70*Tu)*Degree_to_Radian);
  392.  
  393. /*    SD= 0.2725 * px;    */
  394.  
  395.     rm= 1.0 / sin(px * Degree_to_Radian);
  396.  
  397.     lst= (100.46 + 36000.77*Tu) * Degree_to_Radian
  398.         + ((tim % Sec_per_day) + Local) * Tsec_to_Radian;
  399.  
  400. /* form geocentric direction cosines */
  401.  
  402.     sl= sin(Lm * Degree_to_Radian);
  403.     cl= cos(Lm * Degree_to_Radian);
  404.     sb= sin(beta* Degree_to_Radian);
  405.     cb= cos(beta * Degree_to_Radian);
  406.  
  407.     l= cb * cl;
  408.     m= 0.9175 * cb * sl - 0.3978 * sb;
  409.     n= 0.3978 * cb * sl + 0.9175 * sb;
  410.  
  411. /* R.A. and Dec of Moon, geocentric*/
  412.  
  413.     am= atan2(m, l);
  414.     dm= asin(n);
  415.  
  416. /* topocentric rectangular coordinates */
  417.  
  418.     cd= cos(dm);
  419.     sd= n;
  420.     ca= cos(am);
  421.     sa= sin(am);
  422.     sl= sin(lst);
  423.     cl= cos(lst);
  424.  
  425.     x= rm * cd *ca - cp * cl;
  426.     y= rm * cd * sa - cp * sl;
  427.     z= rm * sd - sp;
  428.  
  429. /* finally, topocentric Hour-Angle and Dec */
  430.  
  431.     am = lst - atan2(y, x);
  432.     ca = cos(am);
  433.     sa = sin(am);
  434.     am = atan2(sa,ca);
  435.     rm = sqrt(x*x + y*y + z*z);
  436.     dm = asin(z/rm);
  437.     px = asin(1.0 / rm);
  438.     SD = 0.2725 * px;
  439. }
  440.  
  441. timedata(tim)
  442. long    tim;
  443. {
  444.  
  445. /* compute seconds from 2000 Jan 1.5 UT (Ephemeris Epoch) */
  446. /* the VAX Epoch is     1970 Jan 1.0 UT (Midnight on Jan 1) */
  447.  
  448.     Julian_Day = (tim/Sec_per_day) +
  449.                 (double)(tim % Sec_per_day)/Sec_per_day + J1970;
  450.     MJD= Julian_Day -J2000;    /* Julian Days past Epoch */
  451.     Tu = MJD/36525.;        /* Julian Centuries past Epoch */
  452.  
  453. /* compute Sidereal time */
  454.  
  455.     Ru= 24110.54841 + Tu * (8640184.812866
  456.         + Tu * (0.09304 - Tu * 6.2e-6));    /* seconds */
  457.     GMST = (tim % Sec_per_day) + Sec_per_day + fmod(Ru, (double)Sec_per_day);
  458.     LST  = GMST + Local;
  459. }
  460.  
  461. /* time functions, for Regulus */
  462. char *gmctime(t)    /* re-hack for VAX, since ctime gives local */
  463. long *t;
  464. {
  465.     long t1;
  466.  
  467.     t1 = *t - localdst;        /* convert to local time */
  468.     return(ctime(&t1));
  469. }
  470.  
  471. char *localctime(t)
  472. long *t;
  473. {
  474.     long t1;
  475.  
  476.     t1 = *t + localdst;        /* convert to local time */
  477.     return(gmctime(&t1));
  478. }
  479.  
  480. struct tm *Rlocaltime(t)
  481. long *t;
  482. {
  483.     long t1;
  484.  
  485.     t1 = *t + localdst;        /* convert to local time */
  486.     return(gmtime(&t1));
  487. }
  488.  
  489. /* double precision modulus, put in range 0 <= result < m */
  490. double fmod(x, m)
  491. double    x, m;
  492. {
  493.     long i;
  494.  
  495.     i = fabs(x)/m;        /* compute integer part of x/m */
  496.     if (x < 0)    return( x + (i+1)*m);
  497.     else        return( x - i*m);
  498. }
  499. %%
  500.  
  501.  
  502.