source: buchla-68k/orig/DATE/SUN.C@ 66b48e7

Last change on this file since 66b48e7 was 3ae31e9, checked in by Thomas Lopatic <thomas@…>, 7 years ago

Imported original source code.
  • Property mode set to 100755
File size: 11.2 KB
Line 
1/*
2 ============================================================================
3 sun.c -- calculate sunrise, sunset, and position of the sun
4 Version 5 -- 1989-01-12 -- Robert Bond / D.N. Lynx Crowe
5
6 sun <options>
7
8 options: -t hh:mm:ss time (default is current system time)
9 -d mm/dd/yy date (default is current system date)
10 -a lat decimal latitude (default = 45.5333)
11 -o lon decimal longitude (default = 122.8333)
12 -z tz timezone (default = 8, pst)
13 -p show position of sun (azimuth)
14
15 All output is to standard io.
16
17 Compile with cc -O -o sun -lm
18
19 Define BSD for BSD4.2 systems. System V assumed otherwise.
20 This selects between <sys/time.h> and <time.h> below.
21
22 Note that the latitude, longitude, time zone correction and
23 time zone string are all defaulted in the global variable section.
24 The default is for Oakland, CA.
25
26 Most of the code in this program is adapted from algorithms
27 presented in "Practical Astronomy With Your Calculator" by
28 Peter Duffet-Smith.
29
30 The GST and ALT-AZIMUTH algorithms are from Sky and Telescope,
31 June, 1984 by Roger W. Sinnott
32
33 Original author: Robert Bond - Beaverton, Oregon.
34 Mods by: D.N. Lynx Crowe - Oakland, California.
35 ============================================================================
36*/
37
38#include "stdio.h"
39#include "math.h"
40#include "types.h"
41
42#ifdef BSD
43#include <sys/time.h>
44#else
45#include "time.h"
46#endif
47
48#ifndef PI
49#define PI 3.141592654
50#endif
51
52#define EPOCH 1980
53#define JDE 2444238.5 /* Julian date of EPOCH */
54
55double dtor();
56double adj360();
57double adj24();
58double julian_date();
59double hms_to_dh();
60double solar_lon();
61double acos_deg();
62double asin_deg();
63double atan_q_deg();
64double atan_deg();
65double sin_deg();
66double cos_deg();
67double tan_deg();
68double gmst();
69
70int th;
71int tm;
72int ts;
73int mo;
74int day;
75int yr;
76int popt = 0;
77
78int tz=8; /* Default time zone */
79char *tzs = "(PST)"; /* Default time zone string */
80char *dtzs = "(PDT)"; /* Default daylight savings time string */
81
82/* double lat = 45.5333; /* Default latitude (Beaverton Ore.) */
83/* double lon = 122.8333; /* Default Longitude (Degrees west) */
84
85double lat = 37.75; /* Default latitude (Oakland, CA) */
86double lon = 122.25; /* Default longitude (degrees west) */
87
88/*
89
90*/
91
92main(argc,argv)
93int argc;
94char *argv[];
95{
96 double ed, jd;
97 double alpha1, delta1, alpha2, delta2, st1r, st1s, st2r, st2s;
98 double a1r, a1s, a2r, a2s, dt, dh, x, y;
99 double trise, tset, ar, as, alpha, delta, tri, da;
100 double lambda1, lambda2;
101 double alt, az, gst, m1;
102 double hsm, ratio;
103 time_t sec_1970;
104 int h, m;
105 struct tm *pt;
106
107 sec_1970 = time((time_t *)0);
108 pt = localtime(&sec_1970);
109
110 th = pt->tm_hour;
111 tm = pt->tm_min;
112 ts = pt->tm_sec;
113 yr = pt->tm_year + 1900;
114 mo = pt->tm_mon + 1;
115 day = pt->tm_mday;
116
117 if (pt->tm_isdst) { /* convert tz to daylight savings time */
118
119 tz--;
120 tzs = dtzs;
121 }
122
123 initopts(argc,argv);
124
125 printf("On %04d-%02d-%02d at %02d:%02d:%02d\n",
126 yr, mo, day, th, tm, ts);
127
128 jd = julian_date(mo,day,yr);
129 ed = jd - JDE;
130
131 lambda1 = solar_lon(ed);
132 lambda2 = solar_lon(ed + 1.0);
133
134 lon_to_eq(lambda1, &alpha1, &delta1);
135 lon_to_eq(lambda2, &alpha2, &delta2);
136
137 rise_set(alpha1, delta1, &st1r, &st1s, &a1r, &a1s);
138 rise_set(alpha2, delta2, &st2r, &st2s, &a2r, &a2s);
139
140 m1 = adj24(gmst(jd - 0.5, 0.5 + tz / 24.0) - lon / 15); /* lst midnight */
141 hsm = adj24(st1r - m1);
142 ratio = hsm / 24.07;
143
144 if (fabs(st2r - st1r) > 1.0) {
145
146 st2r += 24.0;
147 }
148
149 trise = adj24((1.0 - ratio) * st1r + ratio * st2r);
150
151 hsm = adj24(st1s - m1);
152 ratio = hsm / 24.07;
153
154 if (fabs(st2s - st1s) > 1.0) {
155
156 st2s += 24.0;
157 }
158
159 tset = adj24((1.0 - ratio) * st1s + ratio * st2s);
160
161 ar = a1r * 360.0 / (360.0 + a1r - a2r);
162 as = a1s * 360.0 / (360.0 + a1s - a2s);
163
164 delta = (delta1 + delta2) / 2.0;
165 tri = acos_deg(sin_deg(lat)/cos_deg(delta));
166
167 x = 0.835608; /* correction for refraction, parallax, size of sun */
168 y = asin_deg(sin_deg(x)/sin_deg(tri));
169 da = asin_deg(tan_deg(x)/tan_deg(tri));
170 dt = 240.0 * y / cos_deg(delta) / 3600;
171
172 lst_to_hm(trise - dt, jd, &h, &m);
173 printf("Sunrise: %2d:%02d ", h, m);
174
175 if (popt) {
176
177 dh_to_hm(ar - da, &h, &m);
178 printf("Azimuth: %2d deg %02d min \n", h, m);
179 }
180
181 lst_to_hm(tset + dt, jd, &h, &m);
182 printf("Sunset: %2d:%02d ", h, m);
183
184 if (popt) {
185
186 dh_to_hm(as + da, &h, &m);
187 printf("Azimuth: %2d deg %02d min \n", h, m);
188
189 } else
190 printf("%s \n",tzs);
191
192 if (popt) {
193
194 if (alpha1 < alpha2)
195 alpha = (alpha1 + alpha2) / 2.0;
196 else
197 alpha = (alpha1 + 24.0 + alpha2) / 2.0;
198
199 if (alpha > 24.0)
200 alpha -= 24.0;
201
202 dh = (hms_to_dh(th, tm, ts) + tz) / 24.0;
203
204 if (dh > 0.5) {
205
206 dh -= 0.5;
207 jd += 0.5;
208
209 } else {
210
211 dh += 0.5;
212 jd -= 0.5;
213 }
214
215 gst = gmst(jd, dh);
216
217 eq_to_altaz(alpha, delta, gst, &alt, &az);
218
219 printf("The sun is at ");
220 dh_to_hm(alt, &h, &m);
221 printf("Altitude: %2d deg %02d min", h, m);
222 dh_to_hm(az, &h, &m);
223 printf(" Azimuth: %2d deg %02d min. \n", h, m);
224 }
225}
226
227/*
228
229*/
230
231double
232dtor(deg)
233double deg;
234{
235 return (deg * PI / 180.0);
236}
237
238double
239rtod(deg)
240double deg;
241{
242 return (deg * 180.0 / PI);
243}
244
245
246double
247adj360(deg)
248double deg;
249{
250 while (deg < 0.0)
251 deg += 360.0;
252
253 while (deg > 360.0)
254 deg -= 360.0;
255
256 return(deg);
257}
258
259double
260adj24(hrs)
261double hrs;
262{
263 while (hrs < 0.0)
264 hrs += 24.0;
265
266 while (hrs > 24.0)
267 hrs -= 24.0;
268
269 return(hrs);
270}
271
272/*
273
274*/
275
276initopts(argc,argv)
277int argc;
278char *argv[];
279{
280 int ai;
281 char *ca;
282 char *str;
283
284 while (--argc) {
285
286 if ((*++argv)[0] == '-') {
287
288 ca = *argv;
289
290 for(ai = 1; ca[ai] != '\0'; ai++)
291 switch (ca[ai]) {
292
293 case 'p':
294 popt++;
295 break;
296
297 case 'a':
298 str = *++argv;
299
300 if (sscanf(str, "%lf", &lat) != 1)
301 usage();
302
303 argc--;
304 break;
305
306 case 'o':
307 str = *++argv;
308
309 if (sscanf(str, "%lf", &lon) != 1)
310 usage();
311
312 argc--;
313 break;
314
315 case 'z':
316 str = *++argv;
317
318 if (sscanf(str, "%d", &tz) != 1)
319 usage();
320
321 tzs = " ";
322 argc--;
323 break;
324
325 case 't':
326 str = *++argv;
327
328 if (sscanf(str, "%d:%d:%d", &th, &tm, &ts) != 3)
329 usage();
330
331 argc--;
332 break;
333
334 case 'd':
335 str = *++argv;
336
337 if (sscanf(str, "%d/%d/%d", &mo, &day, &yr) != 3)
338 usage();
339
340 argc--;
341 break;
342
343 default:
344 usage();
345 }
346
347 } else
348 usage();
349 }
350}
351
352/*
353
354*/
355
356usage()
357{
358 printf("Usage: sun [-p] [-t h:m:s] [-d m/d/y] [-a lat] [-o lon] [-z tz]\n");
359
360 printf("-p\t\tshow position of sun (azimuth)\n");
361 printf("-t hh:mm:ss\ttime (default is current system time)\n");
362 printf("-d mm/dd/yy\tdate (default is current system date)\n");
363 printf("-a lat\t\tdecimal latitude (default = 37.75)\n");
364 printf("-o lon\t\tdecimal longitude (default = 122.25) \n");
365 printf("-z tz\t\ttimezone (default = 8, pst)\n");
366 printf("\nDefaults are for Oakland, CA\n");
367 exit(1);
368}
369
370/*
371
372*/
373
374double
375julian_date(m, d, y)
376{
377 int a, b;
378 double jd;
379
380 if (m == 1 || m == 2) {
381
382 --y;
383 m += 12;
384 }
385
386 if (y < 1583) {
387
388 printf("Can't handle dates before 1583 \n");
389 exit(1);
390 }
391
392 a = y/100;
393 b = 2 - a + a/4;
394 b += (int)((double)y * 365.25);
395 b += (int)(30.6001 * ((double)m + 1.0));
396 jd = (double)d + (double)b + 1720994.5;
397 return(jd);
398}
399
400double
401hms_to_dh(h, m, s)
402{
403 double rv;
404
405 rv = h + m / 60.0 + s / 3600.0;
406 return rv;
407}
408
409/*
410
411*/
412
413double
414solar_lon(ed)
415double ed;
416{
417 double n, m, e, ect, errt, v;
418
419 n = 360.0 * ed / 365.2422;
420 n = adj360(n);
421 m = n + 278.83354 - 282.596403;
422 m = adj360(m);
423 m = dtor(m);
424 e = m; ect = 0.016718;
425
426 while ((errt = e - ect * sin(e) - m) > 0.0000001)
427 e = e - errt / (1 - ect * cos(e));
428
429 v = 2 * atan(1.0168601 * tan(e/2));
430 v = adj360(v * 180.0 / PI + 282.596403);
431 return(v);
432}
433
434/*
435
436*/
437
438double
439acos_deg(x)
440double x;
441{
442 return rtod(acos(x));
443}
444
445double
446asin_deg(x)
447double x;
448{
449 return rtod(asin(x));
450}
451
452double
453atan_q_deg(y,x)
454double y,x;
455{
456 double rv;
457
458 if (y == 0)
459 rv = 0;
460 else if (x == 0)
461 rv = y>0 ? 90.0 : -90.0;
462 else
463 rv = atan_deg(y/x);
464
465 if (x<0)
466 return rv+180.0;
467
468 if (y<0)
469 return rv+360.0;
470
471 return(rv);
472}
473
474/*
475
476*/
477
478double
479atan_deg(x)
480double x;
481{
482 return rtod(atan(x));
483}
484
485double
486sin_deg(x)
487double x;
488{
489 return sin(dtor(x));
490}
491
492double
493cos_deg(x)
494double x;
495{
496 return cos(dtor(x));
497}
498
499double
500tan_deg(x)
501double x;
502{
503 return tan(dtor(x));
504}
505
506/*
507
508*/
509
510lon_to_eq(lambda, alpha, delta)
511double lambda;
512double *alpha;
513double *delta;
514{
515 double tlam,epsilon;
516
517 tlam = dtor(lambda);
518 epsilon = dtor((double)23.441884);
519 *alpha = atan_q_deg((sin(tlam))*cos(epsilon),cos(tlam)) / 15.0;
520 *delta = asin_deg(sin(epsilon)*sin(tlam));
521}
522
523rise_set(alpha, delta, lstr, lsts, ar, as)
524double alpha, delta, *lstr, *lsts, *ar, *as;
525{
526 double tar;
527 double h;
528
529 tar = sin_deg(delta)/cos_deg(lat);
530
531 if (tar < -1.0 || tar > 1.0) {
532
533 printf("The object is circumpolar \n");
534 exit (1);
535 }
536
537 *ar = acos_deg(tar);
538 *as = 360.0 - *ar;
539
540 h = acos_deg(-tan_deg(lat) * tan_deg(delta)) / 15.0;
541 *lstr = 24.0 + alpha - h;
542
543 if (*lstr > 24.0)
544 *lstr -= 24.0;
545
546 *lsts = alpha + h;
547
548 if (*lsts > 24.0)
549 *lsts -= 24.0;
550}
551
552/*
553
554*/
555
556lst_to_hm(lst, jd, h, m)
557double lst, jd;
558int *h, *m;
559{
560 double ed, gst, jzjd, t, r, b, t0, gmt;
561
562 gst = lst + lon / 15.0;
563
564 if (gst > 24.0)
565 gst -= 24.0;
566
567 jzjd = julian_date(1,0,yr);
568 ed = jd-jzjd;
569 t = (jzjd -2415020.0)/36525.0;
570 r = 6.6460656+2400.05126*t+2.58E-05*t*t;
571 b = 24-(r-24*(yr-1900));
572 t0 = ed * 0.0657098 - b;
573
574 if (t0 < 0.0)
575 t0 += 24;
576
577 gmt = gst-t0;
578
579 if (gmt<0)
580 gmt += 24.0;
581
582 gmt = gmt * 0.99727 - tz;;
583
584 if (gmt < 0)
585 gmt +=24.0;
586
587 dh_to_hm(gmt, h, m);
588}
589
590/*
591
592*/
593
594dh_to_hm(dh, h, m)
595double dh;
596int *h, *m;
597{
598 double tempsec;
599
600 *h = dh;
601 *m = (dh - *h) * 60;
602 tempsec = (dh - *h) * 60 - *m;
603 tempsec = tempsec * 60 + 0.5;
604
605 if (tempsec > 30)
606 (*m)++;
607
608 if (*m == 60) {
609
610 *m = 0;
611 (*h)++;
612 }
613}
614
615/*
616
617*/
618
619eq_to_altaz(r, d, t, alt, az)
620double r, d, t;
621double *alt, *az;
622{
623 double p = 3.14159265;
624 double r1 = p / 180.0;
625 double b = lat * r1;
626 double l = (360 - lon) * r1;
627 double t5, s1, c1, c2, s2, a, h;
628
629 r = r * 15.0 * r1;
630 d = d * r1;
631 t = t * 15.0 * r1;
632 t5 = t - r + l;
633 s1 = sin(b) * sin(d) + cos(b) * cos(d) * cos(t5);
634 c1 = 1 - s1 * s1;
635
636 if (c1 > 0) {
637
638 c1 = sqrt(c1);
639 h = atan(s1 / c1);
640
641 } else {
642
643 h = (s1 / fabs(s1)) * (p / 2.0);
644 }
645
646 c2 = cos(b) * sin(d) - sin(b) * cos(d) * cos(t5);
647 s2 = -cos(d) * sin(t5);
648
649 if (c2 == 0)
650 a = (s2/fabs(s2)) * (p/2);
651 else {
652
653 a = atan(s2/c2);
654
655 if (c2 < 0)
656 a=a+p;
657 }
658
659 if (a<0)
660 a=a+2*p;
661
662 *alt = h / r1;
663 *az = a / r1;
664}
665
666/*
667
668*/
669
670double
671gmst(j, f)
672double j,f;
673{
674 double d, j0, t, t1, t2, s;
675
676 d = j - 2451545.0;
677 t = d / 36525.0;
678 t1 = floor(t);
679 j0 = t1 * 36525.0 + 2451545.0;
680 t2 = (j - j0 + 0.5)/36525.0;
681 s = 24110.54841 + 184.812866 * t1;
682 s += 8640184.812866 * t2;
683 s += 0.093104 * t * t;
684 s -= 0.0000062 * t * t * t;
685 s /= 86400.0;
686 s -= floor(s);
687 s = 24 * (s + (f - 0.5) * 1.002737909);
688
689 if (s < 0)
690 s += 24.0;
691
692 if (s > 24.0)
693 s -= 24.0;
694
695 return(s);
696}
Note: See TracBrowser for help on using the repository browser.