1 // $Id: Skycal.cc 5748 2014-10-11 19:38:53Z flaterco $
2
3 // Skycal.cc -- Functions for sun and moon events.
4
5 // Prediction of moon phases, sun and moon rises and sets has nothing
6 // to do with tide prediction. There is no overlap between this code
7 // and the tide prediction code.
8
9 // This source file began its life as skycalendar.c and skycalc.c in
10 // John Thorstensen's skycal distribution (version 4.1, 1994-09) at
11 // ftp://iraf.noao.edu/contrib/skycal.tar.Z. Those portions that are
12 // unchanged from the original sources are covered by the original
13 // license statement, included below. The new portions and "value
14 // added" by David Flater are covered under the GNU General Public
15 // License.
16
17 // 2013-05-26
18 //
19 // Patch from James Ashton to replace etcorr.
20
21 // 2003-02-04
22 //
23 // The release notes for Skycal V5 indicated that some code in 4.1 had
24 // copyright problems. Reviewed all code snippets here and found that
25 // none were impacted by copyright-related changes.
26 //
27 // Harmonized atan_circ with slightly improved V5 version.
28
29 /*
30 Copyright (C) 1998 David Flater.
31
32 This program is free software: you can redistribute it and/or modify
33 it under the terms of the GNU General Public License as published by
34 the Free Software Foundation, either version 3 of the License, or
35 (at your option) any later version.
36
37 This program is distributed in the hope that it will be useful,
38 but WITHOUT ANY WARRANTY; without even the implied warranty of
39 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
40 GNU General Public License for more details.
41
42 You should have received a copy of the GNU General Public License
43 along with this program. If not, see <http://www.gnu.org/licenses/>.
44 */
45
46 // The relevant portions of the Skycal 5 comments and license statement
47 // are as follows:
48
49 /* SKY CALCULATOR PROGRAM
50 John Thorstensen, Dartmouth College.
51 This program computes many quantities frequently needed by the
52 observational astronomer. It is written as a completely
53 self-contained program in standard c, so it should be
54 very transportable; the only issue I know of that really affects
55 portability is the adequacy of the double-precision floating
56 point accuracy on the machine. Experience shows that c compilers
57 on various systems have idiosyncracies, though, so be sure
58 to check carefully.
59
60 This is intended as an observatory utility program; I assume the
61 user is familiar with astronomical coordinates and nomenclature.
62 While the code should be very transportable, I also
63 assume it will be installed by a conscientious person who
64 will run critical tests before it is released at a new site.
65 Experience shows that some c compilers generate unforseen errors
66 when the code is ported, so the output should be checked meticulously
67 against data from other sites.
68
69 [...]
70
71 The program is self-contained. It has been developed primarily on
72 UNIX and Linux machines, and should adapt easily to any system with
73 a c compiler.
74
75 ** BUT CAUTION ... **
76 Because many of the routines take a double-precision floating point
77 Julian Date as their time argument, one must be sure that the machine
78 and compiler carry sufficient mantissa to reach the desired accuracy.
79 On most architectures the double-precision floating point julian date
80 has an accuracy of order 0.01 seconds of time, which is just adequate.
81
82 LEGALITIES:
83
84 I make no guarantee as to the accuracy, reliability, or
85 appropriateness of this program, though I have found it to be
86 reasonably accurate and quite useful to the working astronomer.
87
88 The program is COPYRIGHT 2000 BY JOHN THORSTENSEN.
89 Permission is hereby granted for non-profit scientific or educational use.
90 For-profit use (e. g., by astrologers!) must be through negotiated
91 license. The author requests that observatories and astronomy
92 departments which install this as a utility notify the author
93 by paper mail, just so I know how widely it is used.
94
95 Credits:
96 * The julian date and sidereal time routines were
97 originally coded in PL/I by Steve Maker of Dartmouth College.
98 They were based on routines in the old American Ephemeris.
99 Many of the routines were coded from Jean Meeus' "Astronomical
100 Formulae for Calculators", published by Willman-Bell. This is
101 an extraordinarily helpful little book!
102 */
103
104 // The Skycal 4.1 comments and license statement are as follows:
105
106 /*
107 This is a self-contained c-language program to print a nighttime
108 astronomical calendar for use in planning observations.
109 It prints to standard output (usually the terminal); the
110 operator should capture this output (e. g., using redirection
111 in UNIX or the /out= switch in VMS) and then print it on an
112 appropriate output device. The table which is printed is some
113 125 columns wide, so a wide device is required (either a line
114 printer or a laserprinter in LANDSCAPE mode.) It is assumed that
115 the ASCII form-feed character will actually begin a new page.
116 The original program was to run on VMS, but it should be very
117 transportable. Non-vms users will probably want to change
118 'unixio.h' to 'stdio.h' in the first line.
119 An explanatory text is printed at the beginning of the output, which
120 includes the appropriate CAUTIONS regarding accuracy and applicability.
121
122 A number of 'canned site' parameters have been included. Be
123 careful of time zones, DST etc. for foreign sites.
124 To customize to your own site, install an option in the
125 routine 'load_site'. The code is very straightforward; just do
126 it exactly the same as the others. You might also want to erase
127 some seldom-used choices. One can also specify new site parameters
128 at run time.
129
130 This program may be used freely by anyone for scientific or educational
131 purposes. If you use it for profit, I want a cut, and claim
132 a copyright herewith. In any case please acknowledge the source:
133
134 John Thorstensen
135 Dept. of Physics and Astronomy
136 Dartmouth College
137 Hanover, NH 03755
138 John.Thorstensen@dartmouth.edu
139
140 May 26, 1993.
141 */
142
143 #include "libxtide.hh"
144 #include "Skycal.hh"
145
146 // To avoid creating a clash on every reference of a math.h trig function,
147 // namespace libxtide is applied only to the non-static definitions.
148
149
150 // DWF: This affects the rise/set predictions. Normally you would
151 // need to adjust it for the elevation of the location, but since this
152 // is a tide prediction program, we can safely assume that we are
153 // always at sea level :-)
154 static const double riseAltitude (-0.83);
155
156 #define DEG_IN_RADIAN 57.2957795130823
157 #define HRS_IN_RADIAN 3.819718634
158 #define SEC_IN_DAY 86400.
159 #define GREG_DAYS_IN_YEAR 365.2425 /* 365 + 97 / 400 */
160 #define FLATTEN 0.003352813 /* flattening of earth, 1/298.257 */
161 #define EQUAT_RAD 6378137. /* equatorial radius of earth, meters */
162 #define J2000 2451545. /* Julian date at standard epoch */
163
164
165 // Harmonized with Skycal V5 2003-02-04
atan_circ(double x,double y)166 static double atan_circ(double x, double y)
167 {
168 /* returns radian angle 0 to 2pi for coords x, y --
169 get that quadrant right !! */
170
171 double theta;
172
173 if((x == 0.) && (y == 0.)) return(0.); /* guard ... */
174
175 theta = atan2(y,x); /* turns out there is such a thing in math.h */
176 while(theta < 0.) theta += 2.0 * M_PI;
177 return(theta);
178 }
179
180
181 // This is stripped down from Skycal V4.1 and even more stripped down
182 // from Skycal V5.
altit(double dec,double ha,double lat)183 static double altit (double dec, double ha, double lat)
184 /* dec deg, dec hrs, dec deg */
185 {
186 double x;
187 dec = dec / DEG_IN_RADIAN;
188 ha = ha / HRS_IN_RADIAN;
189 lat = lat / DEG_IN_RADIAN; /* thank heavens for pass-by-value */
190 x = DEG_IN_RADIAN * asin(cos(dec)*cos(ha)*cos(lat) + sin(dec)*sin(lat));
191 return(x);
192 }
193
194
195 // No important changes in Skycal V5.
lst(double jd,double longit)196 static double lst(double jd, double longit)
197
198 {
199 /* returns the local MEAN sidereal time (dec hrs) at julian date jd
200 at west longitude long (decimal hours). Follows
201 definitions in 1992 Astronomical Almanac, pp. B7 and L2.
202 Expression for GMST at 0h ut referenced to Aoki et al, A&A 105,
203 p.359, 1982. */
204
205 double t, ut, jdmid, jdint, jdfrac, sid_g;
206 long jdin, sid_int;
207
208 jdin = (long)jd; /* fossil code from earlier package which
209 split jd into integer and fractional parts ... */
210 jdint = jdin;
211 jdfrac = jd - jdint;
212 if(jdfrac < 0.5) {
213 jdmid = jdint - 0.5;
214 ut = jdfrac + 0.5;
215 }
216 else {
217 jdmid = jdint + 0.5;
218 ut = jdfrac - 0.5;
219 }
220 t = (jdmid - J2000)/36525;
221 sid_g = (24110.54841+8640184.812866*t+0.093104*t*t-6.2e-6*t*t*t)/86400.;
222 sid_int = (long)sid_g;
223 sid_g = sid_g - (double) sid_int;
224 sid_g = sid_g + 1.0027379093 * ut - longit/24.;
225 sid_int = (long)sid_g;
226 sid_g = (sid_g - (double) sid_int) * 24.;
227 if(sid_g < 0.) sid_g = sid_g + 24.;
228 return(sid_g);
229 }
230
231
232 // No important changes in Skycal V5.
233 static void
lpsun(double jd,double * ra,double * dec)234 lpsun(double jd, double *ra, double *dec)
235
236 /* Low precision formulae for the sun, from Almanac p. C24 (1990) */
237 /* ra and dec are returned as decimal hours and decimal degrees. */
238
239 {
240 double n, L, g, lambda,epsilon,x,y,z;
241
242 n = jd - J2000;
243 L = 280.460 + 0.9856474 * n;
244 g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
245 lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
246 epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
247
248 x = cos(lambda);
249 y = cos(epsilon) * sin(lambda);
250 z = sin(epsilon)*sin(lambda);
251
252 *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
253 *dec = (asin(z))*DEG_IN_RADIAN;
254 }
255
256
257 // Fwd Decl
258 static void accumoon (double jd, double geolat, double lst, double elevsea,
259 double *topora, double *topodec, double *topodist);
260
261
262 // Not in Skycal.
263 // This function combines a few steps that are always used together to
264 // find the altitude of the sun or moon.
altitude(double jd,double lat,double longit,bool lunar)265 static double altitude (double jd, double lat, double longit, bool lunar) {
266 if (lunar) {
267 double ra, dec, dist, sid;
268 sid = lst(jd,longit);
269 accumoon (jd, lat, sid, 0, &ra, &dec, &dist);
270 return altit(dec, sid-ra, lat);
271 } else {
272 double ra, dec;
273 lpsun(jd, &ra, &dec);
274 return altit(dec, lst(jd,longit)-ra, lat);
275 }
276 }
277
278
279 // I messed with this a good bit.
280 //
281 // * It now converges to within 1 minute.
282 // * It converges better from bad initial guesses. (Deriv is now
283 // updated inside of the loop.)
284 // * It won't roam more than half a day in either direction.
285 // * Max iterations chosen conservatively.
286 // * It finishes with a check to determine what it found.
287
288 // 2003-02-04
289 // Expanded to handle moonrise/moonset (replacing jd_moon_alt too).
290
291 // No important changes in Skycal V5.
jd_alt(double alt,double jdorig,double lat,double longit,bool lunar,bool & is_rise)292 static double jd_alt (double alt, double jdorig, double lat, double longit,
293 bool lunar, bool &is_rise)
294 {
295 /* returns jd at which sun/moon is at a given
296 altitude, given jdguess as a starting point. */
297
298 double jdguess = jdorig;
299 double jdout, adj = 1.0;
300 double deriv, err, del = 0.002;
301 double alt2,alt3;
302 short i = 0;
303
304 /* first guess */
305
306 alt2 = altitude (jdguess, lat, longit, lunar);
307 jdguess = jdguess + del;
308 alt3 = altitude (jdguess, lat, longit, lunar);
309 err = alt3 - alt;
310 deriv = (alt3 - alt2) / del;
311 if (deriv == 0.0)
312 return (-1.0e10); // Found dead end.
313 adj = -err/deriv;
314 while(fabs(adj) >= libxtide::Global::eventPrecisionJD) {
315 if (i++ == 12)
316 return(-1.0e10); // Exceeded max iterations.
317 jdguess += adj;
318 if (fabs (jdguess - jdorig) > 0.5)
319 return (-1.0e10); // Ran out of bounds.
320 alt2 = alt3;
321 alt3 = altitude (jdguess, lat, longit, lunar);
322 err = alt3 - alt;
323 deriv = (alt3 - alt2) / adj;
324 if (deriv == 0.0)
325 return (-1.0e10); // Found dead end.
326 adj = -err/deriv;
327 }
328 jdout = jdguess;
329
330 // Figure out whether this is a rise or a set by shifting
331 // by 1 second.
332 {
333 jdguess -= 1.0 / SEC_IN_DAY;
334 alt2 = altitude (jdguess, lat, longit, lunar);
335 is_rise = (alt2 < alt3);
336 }
337
338 return(jdout);
339 }
340
341
342 // No important changes in Skycal V5.
flmoon(int n,int nph,double * jdout)343 static void flmoon(int n, int nph, double *jdout)
344
345 /* Gives jd (+- 2 min) of phase nph on lunation n; replaces
346 less accurate Numerical Recipes routine. This routine
347 implements formulae found in Jean Meeus' *Astronomical Formulae
348 for Calculators*, 2nd edition, Willman-Bell. A very useful
349 book!! */
350
351 {
352 double jd, cor;
353 double M, Mpr, F;
354 double T;
355 double lun;
356
357 lun = (double) n + (double) nph / 4.;
358 T = lun / 1236.85;
359 jd = 2415020.75933 + 29.53058868 * lun
360 + 0.0001178 * T * T
361 - 0.000000155 * T * T * T
362 + 0.00033 * sin((166.56 + 132.87 * T - 0.009173 * T * T)/DEG_IN_RADIAN);
363 M = 359.2242 + 29.10535608 * lun - 0.0000333 * T * T - 0.00000347 * T * T * T;
364 M = M / DEG_IN_RADIAN;
365 Mpr = 306.0253 + 385.81691806 * lun + 0.0107306 * T * T + 0.00001236 * T * T * T;
366 Mpr = Mpr / DEG_IN_RADIAN;
367 F = 21.2964 + 390.67050646 * lun - 0.0016528 * T * T - 0.00000239 * T * T * T;
368 F = F / DEG_IN_RADIAN;
369 if((nph == 0) || (nph == 2)) {/* new or full */
370 cor = (0.1734 - 0.000393*T) * sin(M)
371 + 0.0021 * sin(2*M)
372 - 0.4068 * sin(Mpr)
373 + 0.0161 * sin(2*Mpr)
374 - 0.0004 * sin(3*Mpr)
375 + 0.0104 * sin(2*F)
376 - 0.0051 * sin(M + Mpr)
377 - 0.0074 * sin(M - Mpr)
378 + 0.0004 * sin(2*F+M)
379 - 0.0004 * sin(2*F-M)
380 - 0.0006 * sin(2*F+Mpr)
381 + 0.0010 * sin(2*F-Mpr)
382 + 0.0005 * sin(M+2*Mpr);
383 jd = jd + cor;
384 }
385 else {
386 cor = (0.1721 - 0.0004*T) * sin(M)
387 + 0.0021 * sin(2 * M)
388 - 0.6280 * sin(Mpr)
389 + 0.0089 * sin(2 * Mpr)
390 - 0.0004 * sin(3 * Mpr)
391 + 0.0079 * sin(2*F)
392 - 0.0119 * sin(M + Mpr)
393 - 0.0047 * sin(M - Mpr)
394 + 0.0003 * sin(2 * F + M)
395 - 0.0004 * sin(2 * F - M)
396 - 0.0006 * sin(2 * F + Mpr)
397 + 0.0021 * sin(2 * F - Mpr)
398 + 0.0003 * sin(M + 2 * Mpr)
399 + 0.0004 * sin(M - 2 * Mpr)
400 - 0.0003 * sin(2*M + Mpr);
401 if(nph == 1) cor = cor + 0.0028 -
402 0.0004 * cos(M) + 0.0003 * cos(Mpr);
403 if(nph == 3) cor = cor - 0.0028 +
404 0.0004 * cos(M) - 0.0003 * cos(Mpr);
405 jd = jd + cor;
406
407 }
408 *jdout = jd;
409 }
410
411
412 // Added 2003-02-04 from Skycal 5 for moonrise/moonset
circulo(double x)413 static double circulo (double x) {
414 /* assuming x is an angle in degrees, returns
415 modulo 360 degrees. */
416
417 int n;
418
419 n = (int)(x / 360.);
420 return(x - 360. * n);
421 }
422
423
424 // Added 2003-02-04 from Skycal 5 for moonrise/moonset
geocent(double geolong,double geolat,double height,double * x_geo,double * y_geo,double * z_geo)425 static void geocent (double geolong, double geolat, double height,
426 double *x_geo, double *y_geo, double *z_geo)
427
428 /* computes the geocentric coordinates from the geodetic
429 (standard map-type) longitude, latitude, and height.
430 These are assumed to be in decimal hours, decimal degrees, and
431 meters respectively. Notation generally follows 1992 Astr Almanac,
432 p. K11 */
433
434 {
435
436 double denom, C_geo, S_geo;
437
438 geolat = geolat / DEG_IN_RADIAN;
439 geolong = geolong / HRS_IN_RADIAN;
440 denom = (1. - FLATTEN) * sin(geolat);
441 denom = cos(geolat) * cos(geolat) + denom*denom;
442 C_geo = 1. / sqrt(denom);
443 S_geo = (1. - FLATTEN) * (1. - FLATTEN) * C_geo;
444 C_geo = C_geo + height / EQUAT_RAD; /* deviation from almanac
445 notation -- include height here. */
446 S_geo = S_geo + height / EQUAT_RAD;
447 *x_geo = C_geo * cos(geolat) * cos(geolong);
448 *y_geo = C_geo * cos(geolat) * sin(geolong);
449 *z_geo = S_geo * sin(geolat);
450 }
451
452
453 // Added 2003-02-04 from Skycal 5 for moonrise/moonset
eclrot(double jd,double * x unusedParameter,double * y,double * z)454 static void eclrot (double jd, double *x unusedParameter, double *y, double *z)
455
456 /* rotates ecliptic rectangular coords x, y, z to
457 equatorial (all assumed of date.) */
458
459 {
460 double incl;
461 double ypr,zpr;
462 double T;
463
464 T = (jd - J2000) / 36525; /* centuries since J2000 */
465
466 incl = (23.439291 + T * (-0.0130042 - 0.00000016 * T))/DEG_IN_RADIAN;
467 /* 1992 Astron Almanac, p. B18, dropping the
468 cubic term, which is 2 milli-arcsec! */
469 ypr = cos(incl) * *y - sin(incl) * *z;
470 zpr = sin(incl) * *y + cos(incl) * *z;
471 *y = ypr;
472 *z = zpr;
473 /* x remains the same. */
474 }
475
476
477 // Code by James.Ashton@anu.edu.au 2013-05-24
478 // Use the polynomial approximations from
479 // http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
480 // good for 1999 BC to 3000 AD, derived by NASA from
481 // Morrison, L. and Stephenson, F. R., "Historical Values of the Earth's
482 // Clock Error delt and the Calculation of Eclipses", J. Hist. Astron.,
483 // Vol. 35 Part 3, August 2004, No. 120, pp 327-336 (2004).
etcorr(double jd)484 static double etcorr (double jd) {
485 double y = (jd - J2000) / GREG_DAYS_IN_YEAR + 2000;
486 double t, delt;
487
488 if (y < -500) {
489 t = (y-1820)/100;
490 delt = -20 + 32 * t * t;
491 } else if (y < 500) {
492 t = y/100;
493 delt = 10583.6 - 1014.41 * t + 33.78311 * t * t - 5.952053 * t * t * t - 0.1798452 * t * t * t * t + 0.022174192 * t * t * t * t * t + 0.0090316521 * t * t * t * t * t * t;
494 } else if (y < 1600) {
495 t = (y-1000)/100;
496 delt = 1574.2 - 556.01 * t + 71.23472 * t * t + 0.319781 * t * t * t - 0.8503463 * t * t * t * t - 0.005050998 * t * t * t * t * t + 0.0083572073 * t * t * t * t * t * t;
497 } else if (y < 1700) {
498 t = y - 1600;
499 delt = 120 - 0.9808 * t - 0.01532 * t * t + t * t * t / 7129;
500 } else if (y < 1800) {
501 t = y - 1700;
502 delt = 8.83 + 0.1603 * t - 0.0059285 * t * t + 0.00013336 * t * t * t - t * t * t * t / 1174000;
503 } else if (y < 1860) {
504 t = y - 1800;
505 delt = 13.72 - 0.332447 * t + 0.0068612 * t * t + 0.0041116 * t * t * t - 0.00037436 * t * t * t * t + 0.0000121272 * t * t * t * t * t - 0.0000001699 * t * t * t * t * t * t + 0.000000000875 * t * t * t * t * t * t * t;
506 } else if (y < 1900) {
507 t = y - 1860;
508 delt = 7.62 + 0.5737 * t - 0.251754 * t * t + 0.01680668 * t * t * t -0.0004473624 * t * t * t * t + t * t * t * t * t / 233174;
509 } else if (y < 1920) {
510 t = y - 1900;
511 delt = -2.79 + 1.494119 * t - 0.0598939 * t * t + 0.0061966 * t * t * t - 0.000197 * t * t * t * t;
512 } else if (y < 1941) {
513 t = y - 1920;
514 delt = 21.20 + 0.84493 * t - 0.076100 * t * t + 0.0020936 * t * t * t;
515 } else if (y < 1961) {
516 t = y - 1950;
517 delt = 29.07 + 0.407 * t - t * t/233 + t * t * t / 2547;
518 } else if (y < 1986) {
519 t = y - 1975;
520 delt = 45.45 + 1.067 * t - t * t/260 - t * t * t / 718;
521 } else if (y < 2005) {
522 t = y - 2000;
523 delt = 63.86 + 0.3345 * t - 0.060374 * t * t + 0.0017275 * t * t * t + 0.000651814 * t * t * t * t + 0.00002373599 * t * t * t * t * t;
524 } else if (y < 2050) {
525 t = y - 2000;
526 delt = 62.92 + 0.32217 * t + 0.005589 * t * t;
527 } else if (y < 2150) {
528 t = (y - 1820)/100;
529 delt = -20 + 32 * t * t - 0.5628 * (2150 - y);
530 } else {
531 t = (y - 1820) / 100;
532 delt = -20 + 32 * t * t;
533 }
534
535 return(jd + delt/SEC_IN_DAY);
536 }
537
538
539 // Added 2003-02-04 from Skycal 5 for moonrise/moonset
accumoon(double jd,double geolat,double lst,double elevsea,double * topora,double * topodec,double * topodist)540 static void accumoon (double jd, double geolat, double lst, double elevsea,
541 double *topora, double *topodec, double *topodist)
542
543 // double jd,geolat,lst,elevsea; /* jd, dec. degr., dec. hrs., meters */
544
545 /* More accurate (but more elaborate and slower) lunar
546 ephemeris, from Jean Meeus' *Astronomical Formulae For Calculators*,
547 pub. Willman-Bell. Includes all the terms given there. */
548
549 {
550 /* double *eclatit,*eclongit, *pie,*ra,*dec,*dist; geocent quantities,
551 formerly handed out but not in this version */
552 double pie, dist; /* horiz parallax */
553 double Lpr,M,Mpr,D,F,Om,T,Tsq,Tcb;
554 double e,lambda,B,beta,om1,om2;
555 double sinx, x, y, z, l, m, n;
556 double x_geo, y_geo, z_geo; /* geocentric position of *observer* */
557
558 jd = etcorr(jd); /* approximate correction to ephemeris time */
559 T = (jd - 2415020.) / 36525.; /* this based around 1900 ... */
560 Tsq = T * T;
561 Tcb = Tsq * T;
562
563 Lpr = 270.434164 + 481267.8831 * T - 0.001133 * Tsq
564 + 0.0000019 * Tcb;
565 M = 358.475833 + 35999.0498*T - 0.000150*Tsq
566 - 0.0000033*Tcb;
567 Mpr = 296.104608 + 477198.8491*T + 0.009192*Tsq
568 + 0.0000144*Tcb;
569 D = 350.737486 + 445267.1142*T - 0.001436 * Tsq
570 + 0.0000019*Tcb;
571 F = 11.250889 + 483202.0251*T -0.003211 * Tsq
572 - 0.0000003*Tcb;
573 Om = 259.183275 - 1934.1420*T + 0.002078*Tsq
574 + 0.0000022*Tcb;
575
576 Lpr = circulo(Lpr);
577 Mpr = circulo(Mpr);
578 M = circulo(M);
579 D = circulo(D);
580 F = circulo(F);
581 Om = circulo(Om);
582
583
584 sinx = sin((51.2 + 20.2 * T)/DEG_IN_RADIAN);
585 Lpr = Lpr + 0.000233 * sinx;
586 M = M - 0.001778 * sinx;
587 Mpr = Mpr + 0.000817 * sinx;
588 D = D + 0.002011 * sinx;
589
590 sinx = 0.003964 * sin((346.560+132.870*T -0.0091731*Tsq)/DEG_IN_RADIAN);
591
592 Lpr = Lpr + sinx;
593 Mpr = Mpr + sinx;
594 D = D + sinx;
595 F = F + sinx;
596
597 sinx = sin(Om/DEG_IN_RADIAN);
598 Lpr = Lpr + 0.001964 * sinx;
599 Mpr = Mpr + 0.002541 * sinx;
600 D = D + 0.001964 * sinx;
601 F = F - 0.024691 * sinx;
602 F = F - 0.004328 * sin((Om + 275.05 -2.30*T)/DEG_IN_RADIAN);
603
604 e = 1 - 0.002495 * T - 0.00000752 * Tsq;
605
606 M = M / DEG_IN_RADIAN; /* these will all be arguments ... */
607 Mpr = Mpr / DEG_IN_RADIAN;
608 D = D / DEG_IN_RADIAN;
609 F = F / DEG_IN_RADIAN;
610
611 lambda = Lpr + 6.288750 * sin(Mpr)
612 + 1.274018 * sin(2*D - Mpr)
613 + 0.658309 * sin(2*D)
614 + 0.213616 * sin(2*Mpr)
615 - e * 0.185596 * sin(M)
616 - 0.114336 * sin(2*F)
617 + 0.058793 * sin(2*D - 2*Mpr)
618 + e * 0.057212 * sin(2*D - M - Mpr)
619 + 0.053320 * sin(2*D + Mpr)
620 + e * 0.045874 * sin(2*D - M)
621 + e * 0.041024 * sin(Mpr - M)
622 - 0.034718 * sin(D)
623 - e * 0.030465 * sin(M+Mpr)
624 + 0.015326 * sin(2*D - 2*F)
625 - 0.012528 * sin(2*F + Mpr)
626 - 0.010980 * sin(2*F - Mpr)
627 + 0.010674 * sin(4*D - Mpr)
628 + 0.010034 * sin(3*Mpr)
629 + 0.008548 * sin(4*D - 2*Mpr)
630 - e * 0.007910 * sin(M - Mpr + 2*D)
631 - e * 0.006783 * sin(2*D + M)
632 + 0.005162 * sin(Mpr - D);
633
634 /* And furthermore.....*/
635
636 lambda = lambda + e * 0.005000 * sin(M + D)
637 + e * 0.004049 * sin(Mpr - M + 2*D)
638 + 0.003996 * sin(2*Mpr + 2*D)
639 + 0.003862 * sin(4*D)
640 + 0.003665 * sin(2*D - 3*Mpr)
641 + e * 0.002695 * sin(2*Mpr - M)
642 + 0.002602 * sin(Mpr - 2*F - 2*D)
643 + e * 0.002396 * sin(2*D - M - 2*Mpr)
644 - 0.002349 * sin(Mpr + D)
645 + e * e * 0.002249 * sin(2*D - 2*M)
646 - e * 0.002125 * sin(2*Mpr + M)
647 - e * e * 0.002079 * sin(2*M)
648 + e * e * 0.002059 * sin(2*D - Mpr - 2*M)
649 - 0.001773 * sin(Mpr + 2*D - 2*F)
650 - 0.001595 * sin(2*F + 2*D)
651 + e * 0.001220 * sin(4*D - M - Mpr)
652 - 0.001110 * sin(2*Mpr + 2*F)
653 + 0.000892 * sin(Mpr - 3*D)
654 - e * 0.000811 * sin(M + Mpr + 2*D)
655 + e * 0.000761 * sin(4*D - M - 2*Mpr)
656 + e * e * 0.000717 * sin(Mpr - 2*M)
657 + e * e * 0.000704 * sin(Mpr - 2 * M - 2*D)
658 + e * 0.000693 * sin(M - 2*Mpr + 2*D)
659 + e * 0.000598 * sin(2*D - M - 2*F)
660 + 0.000550 * sin(Mpr + 4*D)
661 + 0.000538 * sin(4*Mpr)
662 + e * 0.000521 * sin(4*D - M)
663 + 0.000486 * sin(2*Mpr - D);
664
665 /* *eclongit = lambda; */
666
667 B = 5.128189 * sin(F)
668 + 0.280606 * sin(Mpr + F)
669 + 0.277693 * sin(Mpr - F)
670 + 0.173238 * sin(2*D - F)
671 + 0.055413 * sin(2*D + F - Mpr)
672 + 0.046272 * sin(2*D - F - Mpr)
673 + 0.032573 * sin(2*D + F)
674 + 0.017198 * sin(2*Mpr + F)
675 + 0.009267 * sin(2*D + Mpr - F)
676 + 0.008823 * sin(2*Mpr - F)
677 + e * 0.008247 * sin(2*D - M - F)
678 + 0.004323 * sin(2*D - F - 2*Mpr)
679 + 0.004200 * sin(2*D + F + Mpr)
680 + e * 0.003372 * sin(F - M - 2*D)
681 + 0.002472 * sin(2*D + F - M - Mpr)
682 + e * 0.002222 * sin(2*D + F - M)
683 + e * 0.002072 * sin(2*D - F - M - Mpr)
684 + e * 0.001877 * sin(F - M + Mpr)
685 + 0.001828 * sin(4*D - F - Mpr)
686 - e * 0.001803 * sin(F + M)
687 - 0.001750 * sin(3*F)
688 + e * 0.001570 * sin(Mpr - M - F)
689 - 0.001487 * sin(F + D)
690 - e * 0.001481 * sin(F + M + Mpr)
691 + e * 0.001417 * sin(F - M - Mpr)
692 + e * 0.001350 * sin(F - M)
693 + 0.001330 * sin(F - D)
694 + 0.001106 * sin(F + 3*Mpr)
695 + 0.001020 * sin(4*D - F)
696 + 0.000833 * sin(F + 4*D - Mpr);
697 /* not only that, but */
698 B = B + 0.000781 * sin(Mpr - 3*F)
699 + 0.000670 * sin(F + 4*D - 2*Mpr)
700 + 0.000606 * sin(2*D - 3*F)
701 + 0.000597 * sin(2*D + 2*Mpr - F)
702 + e * 0.000492 * sin(2*D + Mpr - M - F)
703 + 0.000450 * sin(2*Mpr - F - 2*D)
704 + 0.000439 * sin(3*Mpr - F)
705 + 0.000423 * sin(F + 2*D + 2*Mpr)
706 + 0.000422 * sin(2*D - F - 3*Mpr)
707 - e * 0.000367 * sin(M + F + 2*D - Mpr)
708 - e * 0.000353 * sin(M + F + 2*D)
709 + 0.000331 * sin(F + 4*D)
710 + e * 0.000317 * sin(2*D + F - M + Mpr)
711 + e * e * 0.000306 * sin(2*D - 2*M - F)
712 - 0.000283 * sin(Mpr + 3*F);
713
714 om1 = 0.0004664 * cos(Om/DEG_IN_RADIAN);
715 om2 = 0.0000754 * cos((Om + 275.05 - 2.30*T)/DEG_IN_RADIAN);
716
717 beta = B * (1. - om1 - om2);
718 /* *eclatit = beta; */
719
720 pie = 0.950724
721 + 0.051818 * cos(Mpr)
722 + 0.009531 * cos(2*D - Mpr)
723 + 0.007843 * cos(2*D)
724 + 0.002824 * cos(2*Mpr)
725 + 0.000857 * cos(2*D + Mpr)
726 + e * 0.000533 * cos(2*D - M)
727 + e * 0.000401 * cos(2*D - M - Mpr)
728 + e * 0.000320 * cos(Mpr - M)
729 - 0.000271 * cos(D)
730 - e * 0.000264 * cos(M + Mpr)
731 - 0.000198 * cos(2*F - Mpr)
732 + 0.000173 * cos(3*Mpr)
733 + 0.000167 * cos(4*D - Mpr)
734 - e * 0.000111 * cos(M)
735 + 0.000103 * cos(4*D - 2*Mpr)
736 - 0.000084 * cos(2*Mpr - 2*D)
737 - e * 0.000083 * cos(2*D + M)
738 + 0.000079 * cos(2*D + 2*Mpr)
739 + 0.000072 * cos(4*D)
740 + e * 0.000064 * cos(2*D - M + Mpr)
741 - e * 0.000063 * cos(2*D + M - Mpr)
742 + e * 0.000041 * cos(M + D)
743 + e * 0.000035 * cos(2*Mpr - M)
744 - 0.000033 * cos(3*Mpr - 2*D)
745 - 0.000030 * cos(Mpr + D)
746 - 0.000029 * cos(2*F - 2*D)
747 - e * 0.000029 * cos(2*Mpr + M)
748 + e * e * 0.000026 * cos(2*D - 2*M)
749 - 0.000023 * cos(2*F - 2*D + Mpr)
750 + e * 0.000019 * cos(4*D - M - Mpr);
751
752 beta = beta/DEG_IN_RADIAN;
753 lambda = lambda/DEG_IN_RADIAN;
754 l = cos(lambda) * cos(beta);
755 m = sin(lambda) * cos(beta);
756 n = sin(beta);
757 eclrot(jd,&l,&m,&n);
758
759 dist = 1/sin((pie)/DEG_IN_RADIAN);
760 x = l * dist;
761 y = m * dist;
762 z = n * dist;
763
764 /* *ra = atan_circ(l,m) * DEG_IN_RADIAN;
765 *dec = asin(n) * DEG_IN_RADIAN; */
766
767 geocent(lst,geolat,elevsea,&x_geo,&y_geo,&z_geo);
768
769 x = x - x_geo; /* topocentric correction using elliptical earth fig. */
770 y = y - y_geo;
771 z = z - z_geo;
772
773 *topodist = sqrt(x*x + y*y + z*z);
774
775 l = x / (*topodist);
776 m = y / (*topodist);
777 n = z / (*topodist);
778
779 *topora = atan_circ(l,m) * HRS_IN_RADIAN;
780 *topodec = asin(n) * DEG_IN_RADIAN;
781 }
782
783
784 // This began as print_phase in skycalc.c.
785 // No important changes in Skycal V5.
786 static void
find_next_moon_phase(double & jd,int & phase)787 find_next_moon_phase (double &jd, int &phase) {
788 double newjd, lastnewjd, nextjd;
789 short kount=0;
790
791 // Originally, there was no problem with getting snagged, but since
792 // I introduced the roundoff error going back and forth with Timestamp,
793 // now it's a problem.
794 // Move ahead by 1 second to avoid snagging.
795 jd += 1.0 / SEC_IN_DAY;
796
797 // Find current lunation. I doubled the safety margin since it
798 // seemed biased for forwards search. Backwards search has since
799 // been deleted, but little reason to mess with it.
800 int nlast = (int)((jd - 2415020.5) / 29.5307 - 2);
801
802 flmoon(nlast,0,&lastnewjd);
803 flmoon(++nlast,0,&newjd);
804 while (newjd <= jd) {
805 lastnewjd = newjd;
806 flmoon(++nlast,0,&newjd);
807 require (kount++ < 5); // Original limit was 35 (!)
808 }
809
810 // We might save some work here by estimating, i.e.:
811 // x = jd - lastnewjd;
812 // noctiles = (int)(x / 3.69134); /* 3.69134 = 1/8 month; truncate. */
813 // However....
814
815 assert (lastnewjd <= jd && newjd > jd);
816 phase = 1;
817 // Lunation is lastnewjd's lunation
818 flmoon (--nlast, phase, &nextjd); // Phase = 1
819 if (nextjd <= jd) {
820 flmoon (nlast, ++phase, &nextjd); // Phase = 2
821 if (nextjd <= jd) {
822 flmoon (nlast, ++phase, &nextjd); // Phase = 3
823 if (nextjd <= jd) {
824 phase = 0;
825 nextjd = newjd;
826 }
827 }
828 }
829 jd = nextjd;
830 }
831
832
findNextMoonPhase(Timestamp t,TideEvent & tideEvent_out)833 void libxtide::Skycal::findNextMoonPhase (Timestamp t,
834 TideEvent &tideEvent_out) {
835 int phase;
836 double jd (t.jd());
837 find_next_moon_phase (jd, phase);
838 tideEvent_out.eventTime = jd;
839 switch (phase) {
840 case 0:
841 tideEvent_out.eventType = TideEvent::newmoon;
842 break;
843 case 1:
844 tideEvent_out.eventType = TideEvent::firstquarter;
845 break;
846 case 2:
847 tideEvent_out.eventType = TideEvent::fullmoon;
848 break;
849 case 3:
850 tideEvent_out.eventType = TideEvent::lastquarter;
851 break;
852 default:
853 assert (false);
854 }
855 }
856
857
858 // Here's another opportunity for Jeff Dairiki to write a better root
859 // finder :-)
860 //
861 // jd_sun_alt needed good initial guesses to find sunrises and
862 // sunsets. This was not a problem since good guesses were easy to
863 // come by. The original skycalendar did this with estimates based on
864 // the local midnight:
865 //
866 // jd = date_to_jd(date); /* local midnight */
867 // jdmid = jd + zone(use_dst,stdz,jd,jdbdst,jdedst) / 24.; /* corresponding ut */
868 // stmid = lst(jdmid,longit);
869 // lpsun(jdmid,&rasun,&decsun);
870 // hasunset = ha_alt(decsun,lat,-(0.83+horiz));
871 // jdsunset = jdmid + adj_time(rasun+hasunset-stmid)/24.; /* initial guess */
872 // jdsunset = jd_sun_alt(-(0.83+horiz), jdsunset,lat,longit); /* refinement */
873 // jdsunrise = jdmid + adj_time(rasun-hasunset-stmid)/24.;
874 // jdsunrise = jd_sun_alt(-(0.83+horiz),jdsunrise,lat,longit);
875 //
876 // While efficient, this is an inconvenient way to go about it when
877 // I'm looking for the next event from time t, and don't even know
878 // when midnight is. So I messed with jd_sun_alt to make it converge
879 // better from bad initial guesses, and substituted three bad guesses
880 // for one good one. Normally, two would suffice, but I wanted to
881 // add a safety margin in case one of them happens to land at a point
882 // that nukes jd_sun_alt.
883
884 // 2003-02-04
885 // Expanded to handle moonrise/moonset as well.
886
887 // Set lunar to true for moonrise/set.
888 static void
find_next_rise_or_set(double & jd,double lat,double longit,bool lunar,bool & is_rise)889 find_next_rise_or_set (double &jd, double lat, double longit, bool lunar,
890 bool &is_rise) {
891 // Move ahead by precision interval to avoid snagging.
892 jd += libxtide::Global::eventPrecisionJD;
893
894 double jdorig = jd;
895 double inc = 1.0 / 6.0; // 4 hours
896
897 // First we want to know what we are looking for.
898 bool looking_for = (altitude (jdorig, lat, longit, lunar) < riseAltitude);
899
900 // Now give it a decent try. Because jd_alt is so unpredictable,
901 // we can even find things out of order (which is one reason we need
902 // to know what we're looking for).
903 double jdlooper = jdorig;
904 do {
905 jd = jd_alt (riseAltitude, jdlooper, lat, longit, lunar, is_rise);
906 jdlooper += inc;
907 // Loop either on error return (which is a negative number), or if we
908 // found an event in the wrong direction, or the wrong kind of event.
909 } while ((jd < 0.0) ||
910 (jd <= jdorig) ||
911 (is_rise != looking_for));
912 }
913
914
findNextRiseOrSet(Timestamp t,const Coordinates & c,RiseSetType riseSetType,TideEvent & tideEvent_out)915 void libxtide::Skycal::findNextRiseOrSet (Timestamp t,
916 const Coordinates &c,
917 RiseSetType riseSetType,
918 TideEvent &tideEvent_out) {
919 assert (!(c.isNull()));
920 bool isRise;
921 double jd = t.jd();
922 // skycal "longit" is measured in HOURS WEST, not degrees east.
923 // (lat is unchanged)
924 find_next_rise_or_set (jd,
925 c.lat(),
926 -(c.lng())/15.0,
927 (riseSetType==lunar),
928 isRise);
929 tideEvent_out.eventTime = jd;
930 if (isRise)
931 tideEvent_out.eventType = ((riseSetType == lunar) ? TideEvent::moonrise
932 : TideEvent::sunrise);
933 else
934 tideEvent_out.eventType = ((riseSetType == lunar) ? TideEvent::moonset
935 : TideEvent::sunset);
936 }
937
938
939 // Simple question deserving a simple answer...
sunIsUp(Timestamp t,const Coordinates & c)940 const bool libxtide::Skycal::sunIsUp (Timestamp t, const Coordinates &c) {
941 assert (!(c.isNull()));
942 return (altitude (t.jd(), c.lat(), -(c.lng())/15.0, 0) >= riseAltitude);
943 }
944
945
946 #ifdef EXPERIMENTAL_MOON_AGE_NOT_PHASE
findNewMoons(Timestamp t,Timestamp & prev_out,Timestamp & next_out)947 void libxtide::Skycal::findNewMoons (Timestamp t,
948 Timestamp &prev_out,
949 Timestamp &next_out) {
950 double jd(t.jd()), newjd, lastnewjd;
951 short kount=0;
952 // Duplicated from find_next_moon_phase "find current lunation" block.
953 // IDK whether 29.5307 really wants to be 29.530588853?
954 int nlast = (int)((jd - 2415020.5) / 29.5307 - 2);
955 flmoon(nlast,0,&lastnewjd);
956 flmoon(++nlast,0,&newjd);
957 while (newjd <= jd) {
958 lastnewjd = newjd;
959 flmoon(++nlast,0,&newjd);
960 require (kount++ < 5); // Original limit was 35 (!)
961 }
962 // And presto, we're done.
963 prev_out = lastnewjd;
964 next_out = newjd;
965 }
966 #endif
967