1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4
5 #define RADS 0.0174532925199433
6 #define DEGS 57.2957795130823
7 #define TPI 6.28318530717959
8 #define PI 3.1415927
9
10 /* ratio of earth radius to astronomical unit */
11 #define ER_OVER_AU 0.0000426352325194252
12
13 /* all prototypes here */
14
15 double getcoord(int coord);
16 void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
17 double range(double y);
18 double rangerad(double y);
19 double days(int y, int m, int dn, double hour);
20 double days_(int *y, int *m, int *dn, double *hour);
21 void moonpos(double, double *, double *, double *);
22 void sunpos(double , double *, double *, double *);
23 double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
24 double atan22(double y, double x);
25 double epsilon(double d);
26 void equatorial(double d, double *lon, double *lat, double *r);
27 void ecliptic(double d, double *lon, double *lat, double *r);
28 double gst(double d);
29 void topo(double lst, double glat, double *alp, double *dec, double *r);
30 double alt(double glat, double ha, double dec);
31 void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
32 void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
33 int daysinmonth(int y, int m);
34 int isleap(int y);
35 void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
36 double *mrv, double *l, double *b, double *paxis);
37
38 static const char
39 usage[] = " Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
40 "example: tmoon 200009 0 -00155 5230\n";
41
42 /*
43 getargs() gets the arguments from the command line, does some basic error
44 checking, and converts arguments into numerical form. Arguments are passed
45 back in pointers. Error messages print to stderr so re-direction of output
46 to file won't leave users blind. Error checking prints list of all errors
47 in a command line before quitting.
48 */
getargs(int argc,char * argv[],int * y,int * m,double * tz,double * glong,double * glat)49 void getargs(int argc, char *argv[], int *y, int *m, double *tz,
50 double *glong, double *glat) {
51
52 int date, latitude, longitude;
53 int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
54 int longminflag = 0, latminflag = 0, dflag = 0;
55
56 /* if not right number of arguments, then print example command line */
57
58 if (argc !=5) {
59 fprintf(stderr, usage);
60 exit(EXIT_FAILURE);
61 }
62
63 date = atoi(argv[1]);
64 *y = date / 100;
65 *m = date - *y * 100;
66 *tz = (double) atof(argv[2]);
67 longitude = atoi(argv[3]);
68 latitude = atoi(argv[4]);
69 *glong = RADS * getcoord(longitude);
70 *glat = RADS * getcoord(latitude);
71
72 /* set a flag for each error found */
73
74 if (*m > 12 || *m < 1) mflag = 1;
75 if (*y > 2500) yflag = 1;
76 if (date < 150001) dflag = 1;
77 if (fabs((float) *glong) > 180 * RADS) longflag = 1;
78 if (abs(longitude) % 100 > 59) longminflag = 1;
79 if (fabs((float) *glat) > 90 * RADS) latflag = 1;
80 if (abs(latitude) % 100 > 59) latminflag = 1;
81 if (fabs((float) *tz) > 12) tzflag = 1;
82
83 /* print all the errors found */
84
85 if (dflag == 1) {
86 fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
87 }
88 if (yflag == 1) {
89 fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
90 }
91 if (mflag == 1) {
92 fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
93 }
94 if (tzflag == 1) {
95 fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
96 }
97 if (longflag == 1) {
98 fprintf(stderr, "long: must be in range +/- 180 degrees\n");
99 }
100 if (longminflag == 1) {
101 fprintf(stderr, "long: last two digits are arcmin - max 59\n");
102 }
103 if (latflag == 1) {
104 fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
105 }
106 if (latminflag == 1) {
107 fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
108 }
109
110 /* quits if one or more flags set */
111
112 if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
113 exit(EXIT_FAILURE);
114 }
115
116 }
117
118 /*
119 returns coordinates in decimal degrees given the
120 coord as a ddmm value stored in an integer.
121 */
getcoord(int coord)122 double getcoord(int coord) {
123 int west = 1;
124 double glg, deg;
125 if (coord < 0) west = -1;
126 glg = fabs((double) coord/100);
127 deg = floor(glg);
128 glg = west* (deg + (glg - deg)*100 / 60);
129 return(glg);
130 }
131
132 /*
133 days() takes the year, month, day in the month and decimal hours
134 in the day and returns the number of days since J2000.0.
135 Assumes Gregorian calendar.
136 */
days(int y,int m,int d,double h)137 double days(int y, int m, int d, double h) {
138 int a, b;
139 double day;
140
141 /*
142 The lines below work from 1900 march to feb 2100
143 a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
144 day = (double)a - 730531.5 + hour / 24;
145 */
146
147 /* These lines work for any Gregorian date since 0 AD */
148 if (m ==1 || m==2) {
149 m +=12;
150 y -= 1;
151 }
152 a = y / 100;
153 b = 2 - a + a/4;
154 day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
155 + d + b - 1524.5 - 2451545 + h/24;
156 return(day);
157 }
days_(int * y0,int * m0,int * d0,double * h0)158 double days_(int *y0, int *m0, int *d0, double *h0)
159 {
160 return days(*y0,*m0,*d0,*h0);
161 }
162
163 /*
164 Returns 1 if y a leap year, and 0 otherwise, according
165 to the Gregorian calendar
166 */
isleap(int y)167 int isleap(int y) {
168 int a = 0;
169 if(y % 4 == 0) a = 1;
170 if(y % 100 == 0) a = 0;
171 if(y % 400 == 0) a = 1;
172 return(a);
173 }
174
175 /*
176 Given the year and the month, function returns the
177 number of days in the month. Valid for Gregorian
178 calendar.
179 */
daysinmonth(int y,int m)180 int daysinmonth(int y, int m) {
181 int b = 31;
182 if(m == 2) {
183 if(isleap(y) == 1) b= 29;
184 else b = 28;
185 }
186 if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
187 return(b);
188 }
189
190 /*
191 moonpos() takes days from J2000.0 and returns ecliptic coordinates
192 of moon in the pointers. Note call by reference.
193 This function is within a couple of arcminutes most of the time,
194 and is truncated from the Meeus Ch45 series, themselves truncations of
195 ELP-2000. Returns moon distance in earth radii.
196 Terms have been written out explicitly rather than using the
197 table based method as only a small number of terms is
198 retained.
199 */
moonpos(double d,double * lambda,double * beta,double * rvec)200 void moonpos(double d, double *lambda, double *beta, double *rvec) {
201 double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;
202
203 t = d / 36525;
204
205 L = range(218.3164591 + 481267.88134236 * t) * RADS;
206 D = range(297.8502042 + 445267.1115168 * t) * RADS;
207 M = range(357.5291092 + 35999.0502909 * t) * RADS;
208 M1 = range(134.9634114 + 477198.8676313 * t - .008997 * t * t) * RADS;
209 F = range(93.27209929999999 + 483202.0175273 * t - .0034029*t*t)*RADS;
210 e = 1 - .002516 * t;
211
212 dl = 6288774 * sin(M1);
213 dl += 1274027 * sin(2 * D - M1);
214 dl += 658314 * sin(2 * D);
215 dl += 213618 * sin(2 * M1);
216 dl -= e * 185116 * sin(M);
217 dl -= 114332 * sin(2 * F) ;
218 dl += 58793 * sin(2 * D - 2 * M1);
219 dl += e * 57066 * sin(2 * D - M - M1) ;
220 dl += 53322 * sin(2 * D + M1);
221 dl += e * 45758 * sin(2 * D - M);
222 dl -= e * 40923 * sin(M - M1);
223 dl -= 34720 * sin(D) ;
224 dl -= e * 30383 * sin(M + M1) ;
225 dl += 15327 * sin(2 * D - 2 * F) ;
226 dl -= 12528 * sin(M1 + 2 * F);
227 dl += 10980 * sin(M1 - 2 * F);
228 lm = rangerad(L + dl / 1000000 * RADS);
229
230 dB = 5128122 * sin(F);
231 dB += 280602 * sin(M1 + F);
232 dB += 277693 * sin(M1 - F);
233 dB += 173237 * sin(2 * D - F);
234 dB += 55413 * sin(2 * D - M1 + F);
235 dB += 46271 * sin(2 * D - M1 - F);
236 dB += 32573 * sin(2 * D + F);
237 dB += 17198 * sin(2 * M1 + F);
238 dB += 9266 * sin(2 * D + M1 - F);
239 dB += 8822 * sin(2 * M1 - F);
240 dB += e * 8216 * sin(2 * D - M - F);
241 dB += 4324 * sin(2 * D - 2 * M1 - F);
242 bm = dB / 1000000 * RADS;
243
244 dR = -20905355 * cos(M1);
245 dR -= 3699111 * cos(2 * D - M1);
246 dR -= 2955968 * cos(2 * D);
247 dR -= 569925 * cos(2 * M1);
248 dR += e * 48888 * cos(M);
249 dR -= 3149 * cos(2 * F);
250 dR += 246158 * cos(2 * D - 2 * M1);
251 dR -= e * 152138 * cos(2 * D - M - M1) ;
252 dR -= 170733 * cos(2 * D + M1);
253 dR -= e * 204586 * cos(2 * D - M);
254 dR -= e * 129620 * cos(M - M1);
255 dR += 108743 * cos(D);
256 dR += e * 104755 * cos(M + M1);
257 dR += 79661 * cos(M1 - 2 * F);
258 rm = 385000.56 + dR / 1000;
259
260 *lambda = lm;
261 *beta = bm;
262 /* distance to Moon must be in Earth radii */
263 *rvec = rm / 6378.14;
264 }
265
266 /*
267 topomoon() takes the local siderial time, the geographical
268 latitude of the observer, and pointers to the geocentric
269 equatorial coordinates. The function overwrites the geocentric
270 coordinates with topocentric coordinates on a simple spherical
271 earth model (no polar flattening). Expects Moon-Earth distance in
272 Earth radii. Formulas scavenged from Astronomical Almanac 'low
273 precision formulae for Moon position' page D46.
274 */
275
topo(double lst,double glat,double * alp,double * dec,double * r)276 void topo(double lst, double glat, double *alp, double *dec, double *r) {
277 double x, y, z, r1;
278 x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
279 y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
280 z = *r * sin(*dec) - sin(glat);
281 r1 = sqrt(x*x + y*y + z*z);
282 *alp = atan22(y, x);
283 *dec = asin(z / r1);
284 *r = r1;
285 }
286
287 /*
288 moontransit() takes date, the time zone and geographic longitude
289 of observer and returns the time (decimal hours) of lunar transit
290 on that day if there is one, and sets the notransit flag if there
291 isn't. See Explanatory Supplement to Astronomical Almanac
292 section 9.32 and 9.31 for the method.
293 */
294
moontransit(int y,int m,int d,double tz,double glat,double glong,int * notransit)295 double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
296 double hm, ht, ht1, lon, lat, rv, dnew, lst;
297 int itcount;
298
299 ht1 = 180 * RADS;
300 ht = 0;
301 itcount = 0;
302 *notransit = 0;
303 do {
304 ht = ht1;
305 itcount++;
306 dnew = days(y, m, d, ht * DEGS/15) - tz/24;
307 lst = gst(dnew) + glong;
308 /* find the topocentric Moon ra (hence hour angle) and dec */
309 moonpos(dnew, &lon, &lat, &rv);
310 equatorial(dnew, &lon, &lat, &rv);
311 topo(lst, glat, &lon, &lat, &rv);
312 hm = rangerad(lst - lon);
313 ht1 = rangerad(ht - hm);
314 /* if no convergence, then no transit on that day */
315 if (itcount > 30) {
316 *notransit = 1;
317 break;
318 }
319 }
320 while (fabs(ht - ht1) > 0.04 * RADS);
321 return(ht1);
322 }
323
324 /*
325 Calculates the selenographic coordinates of either the sub Earth point
326 (optical libration) or the sub-solar point (selen. coords of centre of
327 bright hemisphere). Based on Meeus chapter 51 but neglects physical
328 libration and nutation, with some simplification of the formulas.
329 */
libration(double day,double lambda,double beta,double alpha,double * l,double * b,double * p)330 void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
331 double i, f, omega, w, y, x, a, t, eps;
332 t = day / 36525;
333 i = 1.54242 * RADS;
334 eps = epsilon(day);
335 f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
336 omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
337 w = lambda - omega;
338 y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
339 x = cos(w) * cos(beta);
340 a = atan22(y, x);
341 *l = a - f;
342
343 /* kludge to catch cases of 'round the back' angles */
344 if (*l < -90 * RADS) *l += TPI;
345 if (*l > 90 * RADS) *l -= TPI;
346 *b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));
347
348 /* pa pole axis - not used for Sun stuff */
349 x = sin(i) * sin(omega);
350 y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
351 w = atan22(x, y);
352 *p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
353 }
354
355 /*
356 Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
357 eq coords Sun
358 Returns: position angle of bright limb wrt NCP, percentage illumination
359 of Sun
360 */
illumination(double day,double lra,double ldec,double dr,double sra,double sdec,double * pabl,double * ill)361 void illumination(double day , double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
362 double x, y, phi, i;
363 (void)day;
364 y = cos(sdec) * sin(sra - lra);
365 x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
366 *pabl = atan22(y, x);
367 phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
368 i = atan22(sin(phi) , (dr - cos(phi)));
369 *ill = 0.5*(1 + cos(i));
370 }
371
372 /*
373 sunpos() takes days from J2000.0 and returns ecliptic longitude
374 of Sun in the pointers. Latitude is zero at this level of precision,
375 but pointer left in for consistency in number of arguments.
376 This function is within 0.01 degree (1 arcmin) almost all the time
377 for a century either side of J2000.0. This is from the 'low precision
378 fomulas for the Sun' from C24 of Astronomical Alamanac
379 */
sunpos(double d,double * lambda,double * beta,double * rvec)380 void sunpos(double d, double *lambda, double *beta, double *rvec) {
381 double L, g, ls, bs, rs;
382
383 L = range(280.461 + .9856474 * d) * RADS;
384 g = range(357.528 + .9856003 * d) * RADS;
385 ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
386 bs = 0;
387 rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
388 *lambda = ls;
389 *beta = bs;
390 *rvec = rs;
391 }
392
393 /*
394 this routine returns the altitude given the days since J2000.0
395 the hour angle and declination of the object and the latitude
396 of the observer. Used to find the Sun's altitude to put a letter
397 code on the transit time, and to find the Moon's altitude at
398 transit just to make sure that the Moon is visible.
399 */
alt(double glat,double ha,double dec)400 double alt(double glat, double ha, double dec) {
401 return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
402 }
403
404 /* returns an angle in degrees in the range 0 to 360 */
range(double x)405 double range(double x) {
406 double a, b;
407 b = x / 360;
408 a = 360 * (b - floor(b));
409 if (a < 0)
410 a = 360 + a;
411 return(a);
412 }
413
414 /* returns an angle in rads in the range 0 to two pi */
rangerad(double x)415 double rangerad(double x) {
416 double a, b;
417 b = x / TPI;
418 a = TPI * (b - floor(b));
419 if (a < 0)
420 a = TPI + a;
421 return(a);
422 }
423
424 /*
425 gets the atan2 function returning angles in the right
426 order and range
427 */
atan22(double y,double x)428 double atan22(double y, double x) {
429 double a;
430
431 a = atan2(y, x);
432 if (a < 0) a += TPI;
433 return(a);
434 }
435
436 /*
437 returns mean obliquity of ecliptic in radians given days since
438 J2000.0.
439 */
epsilon(double d)440 double epsilon(double d) {
441 double t = d/ 36525;
442 return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
443 }
444
445 /*
446 replaces ecliptic coordinates with equatorial coordinates
447 note: call by reference destroys original values
448 R is unchanged.
449 */
equatorial(double d,double * lon,double * lat,double * r)450 void equatorial(double d, double *lon, double *lat, double * r) {
451 double eps, ceps, seps, l, b;
452 (void)r;
453
454 l = *lon;
455 b = * lat;
456 eps = epsilon(d);
457 ceps = cos(eps);
458 seps = sin(eps);
459 *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
460 *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
461 }
462
463 /*
464 replaces equatorial coordinates with ecliptic ones. Inverse
465 of above, but used to find topocentric ecliptic coords.
466 */
ecliptic(double d,double * lon,double * lat,double * r)467 void ecliptic(double d, double *lon, double *lat, double * r) {
468 double eps, ceps, seps, alp, dec;
469 (void)r;
470
471 alp = *lon;
472 dec = *lat;
473 eps = epsilon(d);
474 ceps = cos(eps);
475 seps = sin(eps);
476 *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
477 *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
478 }
479
480 /*
481 returns the siderial time at greenwich meridian as
482 an angle in radians given the days since J2000.0
483 */
gst(double d)484 double gst( double d) {
485 double t = d / 36525;
486 double theta;
487 theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
488 return(theta * RADS);
489 }
490
tmoonsub_(double * day,double * glat,double * glong,double * moonalt,double * mrv,double * l,double * b,double * paxis)491 void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
492 double *mrv, double *l, double *b, double *paxis)
493 {
494 double mlambda, mbeta;
495 double malpha, mdelta;
496 double lst, mhr;
497 double tlambda, tbeta, trv;
498
499 lst = gst(*day) + *glong;
500
501 /* find Moon topocentric coordinates for libration calculations */
502
503 moonpos(*day, &mlambda, &mbeta, mrv);
504 malpha = mlambda;
505 mdelta = mbeta;
506 equatorial(*day, &malpha, &mdelta, mrv);
507 topo(lst, *glat, &malpha, &mdelta, mrv);
508 mhr = rangerad(lst - malpha);
509 *moonalt = alt(*glat, mhr, mdelta);
510
511 /* Optical libration and Position angle of the Pole */
512
513 tlambda = malpha;
514 tbeta = mdelta;
515 trv = *mrv;
516 ecliptic(*day, &tlambda, &tbeta, &trv);
517 libration(*day, tlambda, tbeta, malpha, l, b, paxis);
518 }
519