source: buchla-68k/orig/DATE/TRANSIT.C@ d21fc6f

/*
   ============================================================================ 
        transit.c -- calculate sunrise, sunset, and hours of daylight
        Version 2 -- 1989-01-26 -- Robert Bond / D.N. Lynx Crowe
 
        transit <options>

        options:        -d yyyy-mm-dd   first date (default is current system date)
                        -e yyyy-mm-dd   last date (default is current system date)
                        -a lat          decimal latitude (default = 45.5333)
                        -o lon          decimal longitude (default = 122.8333) 
                        -y yyyy         year  (calculates for the entire year)
                        -z tz           timezone (default = 8, PST)
        
        Calculated output is written to OUTFILE, errors are written to stdout.  

        Define BSD for BSD4.2 systems.  System V assumed otherwise.
        This selects between <sys/time.h> and <time.h> below.

        Note that the latitude, longitude, time zone correction and
        time zone string are all defaulted in the global variable section.
        The default is for Oakland, CA.

        Most of the code in this program is adapted from algorithms
        presented in "Practical Astronomy With Your Calculator" by
        Peter Duffet-Smith.

        The GST and ALT-AZIMUTH algorithms are from Sky and Telescope,
        June, 1984 by Roger W. Sinnott

        Adapted from:  sun.c
        Original author: Robert Bond - Beaverton, Oregon.

        Mods by:  D.N. Lynx Crowe - Oakland, California.
   ============================================================================ 
*/

#include "stdio.h"
#include "math.h"
#include "types.h"

#ifdef BSD
#include <sys/time.h>
#else
#include "time.h"
#endif

#include "stddefs.h"

#define OUTFILE         "daylight.dat"          /* output file name */

#ifndef PI
#define PI       3.141592654
#endif

#define EPOCH    1980
#define JDE      2444238.5      /* Julian date of EPOCH */

FILE    *outfp;                 /* output file pointer */

double dtor();
double adj360();
double adj24();
double julian_date();
double hms_to_dh();
double solar_lon();
double acos_deg();
double asin_deg();
double atan_q_deg();
double atan_deg();
double sin_deg();
double cos_deg();
double tan_deg();
double gmst();

int mo;
int day;
int yr;

int endmo;
int endday;
int endyr;

int popt = 0;                   /* position option flag */

int tz=8;                       /* Default time zone */
char *tzs  = "PST";             /* Default time zone string */
char *dtzs = "PDT";             /* Default daylight savings time string */

/* double lat = 45.5333;                /* Default latitude (Beaverton Ore.) */
/* double lon = 122.8333;               /* Default Longitude (Degrees west) */ 

double lat = 37.75;             /* Default latitude (Oakland, CA) */
double lon = 122.25;            /* Default longitude (degrees west) */

char mdays[] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };

char *
lite(hrz, mrz, hst, mst, nite, day, rise, set)
int hrz, mrz, hst, mst;
char nite, day, rise, set;
{
        register int i, h, m;
        static char niteday[49];

        if (mrz < 30)           /* quantize rise time */
                mrz = 0;
        else
                mrz = 30;

        if (mst < 30)           /* quantize set time */
                mst = 0;
        else
                mst = 30;

        for (i = 0; i < 48; i++) {

                h = i >> 1;
                m = (i & 1) ? 30 : 0;

                if ( (h < hrz) OR ((h EQ hrz) AND (m < mrz)) OR
                     (h > hst) OR ((h EQ hst) AND (m > mst)) )
                        niteday[i] = nite;
                else if ((h EQ hrz) AND (m EQ mrz))
                        niteday[i] = rise;
                else if ((h EQ hst) AND (m EQ mst))
                        niteday[i] = set;
                else
                        niteday[i] = day;
        }

        niteday[48] = '\0';
        return(&niteday[0]);
}

int
days(y)
register int y;
{
        if (y < 1970)
                y += 1900;

        if ( ((y % 4) EQ 0) AND ((y % 100) OR ((y % 400) EQ 0)) )
                y = 366;
        else
                y = 365;

        return(y);
}

/* 

*/

main(argc,argv)
int argc;
char *argv[];
{
        time_t sec_1970;
        int modays;
        struct tm *pt;

        sec_1970 = time((time_t *)0);  
        pt = localtime(&sec_1970);  

        yr  = pt->tm_year + 1900;
        mo  = pt->tm_mon + 1;
        day = pt->tm_mday;

        endyr = yr;
        endmo = mo;
        endday = day;

        if (pt->tm_isdst) {             /* convert tz to daylight savings time */

                tz--;
                tzs = dtzs;     
        }

        initopts(argc, argv);

        if ((FILE *)NULL EQ (outfp = fopen(OUTFILE, "w"))) {

                printf("ERROR:  Unable to open \"%s\" for output\n", OUTFILE);
                exit(1);
        }

        fprintf(outfp, 
                "Sunrise, Sunset (%s) and Hours of Daylight at Latitude %lf, Longitude %lf\n\n",
                tzs, lat, lon);

        while (TRUE) {

                dosun();

                if ((endyr EQ yr) AND (endmo EQ mo) AND (endday EQ day))
                        exit(0);

                modays = mdays[mo - 1];

                if ((days(yr) EQ 366) AND (mo EQ 2))
                        ++modays;

                if (++day > modays) {

                        day = 1;

                        if (++mo > 12) {

                                mo = 1;
                                ++yr;
                        }
                }
        }
}

dosun()
{
        double ed, jd;
        double alpha1, delta1, alpha2, delta2, st1r, st1s, st2r, st2s;
        double a1r, a1s, a2r, a2s, dt, x, y;
        double trise, tset, ar, as, delta, tri, da;
        double lambda1, lambda2;
        double m1;
        double hsm, ratio;
        int h, m, hrz, mrz, hdl, mdl;

        jd = julian_date(mo,day,yr);
        ed = jd - JDE;

        lambda1 = solar_lon(ed);
        lambda2 = solar_lon(ed + 1.0);

        lon_to_eq(lambda1, &alpha1, &delta1);
        lon_to_eq(lambda2, &alpha2, &delta2);

        rise_set(alpha1, delta1, &st1r, &st1s, &a1r, &a1s);
        rise_set(alpha2, delta2, &st2r, &st2s, &a2r, &a2s);

        m1 = adj24(gmst(jd - 0.5, 0.5 + tz / 24.0) - lon / 15); /* lst midnight */
        hsm = adj24(st1r - m1);
        ratio = hsm / 24.07;

        if (fabs(st2r - st1r) > 1.0)
                st2r += 24.0;

        trise = adj24((1.0 - ratio) * st1r + ratio * st2r);

        hsm = adj24(st1s - m1);
        ratio = hsm / 24.07;

        if (fabs(st2s - st1s) > 1.0)
                st2s += 24.0;

        tset = adj24((1.0 - ratio) * st1s + ratio * st2s);

        ar = a1r * 360.0 / (360.0 + a1r - a2r);
        as = a1s * 360.0 / (360.0 + a1s - a2s);

        delta = (delta1 + delta2) / 2.0;
        tri = acos_deg(sin_deg(lat)/cos_deg(delta));

        x = 0.835608;   /* correction for refraction, parallax, size of sun */
        y = asin_deg(sin_deg(x)/sin_deg(tri));
        da = asin_deg(tan_deg(x)/tan_deg(tri));
        dt = 240.0 * y / cos_deg(delta) / 3600;

        fprintf(outfp, "%04d-%02d-%02d  ", yr, mo, day);

        lst_to_hm(trise - dt, jd, &hrz, &mrz);
        fprintf(outfp, "%02d:%02d  ", hrz, mrz);

        lst_to_hm(tset + dt, jd, &h, &m);
        fprintf(outfp, "%02d:%02d  ", h, m);

        mdl = m - mrz;

        if (mdl < 0) {

                mdl += 60;

                h--;
        }

        hdl = h - hrz;

        fprintf(outfp, "%2d:%02d  ", hdl, mdl);

        fprintf(outfp, "%s\n", lite(hrz, mrz, h, m, '*', ' ', 'R', 'S'));
}

/* 

*/

double
dtor(deg)
double deg;
{
        return (deg * PI / 180.0);
}

double
rtod(deg)
double deg;
{
        return (deg * 180.0 / PI);
}


double
adj360(deg)
double deg;
{
        while (deg < 0.0) 
                deg += 360.0;

        while (deg > 360.0)
                deg -= 360.0;

        return(deg);
}

double
adj24(hrs)
double hrs;
{
        while (hrs < 0.0) 
                hrs += 24.0;

        while (hrs > 24.0)
                hrs -= 24.0;

        return(hrs);
}

/* 

*/

initopts(argc,argv)
int argc;
char *argv[];
{
        int ai;
        char *ca;
        char *str;

        while (--argc) {

                if ((*++argv)[0] == '-') {

                        ca = *argv;

                        for(ai = 1; ca[ai] != '\0'; ai++)
                                switch (ca[ai]) {

                                case 'a':
                                        str = *++argv;
                
                                        if (sscanf(str, "%lf", &lat) != 1)
                                                usage();

                                        argc--;
                                        break;

                                case 'o':
                                        str = *++argv;

                                        if (sscanf(str, "%lf", &lon) != 1)
                                                usage();

                                        argc--;
                                        break;

                                case 'z':
                                        str = *++argv;

                                        if (sscanf(str, "%d", &tz) != 1)
                                                usage();

                                        tzs = "   ";
                                        argc--;
                                        break;

                                case 'y':
                                        str = *++argv;

                                        if (sscanf(str, "%d", &yr) != 1)
                                                usage();

                                        mo = 1;
                                        day = 1;
                                        endmo = 12;
                                        endday = 31;
                                        argc--;
                                        break;

                                case 'd':
                                        str = *++argv;

                                        if (sscanf(str, "%d-%d-%d", &yr, &mo, &day) != 3)
                                                usage();

                                        argc--;
                                        break;

                                case 'e':
                                        str = *++argv;

                                        if (sscanf(str, "%d-%d-%d", &endyr, &endmo, &endday) != 3)
                                                usage();

                                        argc--;
                                        break;

                                default:
                                        usage();
                                }

                } else
                        usage();
        }
}

/* 

*/

usage()
{
        printf("Usage: sun [-p] [-t h:m:s] [-d m/d/y] [-a lat] [-o lon] [-z tz]\n");

        printf("-t hh:mm:ss\ttime (default is current system time)\n");
        printf("-d yyyy-mm-dd\tfirst date (default is current system date)\n");
        printf("-e yyyy-mm-dd\tlast date (default is current system date)\n");
        printf("-a lat\t\tdecimal latitude (default = 37.75)\n");
        printf("-o lon\t\tdecimal longitude (default = 122.25) \n");
        printf("-z tz\t\ttimezone (default = 8, PST)\n");
        printf("\nDefaults are for Oakland, CA\n");
        exit(1);
}

/* 

*/

double
julian_date(m, d, y)
{
        int a, b;
        double jd;

        if (m == 1 || m == 2) {

                --y;
                m += 12;
        }

        if (y < 1583) {

                printf("ERROR:  Can't handle dates before 1583 \n");
                exit(1);
        }

        a = y/100;
        b = 2 - a + a/4;
        b += (int)((double)y * 365.25);
        b += (int)(30.6001 * ((double)m + 1.0));
        jd = (double)d + (double)b + 1720994.5;
        return(jd);
}

double
hms_to_dh(h, m, s)
{
        double rv;

        rv = h + m / 60.0 + s / 3600.0;
        return rv;
}

/* 

*/

double
solar_lon(ed)
double ed;
{
        double n, m, e, ect, errt, v;

        n = 360.0 * ed / 365.2422;
        n = adj360(n);
        m = n + 278.83354 - 282.596403;
        m = adj360(m);
        m = dtor(m);
        e = m; ect = 0.016718;

        while ((errt = e - ect * sin(e) - m) > 0.0000001) 
                e = e - errt / (1 - ect * cos(e));

        v = 2 * atan(1.0168601 * tan(e/2));
        v = adj360(v * 180.0 / PI + 282.596403);
        return(v);
}

/* 

*/

double
acos_deg(x)
double x;
{
        return rtod(acos(x));
}

double
asin_deg(x)
double x;
{
        return rtod(asin(x));
}

double
atan_q_deg(y,x)
double y,x;
{
        double rv;

        if (y == 0)
                rv = 0;
        else if (x == 0)
                rv = y>0 ? 90.0 : -90.0;
        else
                rv = atan_deg(y/x);

        if (x<0)
                return rv+180.0;

        if (y<0)
                return rv+360.0;

        return(rv);
}

/* 

*/

double
atan_deg(x)
double x;
{
        return rtod(atan(x));
}

double
sin_deg(x)
double x;
{
        return sin(dtor(x));
}

double
cos_deg(x)
double x;
{
        return cos(dtor(x));
}

double
tan_deg(x)
double x;
{
        return tan(dtor(x));
}

/* 

*/

lon_to_eq(lambda, alpha, delta)
double lambda;
double *alpha;
double *delta;
{
        double tlam,epsilon;

        tlam = dtor(lambda);
        epsilon = dtor((double)23.441884);
        *alpha = atan_q_deg((sin(tlam))*cos(epsilon),cos(tlam)) / 15.0;
        *delta = asin_deg(sin(epsilon)*sin(tlam));
}

rise_set(alpha, delta, lstr, lsts, ar, as)
double alpha, delta, *lstr, *lsts, *ar, *as;
{
        double tar;
        double h;

        tar = sin_deg(delta)/cos_deg(lat);

        if (tar < -1.0 || tar > 1.0) {

                printf("ERROR:  The object is circumpolar \n");
                exit (1);
        }

        *ar = acos_deg(tar);
        *as = 360.0 - *ar;

        h = acos_deg(-tan_deg(lat) * tan_deg(delta)) / 15.0;
        *lstr = 24.0 + alpha - h;

        if (*lstr > 24.0)
                *lstr -= 24.0;

        *lsts = alpha + h;

        if (*lsts > 24.0)
                *lsts -= 24.0;
}

/* 

*/

lst_to_hm(lst, jd, h, m)
double lst, jd;
int *h, *m;
{
        double ed, gst, jzjd, t, r, b, t0, gmt;

        gst = lst + lon / 15.0;

        if (gst > 24.0)
                gst -= 24.0;

        jzjd = julian_date(1,0,yr);
        ed = jd-jzjd;
        t = (jzjd -2415020.0)/36525.0;
        r = 6.6460656+2400.05126*t+2.58E-05*t*t;
        b = 24-(r-24*(yr-1900));
        t0 = ed * 0.0657098 - b;

        if (t0 < 0.0)
                t0 += 24;

        gmt = gst-t0;

        if (gmt<0) 
                gmt += 24.0;

        gmt = gmt * 0.99727 - tz;;

        if (gmt < 0)
                gmt +=24.0;

        dh_to_hm(gmt, h, m);
}

/* 

*/

dh_to_hm(dh, h, m)
double dh;
int *h, *m;
{
        double tempsec;

        *h = dh;
        *m = (dh - *h) * 60;
        tempsec = (dh - *h) * 60 - *m;
        tempsec = tempsec * 60 + 0.5;

        if (tempsec > 30)
                (*m)++;

        if (*m == 60) {

                *m = 0;
                (*h)++;
        }
}

/* 

*/

eq_to_altaz(r, d, t, alt, az)
double r, d, t;
double *alt, *az;
{
        double p = 3.14159265;
        double r1 = p / 180.0;
        double b = lat * r1;
        double l = (360 - lon) * r1;
        double t5, s1, c1, c2, s2, a, h;

        r = r * 15.0 * r1;
        d = d * r1;
        t = t * 15.0 * r1;
        t5 = t - r + l;
        s1 = sin(b) * sin(d) + cos(b) * cos(d) * cos(t5);
        c1 = 1 - s1 * s1;

        if (c1 > 0) {

                c1 = sqrt(c1);
                h = atan(s1 / c1);

        } else {

                h = (s1 / fabs(s1)) * (p / 2.0);
        }

        c2 = cos(b) * sin(d) - sin(b) * cos(d) * cos(t5);
        s2 = -cos(d) * sin(t5);

        if (c2 == 0) 
                a = (s2/fabs(s2)) * (p/2);
        else {

        a = atan(s2/c2);

                if (c2 < 0)
                        a=a+p;
        }

        if (a<0)
                a=a+2*p;

        *alt = h / r1;
        *az = a / r1;
}

/* 

*/

double
gmst(j, f)
double j,f;
{
        double d, j0, t, t1, t2, s;

        d = j - 2451545.0;
        t = d / 36525.0;
        t1 = floor(t);
        j0 = t1 * 36525.0 + 2451545.0;
        t2 = (j - j0 + 0.5)/36525.0;
        s = 24110.54841 + 184.812866 * t1; 
        s += 8640184.812866 * t2;
        s += 0.093104 * t * t;
        s -= 0.0000062 * t * t * t;
        s /= 86400.0;
        s -= floor(s);
        s = 24 * (s + (f - 0.5) * 1.002737909);

        if (s < 0)
                s += 24.0;

        if (s > 24.0)
                s -= 24.0;

        return(s);
}