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

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

Imported original source code.

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