1 #include "astro.h"
2 
3 char*	herefile;
4 
5 int
6 main(int argc, char *argv[])
7 {
8 	int i, j;
9 	double d;
10 
11 	pi = atan(1.0)*4;
12 	pipi = pi*2;
13 	radian = pi/180;
14 	radsec = radian/3600;
15 	converge = 1.0e-14;
16 
17 	startab = unsharp("#9/sky/estartab");
18 	herefile = unsharp("#9/sky/here");
19 
20 	fmtinstall('R', Rconv);
21 	fmtinstall('D', Dconv);
22 
23 	per = PER;
24 	deld = PER/NPTS;
25 	init();
26 	args(argc, argv);
27 	init();
28 
29 loop:
30 	d = day;
31 	pdate(d);
32 	if(flags['p'] || flags['e']) {
33 		print(" ");
34 		ptime(d);
35 		pstime(d);
36 	}
37 	print("\n");
38 	for(i=0; i<=NPTS+1; i++) {
39 		setime(d);
40 
41 		for(j=0; objlst[j]; j++) {
42 			(*objlst[j]->obj)();
43 			setobj(&objlst[j]->point[i]);
44 			if(flags['p']) {
45 				if(flags['m'])
46 					if(strcmp(objlst[j]->name, "Comet"))
47 						continue;
48 				output(objlst[j]->name, &objlst[j]->point[i]);
49 			}
50 		}
51 		if(flags['e']) {
52 			d = dist(&eobj1->point[i], &eobj2->point[i]);
53 			print("dist %s to %s = %.4f\n", eobj1->name, eobj2->name, d);
54 		}
55 /*		if(flags['p']) { */
56 /*			pdate(d); */
57 /*			print(" "); */
58 /*			ptime(d); */
59 /*			print("\n"); */
60 /*		} */
61 		if(flags['p'] || flags['e'])
62 			break;
63 		d += deld;
64 	}
65 	if(!(flags['p'] || flags['e']))
66 		search();
67 	day += per;
68 	nperiods -= 1;
69 	if(nperiods > 0)
70 		goto loop;
71 	exits(0);
72 	return 0;		/* gcc */
73 }
74 
75 void
76 args(int argc, char *argv[])
77 {
78 	char *p;
79 	long t;
80 	int f, i;
81 	Obj2 *q;
82 
83 	memset(flags, 0, sizeof(flags));
84 	ARGBEGIN {
85 	default:
86 		fprint(2, "astro [-adeklmopst] [-c nperiod] [-C tperiod]\n");
87 		exits(0);
88 
89 	case 'c':
90 		nperiods = 1;
91 		p = ARGF();
92 		if(p)
93 			nperiods = atol(p);
94 		flags['c']++;
95 		break;
96 	case 'C':
97 		p = ARGF();
98 		if(p)
99 			per = atof(p);
100 		break;
101 	case 'e':
102 		eobj1 = nil;
103 		eobj2 = nil;
104 		p = ARGF();
105 		if(p) {
106 			for(i=0; q=objlst[i]; i++) {
107 				if(strcmp(q->name, p) == 0)
108 					eobj1 = q;
109 				if(strcmp(q->name1, p) == 0)
110 					eobj1 = q;
111 			}
112 			p = ARGF();
113 			if(p) {
114 				for(i=0; q=objlst[i]; i++) {
115 					if(strcmp(q->name, p) == 0)
116 						eobj2 = q;
117 					if(strcmp(q->name1, p) == 0)
118 						eobj2 = q;
119 				}
120 			}
121 		}
122 		if(eobj1 && eobj2) {
123 			flags['e']++;
124 			break;
125 		}
126 		fprint(2, "cant recognize eclipse objects\n");
127 		exits("eflag");
128 
129 	case 'a':
130 	case 'd':
131 	case 'j':
132 	case 'k':
133 	case 'l':
134 	case 'm':
135 	case 'o':
136 	case 'p':
137 	case 's':
138 	case 't':
139 		flags[ARGC()]++;
140 		break;
141 	} ARGEND
142 	if(*argv){
143 		fprint(2, "usage: astro [-dlpsatokm] [-c nday] [-e obj1 obj2]\n");
144 		exits("usage");
145 	}
146 
147 	t = time(0);
148 	day = t/86400. + 25567.5;
149 	if(flags['d'])
150 		day = readate();
151 	if(flags['j'])
152 		print("jday = %.4f\n", day);
153 	deltat = day * .001704;
154 	if(deltat > 32.184)		/* assume date is utc1 */
155 		deltat = 32.184;	/* correct by leap sec */
156 	if(flags['t'])
157 		deltat = readdt();
158 
159 	if(flags['l']) {
160 		fprint(2, "nlat wlong elev\n");
161 		readlat(0);
162 	} else {
163 		f = open(herefile, OREAD);
164 		if(f < 0) {
165 			fprint(2, "%s?\n", herefile);
166 			/* btl mh */
167 			nlat = (40 + 41.06/60)*radian;
168 			awlong = (74 + 23.98/60)*radian;
169 			elev = 150 * 3.28084;
170 		} else {
171 			readlat(f);
172 			close(f);
173 		}
174 	}
175 }
176 
177 double
178 readate(void)
179 {
180 	int i;
181 	Tim t;
182 
183 	fprint(2, "year mo da hr min\n");
184 	rline(0);
185 	for(i=0; i<5; i++)
186 		t.ifa[i] = atof(skip(i));
187 	return convdate(&t);
188 }
189 
190 double
191 readdt(void)
192 {
193 
194 	fprint(2, "ΔT (sec) (%.3f)\n", deltat);
195 	rline(0);
196 	return atof(skip(0));
197 }
198 
199 double
200 etdate(long year, int mo, double day)
201 {
202 	Tim t;
203 
204 	t.ifa[0] = year;
205 	t.ifa[1] = mo;
206 	t.ifa[2] = day;
207 	t.ifa[3] = 0;
208 	t.ifa[4] = 0;
209 	return convdate(&t) + 2415020;
210 }
211 
212 void
213 readlat(int f)
214 {
215 
216 	rline(f);
217 	nlat = atof(skip(0)) * radian;
218 	awlong = atof(skip(1)) * radian;
219 	elev = atof(skip(2)) * 3.28084;
220 }
221 
222 double
223 fmod(double a, double b)
224 {
225 	return a - floor(a/b)*b;
226 }