| >> Well, after working almost a month I've got a 'C' program that seems to | >> get | >> within 6 minutes of the U.S. Navy's published times. I actually learned a | >> lot about vector operations. For now I'm going to put this thing to rest, | >> but here's the code... | >
| > Could you post it again without losing all the line breaks | > in the second half? | >
| > Thanks | > -- | > Andrew | | In the almost three weeks since I posted this 'C' program I cleaned it up a | little but kept the same logic. I moved some code into subroutines since it | was used several times from different places. Also, my first posting didn't | include the source for the cross-product and dot-product of vectors since | these were pulled from a VCTRSUPPORT.h file I had previously developed. I am | still at a loss to figure out why I've got about 6 minutes of error (seems | to vary from about 5 to 7 depending on the time of year) but am not willing | "fudge" any of my geometrically accurate analysis to eliminate it.
I have a program that I found years ago that computes sunset/sunrise. It has a similar error which I've always wondered about. Just in case somebody wants more/alternate code to look at, here it is:
/* sun.c
*
*
- options: -t hh:mm:ss time (default is current system time)
* -d mm/dd/yy date (default is current system date)
* -z tz timezone (default = 8, pst)
* -p show position of sun (azimuth)
*
- All output is to standard io.
*
* Compile with cc -O -o sun -lm
* Non 4.2 systems may have to change to below.
*
* Note that the latitude, longitude, time zone correction and
* time zone string are all defaulted in the global variable section.
*
* 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
*
* Author Robert Bond - Beaverton Oregon.
*
*/
#include #include #include #include
#define PI 3.141592654 #define EPOCH 1980 #define JDE 2444238.5 /* Julian date of EPOCH */
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 th; int tm; int ts; int mo; int day; int yr; int tz=5; /* Default time zone */ char *tzs = "(EST)"; /* Default time zone string */ char *dtzs = "(EDT)"; /* Default daylight savings time string */ int popt = 0, qopt = 0;
double lat = 42.59985; /* Default latitude */ double lon = 70.63733; /* Default Longitude (Degrees west) */
main(argc,argv) int argc; char *argv[]; { double ed, jd; double alpha1, delta1, alpha2, delta2, st1r, st1s, st2r, st2s; double a1r, a1s, a2r, a2s, dt, dh, x, y; double trise, tset, ar, as, alpha, delta, tri, da; double lambda1, lambda2; double alt, az, gst, m1; double hsm, ratio; time_t sec_1970; int h, m, mrise, mset; struct tm *pt;
sec_1970 = time((time_t *)0); pt = localtime(&sec_1970);
th = pt->tm_hour; tm = pt->tm_min; ts = pt->tm_sec; yr = pt->tm_year + 1900; mo = pt->tm_mon + 1; day = pt->tm_mday; if (pt->tm_isdst) { /* convert tz to daylight savings time */ tz--; tzs = dtzs; }
initopts(argc,argv);
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;
lst_to_hm(trise - dt, jd, &h, &m); mrise = m + 60 * h; if(!qopt) printf("Sunrise: %2d:%02d ", h, m);
if (popt) { dh_to_hm(ar - da, &h, &m); printf("Azimuth: %2d deg %02d min \\n", h, m); }
lst_to_hm(tset + dt, jd, &h, &m); mset = m + 60 * h; if(!qopt) printf("Sunset: %2d:%02d ", h, m);
if (popt) { dh_to_hm(as + da, &h, &m); printf("Azimuth: %2d deg %02d min \\n", h, m); } else if(!qopt) printf("%s \\n",tzs);
if (popt) {
if (alpha1 < alpha2) alpha = (alpha1 + alpha2) / 2.0; else alpha = (alpha1 + 24.0 + alpha2) / 2.0; if (alpha > 24.0) alpha -= 24.0;
dh = (hms_to_dh(th, tm, ts) + tz) / 24.0; if (dh > 0.5) { dh -= 0.5; jd += 0.5; } else { dh += 0.5; jd -= 0.5; }
gst = gmst(jd, dh);
eq_to_altaz(alpha, delta, gst, &alt, &az);
printf("The sun is at "); dh_to_hm(alt, &h, &m); printf(" altitude: %2d deg %02d min", h, m); dh_to_hm(az, &h, &m); printf(" azimuth: %2d deg %02d min. \\n", h, m); } m = pt->tm_min + 60 * pt->tm_hour; if(m > mrise && m < mset) exit(1); exit(0); }
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 'q': qopt++; break; case 'p': popt++; break; 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 't': str = *++argv; if (sscanf(str, "%d:%d:%d", &th, &tm, &ts) != 3) usage(); argc--; break; case 'd': str = *++argv; if (sscanf(str, "%d/%d/%d", &mo, &day, &yr) != 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"); 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("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 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 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 24.0) s -= 24.0; return(s); }