rgb@nsc-pdc.UUCP (Robert Bond) (05/10/85)
I finally got around to a manual page for the sunrise sunset program I posted
a while ago. Since I don't think the program got wide distribution, and
it's small, I threw it in too.
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#-----cut here-----cut here-----cut here-----cut here-----
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# sun.c
# sun.1
# makefile
# This archive created: Fri May 10 11:35:00 1985
# By: Robert Bond (NSC Portland, Orygun)
echo shar: extracting sun.c '(10571 characters)'
sed 's/^XX//' << \SHAR_EOF > sun.c
XX/* sun.c
XX*
XX* sun <options>
XX*
XX* options: -t hh:mm:ss time (default is current system time)
XX* -d mm/dd/yy date (default is current system date)
XX* -a lat decimal latitude (default = 45.5333)
XX* -o lon decimal longitude (default = 122.8333)
XX* -z tz timezone (default = 8, pst)
XX* -p show position of sun (azimuth)
XX*
XX* All output is to standard io.
XX*
XX* Compile with cc -O -o sun -lm
XX* Non 4.2 systems may have to change <sys/time.h> to <time.h> below.
XX*
XX* Note that the latitude, longitude, time zone correction and
XX* time zone string are all defaulted in the global variable section.
XX*
XX* Most of the code in this program is adapted from algorithms
XX* presented in "Practical Astronomy With Your Calculator" by
XX* Peter Duffet-Smith.
XX*
XX* The GST and ALT-AZIMUTH algorithms are from Sky and Telescope,
XX* June, 1984 by Roger W. Sinnott
XX*
XX* Author Robert Bond - Beaverton Oregon.
XX*
XX*/
XX
XX#include <stdio.h>
XX#include <math.h>
XX#include <sys/types.h>
XX#include <sys/time.h>
XX
XX#define PI 3.141592654
XX#define EPOCH 1980
XX#define JDE 2444238.5 /* Julian date of EPOCH */
XX
XXdouble dtor();
XXdouble adj360();
XXdouble adj24();
XXdouble julian_date();
XXdouble hms_to_dh();
XXdouble solar_lon();
XXdouble acos_deg();
XXdouble asin_deg();
XXdouble atan_q_deg();
XXdouble atan_deg();
XXdouble sin_deg();
XXdouble cos_deg();
XXdouble tan_deg();
XXdouble gmst();
XX
XXint th;
XXint tm;
XXint ts;
XXint mo;
XXint day;
XXint yr;
XXint tz=8; /* Default time zone */
XXchar *tzs = "(PST)"; /* Default time zone string */
XXchar *dtzs = "(PDT)"; /* Default daylight savings time string */
XXint popt = 0;
XX
XXdouble lat = 45.5333; /* Default latitude (Beaverton Ore.) */
XXdouble lon = 122.8333; /* Default Longitude (Degrees west) */
XX
XXmain(argc,argv)
XXint argc;
XXchar *argv[];
XX{
XX double ed, jd;
XX double alpha1, delta1, alpha2, delta2, st1r, st1s, st2r, st2s;
XX double a1r, a1s, a2r, a2s, dt, dh, x, y;
XX double trise, tset, ar, as, alpha, delta, tri, da;
XX double lambda1, lambda2;
XX double alt, az, gst, m1;
XX double hsm, ratio;
XX time_t sec_1970;
XX int h, m;
XX struct tm *pt;
XX
XX sec_1970 = time((time_t *)0);
XX pt = localtime(&sec_1970);
XX
XX th = pt->tm_hour;
XX tm = pt->tm_min;
XX ts = pt->tm_sec;
XX yr = pt->tm_year + 1900;
XX mo = pt->tm_mon + 1;
XX day = pt->tm_mday;
XX if (pt->tm_isdst) { /* convert tz to daylight savings time */
XX tz--;
XX tzs = dtzs;
XX }
XX
XX initopts(argc,argv);
XX
XX jd = julian_date(mo,day,yr);
XX ed = jd - JDE;
XX
XX lambda1 = solar_lon(ed);
XX lambda2 = solar_lon(ed + 1.0);
XX
XX lon_to_eq(lambda1, &alpha1, &delta1);
XX lon_to_eq(lambda2, &alpha2, &delta2);
XX
XX rise_set(alpha1, delta1, &st1r, &st1s, &a1r, &a1s);
XX rise_set(alpha2, delta2, &st2r, &st2s, &a2r, &a2s);
XX
XX m1 = adj24(gmst(jd - 0.5, 0.5 + tz / 24.0) - lon / 15); /* lst midnight */
XX hsm = adj24(st1r - m1);
XX ratio = hsm / 24.07;
XX
XX if (fabs(st2r - st1r) > 1.0) {
XX st2r += 24.0;
XX }
XX
XX trise = adj24((1.0 - ratio) * st1r + ratio * st2r);
XX
XX hsm = adj24(st1s - m1);
XX ratio = hsm / 24.07;
XX
XX if (fabs(st2s - st1s) > 1.0) {
XX st2s += 24.0;
XX }
XX
XX tset = adj24((1.0 - ratio) * st1s + ratio * st2s);
XX
XX ar = a1r * 360.0 / (360.0 + a1r - a2r);
XX as = a1s * 360.0 / (360.0 + a1s - a2s);
XX
XX delta = (delta1 + delta2) / 2.0;
XX tri = acos_deg(sin_deg(lat)/cos_deg(delta));
XX
XX x = 0.835608; /* correction for refraction, parallax, size of sun */
XX y = asin_deg(sin_deg(x)/sin_deg(tri));
XX da = asin_deg(tan_deg(x)/tan_deg(tri));
XX dt = 240.0 * y / cos_deg(delta) / 3600;
XX
XX lst_to_hm(trise - dt, jd, &h, &m);
XX printf("Sunrise: %2d:%02d ", h, m);
XX
XX if (popt) {
XX dh_to_hm(ar - da, &h, &m);
XX printf("Azimuth: %2d deg %02d min \n", h, m);
XX }
XX
XX lst_to_hm(tset + dt, jd, &h, &m);
XX printf("Sunset: %2d:%02d ", h, m);
XX
XX if (popt) {
XX dh_to_hm(as + da, &h, &m);
XX printf("Azimuth: %2d deg %02d min \n", h, m);
XX } else
XX printf("%s \n",tzs);
XX
XX if (popt) {
XX
XX if (alpha1 < alpha2)
XX alpha = (alpha1 + alpha2) / 2.0;
XX else
XX alpha = (alpha1 + 24.0 + alpha2) / 2.0;
XX
XX if (alpha > 24.0)
XX alpha -= 24.0;
XX
XX dh = (hms_to_dh(th, tm, ts) + tz) / 24.0;
XX if (dh > 0.5) {
XX dh -= 0.5;
XX jd += 0.5;
XX } else {
XX dh += 0.5;
XX jd -= 0.5;
XX }
XX
XX gst = gmst(jd, dh);
XX
XX eq_to_altaz(alpha, delta, gst, &alt, &az);
XX
XX printf("The sun is at ");
XX dh_to_hm(alt, &h, &m);
XX printf(" altitude: %2d deg %02d min", h, m);
XX dh_to_hm(az, &h, &m);
XX printf(" azimuth: %2d deg %02d min. \n", h, m);
XX }
XX}
XX
XXdouble
XXdtor(deg)
XXdouble deg;
XX{
XXreturn (deg * PI / 180.0);
XX}
XX
XXdouble
XXrtod(deg)
XXdouble deg;
XX{
XX return (deg * 180.0 / PI);
XX}
XX
XX
XXdouble
XXadj360(deg)
XXdouble deg;
XX{
XX while (deg < 0.0)
XX deg += 360.0;
XX while (deg > 360.0)
XX deg -= 360.0;
XX return(deg);
XX}
XX
XXdouble
XXadj24(hrs)
XXdouble hrs;
XX{
XX while (hrs < 0.0)
XX hrs += 24.0;
XX while (hrs > 24.0)
XX hrs -= 24.0;
XX return(hrs);
XX}
XX
XXinitopts(argc,argv)
XXint argc;
XXchar *argv[];
XX{
XX int ai;
XX char *ca;
XX char *str;
XX
XX while (--argc) {
XX if ((*++argv)[0] == '-') {
XX ca = *argv;
XX for(ai = 1; ca[ai] != '\0'; ai++)
XX switch (ca[ai]) {
XX case 'p':
XX popt++;
XX break;
XX case 'a':
XX str = *++argv;
XX if (sscanf(str, "%lf", &lat) != 1)
XX usage();
XX argc--;
XX break;
XX case 'o':
XX str = *++argv;
XX if (sscanf(str, "%lf", &lon) != 1)
XX usage();
XX argc--;
XX break;
XX case 'z':
XX str = *++argv;
XX if (sscanf(str, "%d", &tz) != 1)
XX usage();
XX tzs = " ";
XX argc--;
XX break;
XX case 't':
XX str = *++argv;
XX if (sscanf(str, "%d:%d:%d", &th, &tm, &ts) != 3)
XX usage();
XX argc--;
XX break;
XX case 'd':
XX str = *++argv;
XX if (sscanf(str, "%d/%d/%d", &mo, &day, &yr) != 3)
XX usage();
XX argc--;
XX break;
XX default: usage();
XX }
XX } else usage();
XX }
XX}
XX
XXusage()
XX{
XX printf("Usage: sun [-p] [-t h:m:s] [-d m/d/y] [-a lat] [-o lon] [-z tz]\n");
XX exit(1);
XX}
XX
XXdouble
XXjulian_date(m, d, y)
XX{
XX int a, b;
XX double jd;
XX
XX if (m == 1 || m == 2) {
XX --y;
XX m += 12;
XX }
XX if (y < 1583) {
XX printf("Can't handle dates before 1583 \n");
XX exit(1);
XX }
XX a = y/100;
XX b = 2 - a + a/4;
XX b += (int)((double)y * 365.25);
XX b += (int)(30.6001 * ((double)m + 1.0));
XX jd = (double)d + (double)b + 1720994.5;
XX return(jd);
XX}
XX
XXdouble
XXhms_to_dh(h, m, s)
XX{
XX double rv;
XX
XX rv = h + m / 60.0 + s / 3600.0;
XX return rv;
XX}
XX
XXdouble
XXsolar_lon(ed)
XXdouble ed;
XX{
XX double n, m, e, ect, errt, v;
XX
XX n = 360.0 * ed / 365.2422;
XX n = adj360(n);
XX m = n + 278.83354 - 282.596403;
XX m = adj360(m);
XX m = dtor(m);
XX e = m; ect = 0.016718;
XX while ((errt = e - ect * sin(e) - m) > 0.0000001)
XX e = e - errt / (1 - ect * cos(e));
XX v = 2 * atan(1.0168601 * tan(e/2));
XX v = adj360(v * 180.0 / PI + 282.596403);
XX return(v);
XX}
XX
XXdouble
XXacos_deg(x)
XXdouble x;
XX{
XX return rtod(acos(x));
XX}
XX
XXdouble
XXasin_deg(x)
XXdouble x;
XX{
XX return rtod(asin(x));
XX}
XX
XXdouble
XXatan_q_deg(y,x)
XXdouble y,x;
XX{
XX double rv;
XX
XX if (y == 0)
XX rv = 0;
XX else if (x == 0)
XX rv = y>0 ? 90.0 : -90.0;
XX else rv = atan_deg(y/x);
XX
XX if (x<0) return rv+180.0;
XX if (y<0) return rv+360.0;
XX return(rv);
XX}
XX
XXdouble
XXatan_deg(x)
XXdouble x;
XX{
XX return rtod(atan(x));
XX}
XX
XXdouble
XXsin_deg(x)
XXdouble x;
XX{
XX return sin(dtor(x));
XX}
XX
XXdouble
XXcos_deg(x)
XXdouble x;
XX{
XX return cos(dtor(x));
XX}
XX
XXdouble
XXtan_deg(x)
XXdouble x;
XX{
XX return tan(dtor(x));
XX}
XX
XXlon_to_eq(lambda, alpha, delta)
XXdouble lambda;
XXdouble *alpha;
XXdouble *delta;
XX{
XX double tlam,epsilon;
XX
XX tlam = dtor(lambda);
XX epsilon = dtor((double)23.441884);
XX *alpha = atan_q_deg((sin(tlam))*cos(epsilon),cos(tlam)) / 15.0;
XX *delta = asin_deg(sin(epsilon)*sin(tlam));
XX}
XX
XXrise_set(alpha, delta, lstr, lsts, ar, as)
XXdouble alpha, delta, *lstr, *lsts, *ar, *as;
XX{
XX double tar;
XX double h;
XX
XX tar = sin_deg(delta)/cos_deg(lat);
XX if (tar < -1.0 || tar > 1.0) {
XX printf("The object is circumpolar \n");
XX exit (1);
XX }
XX *ar = acos_deg(tar);
XX *as = 360.0 - *ar;
XX
XX h = acos_deg(-tan_deg(lat) * tan_deg(delta)) / 15.0;
XX *lstr = 24.0 + alpha - h;
XX if (*lstr > 24.0)
XX *lstr -= 24.0;
XX *lsts = alpha + h;
XX if (*lsts > 24.0)
XX *lsts -= 24.0;
XX}
XX
XXlst_to_hm(lst, jd, h, m)
XXdouble lst, jd;
XXint *h, *m;
XX{
XX double ed, gst, jzjd, t, r, b, t0, gmt;
XX
XX gst = lst + lon / 15.0;
XX if (gst > 24.0)
XX gst -= 24.0;
XX jzjd = julian_date(1,0,yr);
XX ed = jd-jzjd;
XX t = (jzjd -2415020.0)/36525.0;
XX r = 6.6460656+2400.05126*t+2.58E-05*t*t;
XX b = 24-(r-24*(yr-1900));
XX t0 = ed * 0.0657098 - b;
XX if (t0 < 0.0)
XX t0 += 24;
XX gmt = gst-t0;
XX if (gmt<0)
XX gmt += 24.0;
XX gmt = gmt * 0.99727 - tz;;
XX if (gmt < 0)
XX gmt +=24.0;
XX dh_to_hm(gmt, h, m);
XX}
XX
XXdh_to_hm(dh, h, m)
XXdouble dh;
XXint *h, *m;
XX{
XX double tempsec;
XX
XX *h = dh;
XX *m = (dh - *h) * 60;
XX tempsec = (dh - *h) * 60 - *m;
XX tempsec = tempsec * 60 + 0.5;
XX if (tempsec > 30)
XX (*m)++;
XX if (*m == 60) {
XX *m = 0;
XX (*h)++;
XX }
XX}
XX
XXeq_to_altaz(r, d, t, alt, az)
XXdouble r, d, t;
XXdouble *alt, *az;
XX{
XX double p = 3.14159265;
XX double r1 = p / 180.0;
XX double b = lat * r1;
XX double l = (360 - lon) * r1;
XX double t5, s1, c1, c2, s2, a, h;
XX
XX r = r * 15.0 * r1;
XX d = d * r1;
XX t = t * 15.0 * r1;
XX t5 = t - r + l;
XX s1 = sin(b) * sin(d) + cos(b) * cos(d) * cos(t5);
XX c1 = 1 - s1 * s1;
XX if (c1 > 0) {
XX c1 = sqrt(c1);
XX h = atan(s1 / c1);
XX } else {
XX h = (s1 / fabs(s1)) * (p / 2.0);
XX }
XX c2 = cos(b) * sin(d) - sin(b) * cos(d) * cos(t5);
XX s2 = -cos(d) * sin(t5);
XX if (c2 == 0)
XX a = (s2/fabs(s2)) * (p/2);
XX else {
XX a = atan(s2/c2);
XX if (c2 < 0)
XX a=a+p;
XX }
XX if (a<0)
XX a=a+2*p;
XX *alt = h / r1;
XX *az = a / r1;
XX}
XX
XXdouble
XXgmst(j, f)
XXdouble j,f;
XX{
XX double d, j0, t, t1, t2, s;
XX
XX d = j - 2451545.0;
XX t = d / 36525.0;
XX t1 = floor(t);
XX j0 = t1 * 36525.0 + 2451545.0;
XX t2 = (j - j0 + 0.5)/36525.0;
XX s = 24110.54841 + 184.812866 * t1;
XX s += 8640184.812866 * t2;
XX s += 0.093104 * t * t;
XX s -= 0.0000062 * t * t * t;
XX s /= 86400.0;
XX s -= floor(s);
XX s = 24 * (s + (f - 0.5) * 1.002737909);
XX if (s < 0)
XX s += 24.0;
XX if (s > 24.0)
XX s -= 24.0;
XX return(s);
XX}
XX
XX
SHAR_EOF
if test 10571 -ne "`wc -c sun.c`"
then
echo shar: error transmitting sun.c '(should have been 10571 characters)'
fi
echo shar: extracting sun.1 '(1625 characters)'
sed 's/^XX//' << \SHAR_EOF > sun.1
XX.TH SUN 1 "9 May 1985"
XX.UC 4
XX.SH NAME
XXsun \- calculate sunrise and sunset times
XX.SH SYNOPSIS
XX.B sun
XX[
XX.B \-p
XX]
XX[
XX.B -t \fIh:m:s\fP
XX]
XX[
XX.B -d \fIm/d/y\fP
XX]
XX[
XX.B -a \fIlat\fP
XX]
XX[
XX.B -o \fIlon\fP
XX]
XX [
XX.B -z \fItz\fP
XX]
XX.LP
XX.SH DESCRIPTION
XX.I Sun
XXCalculates the times of sunrise and sunset for the current system date.
XXOptionally,
XX.I sun
XXcalculates the current position of the sun and its position at sunrise and
XXsunset for other dates, times and places.
XX.PP
XXThe command line options are:
XX.TP
XX.B \-p
XXShow current altitude and azimuth for the sun and its azimuth
XXat sunrise and sunset.
XX.TP
XX.B -t \fIh:m:s\fP
XXUse \fIh:m:s\fP as local hours (24 hour clock), minutes and seconds as
XXthe time instead of the current system time.
XX.TP
XX.B -d \fIm/d/y\fP
XXUse \fIm/d/y\fP as the month, day and year instead of the current system date.
XX.TP
XX.B -a \fIlat\fP
XXUse \fIlat\fP as the latitude of the observer, instead of the system latitude.
XXThe latitude is specified in decimal degrees north. For example, use
XX45.5333 for Beaverton, Oregon.
XX.TP
XX.B -o \fIlon\fP
XXUse \fIlon\fP as the longitude of the observer, instead of the system longitude.
XXThe longitude is specified in decimal degrees west. For example, use
XX122.8333 for Beaverton, Oregon.
XX.TP
XX.B -z \fItz\fP
XXUse \fItz\fP as the time zone of the observer, instead of the current system
XXtimezone. Timezone is specified as a positive number for locations west of
XXGreenwich. For example, use 8 for Pacific Standard Time.
XX.SH AUTHOR
XXRobert Bond
XX.SH BUGS
XXShould figure out timezone and guess at longitude from the system time.
XXTimezone, latitude and longitude are compiled in.
SHAR_EOF
if test 1625 -ne "`wc -c sun.1`"
then
echo shar: error transmitting sun.1 '(should have been 1625 characters)'
fi
echo shar: extracting makefile '(82 characters)'
sed 's/^XX//' << \SHAR_EOF > makefile
XX
XXsun: sun.c
XX cc -O -o sun sun.c -lm
XX
XXmanual: sun.1
XX nroff -man sun.1 >sun.doc
SHAR_EOF
if test 82 -ne "`wc -c makefile`"
then
echo shar: error transmitting makefile '(should have been 82 characters)'
fi
# End of shell archive
exit 0
--
Robert Bond nsc!nsc-pdc!rgb
National Semiconductor tektronix!reed!nsc-pdc!rgb