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 y = cos(sdec) * sin(sra - lra);
364 x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
365 *pabl = atan22(y, x);
366 phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
367 i = atan22(sin(phi) , (dr - cos(phi)));
368 *ill = 0.5*(1 + cos(i));
369 }
370
371 /*
372 sunpos() takes days from J2000.0 and returns ecliptic longitude
373 of Sun in the pointers. Latitude is zero at this level of precision,
374 but pointer left in for consistency in number of arguments.
375 This function is within 0.01 degree (1 arcmin) almost all the time
376 for a century either side of J2000.0. This is from the 'low precision
377 fomulas for the Sun' from C24 of Astronomical Alamanac
378 */
sunpos(double d,double * lambda,double * beta,double * rvec)379 void sunpos(double d, double *lambda, double *beta, double *rvec) {
380 double L, g, ls, bs, rs;
381
382 L = range(280.461 + .9856474 * d) * RADS;
383 g = range(357.528 + .9856003 * d) * RADS;
384 ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
385 bs = 0;
386 rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
387 *lambda = ls;
388 *beta = bs;
389 *rvec = rs;
390 }
391
392 /*
393 this routine returns the altitude given the days since J2000.0
394 the hour angle and declination of the object and the latitude
395 of the observer. Used to find the Sun's altitude to put a letter
396 code on the transit time, and to find the Moon's altitude at
397 transit just to make sure that the Moon is visible.
398 */
alt(double glat,double ha,double dec)399 double alt(double glat, double ha, double dec) {
400 return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
401 }
402
403 /* returns an angle in degrees in the range 0 to 360 */
range(double x)404 double range(double x) {
405 double a, b;
406 b = x / 360;
407 a = 360 * (b - floor(b));
408 if (a < 0)
409 a = 360 + a;
410 return(a);
411 }
412
413 /* returns an angle in rads in the range 0 to two pi */
rangerad(double x)414 double rangerad(double x) {
415 double a, b;
416 b = x / TPI;
417 a = TPI * (b - floor(b));
418 if (a < 0)
419 a = TPI + a;
420 return(a);
421 }
422
423 /*
424 gets the atan2 function returning angles in the right
425 order and range
426 */
atan22(double y,double x)427 double atan22(double y, double x) {
428 double a;
429
430 a = atan2(y, x);
431 if (a < 0) a += TPI;
432 return(a);
433 }
434
435 /*
436 returns mean obliquity of ecliptic in radians given days since
437 J2000.0.
438 */
epsilon(double d)439 double epsilon(double d) {
440 double t = d/ 36525;
441 return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
442 }
443
444 /*
445 replaces ecliptic coordinates with equatorial coordinates
446 note: call by reference destroys original values
447 R is unchanged.
448 */
equatorial(double d,double * lon,double * lat,double * r)449 void equatorial(double d, double *lon, double *lat, double *r) {
450 double eps, ceps, seps, l, b;
451
452 l = *lon;
453 b = * lat;
454 eps = epsilon(d);
455 ceps = cos(eps);
456 seps = sin(eps);
457 *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
458 *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
459 }
460
461 /*
462 replaces equatorial coordinates with ecliptic ones. Inverse
463 of above, but used to find topocentric ecliptic coords.
464 */
ecliptic(double d,double * lon,double * lat,double * r)465 void ecliptic(double d, double *lon, double *lat, double *r) {
466 double eps, ceps, seps, alp, dec;
467 alp = *lon;
468 dec = *lat;
469 eps = epsilon(d);
470 ceps = cos(eps);
471 seps = sin(eps);
472 *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
473 *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
474 }
475
476 /*
477 returns the siderial time at greenwich meridian as
478 an angle in radians given the days since J2000.0
479 */
gst(double d)480 double gst( double d) {
481 double t = d / 36525;
482 double theta;
483 theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
484 return(theta * RADS);
485 }
486
tmoonsub_(double * day,double * glat,double * glong,double * moonalt,double * mrv,double * l,double * b,double * paxis)487 void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
488 double *mrv, double *l, double *b, double *paxis)
489 {
490 double mlambda, mbeta;
491 double malpha, mdelta;
492 double lst, mhr;
493 double tlambda, tbeta, trv;
494
495 lst = gst(*day) + *glong;
496
497 /* find Moon topocentric coordinates for libration calculations */
498
499 moonpos(*day, &mlambda, &mbeta, mrv);
500 malpha = mlambda;
501 mdelta = mbeta;
502 equatorial(*day, &malpha, &mdelta, mrv);
503 topo(lst, *glat, &malpha, &mdelta, mrv);
504 mhr = rangerad(lst - malpha);
505 *moonalt = alt(*glat, mhr, mdelta);
506
507 /* Optical libration and Position angle of the Pole */
508
509 tlambda = malpha;
510 tbeta = mdelta;
511 trv = *mrv;
512 ecliptic(*day, &tlambda, &tbeta, &trv);
513 libration(*day, tlambda, tbeta, malpha, l, b, paxis);
514 }
515