// $Id: Skycal.cc 5748 2014-10-11 19:38:53Z flaterco $ // Skycal.cc -- Functions for sun and moon events. // Prediction of moon phases, sun and moon rises and sets has nothing // to do with tide prediction. There is no overlap between this code // and the tide prediction code. // This source file began its life as skycalendar.c and skycalc.c in // John Thorstensen's skycal distribution (version 4.1, 1994-09) at // ftp://iraf.noao.edu/contrib/skycal.tar.Z. Those portions that are // unchanged from the original sources are covered by the original // license statement, included below. The new portions and "value // added" by David Flater are covered under the GNU General Public // License. // 2013-05-26 // // Patch from James Ashton to replace etcorr. // 2003-02-04 // // The release notes for Skycal V5 indicated that some code in 4.1 had // copyright problems. Reviewed all code snippets here and found that // none were impacted by copyright-related changes. // // Harmonized atan_circ with slightly improved V5 version. /* Copyright (C) 1998 David Flater. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ // The relevant portions of the Skycal 5 comments and license statement // are as follows: /* SKY CALCULATOR PROGRAM John Thorstensen, Dartmouth College. This program computes many quantities frequently needed by the observational astronomer. It is written as a completely self-contained program in standard c, so it should be very transportable; the only issue I know of that really affects portability is the adequacy of the double-precision floating point accuracy on the machine. Experience shows that c compilers on various systems have idiosyncracies, though, so be sure to check carefully. This is intended as an observatory utility program; I assume the user is familiar with astronomical coordinates and nomenclature. While the code should be very transportable, I also assume it will be installed by a conscientious person who will run critical tests before it is released at a new site. Experience shows that some c compilers generate unforseen errors when the code is ported, so the output should be checked meticulously against data from other sites. [...] The program is self-contained. It has been developed primarily on UNIX and Linux machines, and should adapt easily to any system with a c compiler. ** BUT CAUTION ... ** Because many of the routines take a double-precision floating point Julian Date as their time argument, one must be sure that the machine and compiler carry sufficient mantissa to reach the desired accuracy. On most architectures the double-precision floating point julian date has an accuracy of order 0.01 seconds of time, which is just adequate. LEGALITIES: I make no guarantee as to the accuracy, reliability, or appropriateness of this program, though I have found it to be reasonably accurate and quite useful to the working astronomer. The program is COPYRIGHT 2000 BY JOHN THORSTENSEN. Permission is hereby granted for non-profit scientific or educational use. For-profit use (e. g., by astrologers!) must be through negotiated license. The author requests that observatories and astronomy departments which install this as a utility notify the author by paper mail, just so I know how widely it is used. Credits: * The julian date and sidereal time routines were originally coded in PL/I by Steve Maker of Dartmouth College. They were based on routines in the old American Ephemeris. Many of the routines were coded from Jean Meeus' "Astronomical Formulae for Calculators", published by Willman-Bell. This is an extraordinarily helpful little book! */ // The Skycal 4.1 comments and license statement are as follows: /* This is a self-contained c-language program to print a nighttime astronomical calendar for use in planning observations. It prints to standard output (usually the terminal); the operator should capture this output (e. g., using redirection in UNIX or the /out= switch in VMS) and then print it on an appropriate output device. The table which is printed is some 125 columns wide, so a wide device is required (either a line printer or a laserprinter in LANDSCAPE mode.) It is assumed that the ASCII form-feed character will actually begin a new page. The original program was to run on VMS, but it should be very transportable. Non-vms users will probably want to change 'unixio.h' to 'stdio.h' in the first line. An explanatory text is printed at the beginning of the output, which includes the appropriate CAUTIONS regarding accuracy and applicability. A number of 'canned site' parameters have been included. Be careful of time zones, DST etc. for foreign sites. To customize to your own site, install an option in the routine 'load_site'. The code is very straightforward; just do it exactly the same as the others. You might also want to erase some seldom-used choices. One can also specify new site parameters at run time. This program may be used freely by anyone for scientific or educational purposes. If you use it for profit, I want a cut, and claim a copyright herewith. In any case please acknowledge the source: John Thorstensen Dept. of Physics and Astronomy Dartmouth College Hanover, NH 03755 John.Thorstensen@dartmouth.edu May 26, 1993. */ #include "libxtide.hh" #include "Skycal.hh" // To avoid creating a clash on every reference of a math.h trig function, // namespace libxtide is applied only to the non-static definitions. // DWF: This affects the rise/set predictions. Normally you would // need to adjust it for the elevation of the location, but since this // is a tide prediction program, we can safely assume that we are // always at sea level :-) static const double riseAltitude (-0.83); #define DEG_IN_RADIAN 57.2957795130823 #define HRS_IN_RADIAN 3.819718634 #define SEC_IN_DAY 86400. #define GREG_DAYS_IN_YEAR 365.2425 /* 365 + 97 / 400 */ #define FLATTEN 0.003352813 /* flattening of earth, 1/298.257 */ #define EQUAT_RAD 6378137. /* equatorial radius of earth, meters */ #define J2000 2451545. /* Julian date at standard epoch */ // Harmonized with Skycal V5 2003-02-04 static double atan_circ(double x, double y) { /* returns radian angle 0 to 2pi for coords x, y -- get that quadrant right !! */ double theta; if((x == 0.) && (y == 0.)) return(0.); /* guard ... */ theta = atan2(y,x); /* turns out there is such a thing in math.h */ while(theta < 0.) theta += 2.0 * M_PI; return(theta); } // This is stripped down from Skycal V4.1 and even more stripped down // from Skycal V5. static double altit (double dec, double ha, double lat) /* dec deg, dec hrs, dec deg */ { double x; dec = dec / DEG_IN_RADIAN; ha = ha / HRS_IN_RADIAN; lat = lat / DEG_IN_RADIAN; /* thank heavens for pass-by-value */ x = DEG_IN_RADIAN * asin(cos(dec)*cos(ha)*cos(lat) + sin(dec)*sin(lat)); return(x); } // No important changes in Skycal V5. static double lst(double jd, double longit) { /* returns the local MEAN sidereal time (dec hrs) at julian date jd at west longitude long (decimal hours). Follows definitions in 1992 Astronomical Almanac, pp. B7 and L2. Expression for GMST at 0h ut referenced to Aoki et al, A&A 105, p.359, 1982. */ double t, ut, jdmid, jdint, jdfrac, sid_g; long jdin, sid_int; jdin = (long)jd; /* fossil code from earlier package which split jd into integer and fractional parts ... */ jdint = jdin; jdfrac = jd - jdint; if(jdfrac < 0.5) { jdmid = jdint - 0.5; ut = jdfrac + 0.5; } else { jdmid = jdint + 0.5; ut = jdfrac - 0.5; } t = (jdmid - J2000)/36525; sid_g = (24110.54841+8640184.812866*t+0.093104*t*t-6.2e-6*t*t*t)/86400.; sid_int = (long)sid_g; sid_g = sid_g - (double) sid_int; sid_g = sid_g + 1.0027379093 * ut - longit/24.; sid_int = (long)sid_g; sid_g = (sid_g - (double) sid_int) * 24.; if(sid_g < 0.) sid_g = sid_g + 24.; return(sid_g); } // No important changes in Skycal V5. static void lpsun(double jd, double *ra, double *dec) /* Low precision formulae for the sun, from Almanac p. C24 (1990) */ /* ra and dec are returned as decimal hours and decimal degrees. */ { double n, L, g, lambda,epsilon,x,y,z; n = jd - J2000; L = 280.460 + 0.9856474 * n; g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN; lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN; epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN; x = cos(lambda); y = cos(epsilon) * sin(lambda); z = sin(epsilon)*sin(lambda); *ra = (atan_circ(x,y))*HRS_IN_RADIAN; *dec = (asin(z))*DEG_IN_RADIAN; } // Fwd Decl static void accumoon (double jd, double geolat, double lst, double elevsea, double *topora, double *topodec, double *topodist); // Not in Skycal. // This function combines a few steps that are always used together to // find the altitude of the sun or moon. static double altitude (double jd, double lat, double longit, bool lunar) { if (lunar) { double ra, dec, dist, sid; sid = lst(jd,longit); accumoon (jd, lat, sid, 0, &ra, &dec, &dist); return altit(dec, sid-ra, lat); } else { double ra, dec; lpsun(jd, &ra, &dec); return altit(dec, lst(jd,longit)-ra, lat); } } // I messed with this a good bit. // // * It now converges to within 1 minute. // * It converges better from bad initial guesses. (Deriv is now // updated inside of the loop.) // * It won't roam more than half a day in either direction. // * Max iterations chosen conservatively. // * It finishes with a check to determine what it found. // 2003-02-04 // Expanded to handle moonrise/moonset (replacing jd_moon_alt too). // No important changes in Skycal V5. static double jd_alt (double alt, double jdorig, double lat, double longit, bool lunar, bool &is_rise) { /* returns jd at which sun/moon is at a given altitude, given jdguess as a starting point. */ double jdguess = jdorig; double jdout, adj = 1.0; double deriv, err, del = 0.002; double alt2,alt3; short i = 0; /* first guess */ alt2 = altitude (jdguess, lat, longit, lunar); jdguess = jdguess + del; alt3 = altitude (jdguess, lat, longit, lunar); err = alt3 - alt; deriv = (alt3 - alt2) / del; if (deriv == 0.0) return (-1.0e10); // Found dead end. adj = -err/deriv; while(fabs(adj) >= libxtide::Global::eventPrecisionJD) { if (i++ == 12) return(-1.0e10); // Exceeded max iterations. jdguess += adj; if (fabs (jdguess - jdorig) > 0.5) return (-1.0e10); // Ran out of bounds. alt2 = alt3; alt3 = altitude (jdguess, lat, longit, lunar); err = alt3 - alt; deriv = (alt3 - alt2) / adj; if (deriv == 0.0) return (-1.0e10); // Found dead end. adj = -err/deriv; } jdout = jdguess; // Figure out whether this is a rise or a set by shifting // by 1 second. { jdguess -= 1.0 / SEC_IN_DAY; alt2 = altitude (jdguess, lat, longit, lunar); is_rise = (alt2 < alt3); } return(jdout); } // No important changes in Skycal V5. static void flmoon(int n, int nph, double *jdout) /* Gives jd (+- 2 min) of phase nph on lunation n; replaces less accurate Numerical Recipes routine. This routine implements formulae found in Jean Meeus' *Astronomical Formulae for Calculators*, 2nd edition, Willman-Bell. A very useful book!! */ { double jd, cor; double M, Mpr, F; double T; double lun; lun = (double) n + (double) nph / 4.; T = lun / 1236.85; jd = 2415020.75933 + 29.53058868 * lun + 0.0001178 * T * T - 0.000000155 * T * T * T + 0.00033 * sin((166.56 + 132.87 * T - 0.009173 * T * T)/DEG_IN_RADIAN); M = 359.2242 + 29.10535608 * lun - 0.0000333 * T * T - 0.00000347 * T * T * T; M = M / DEG_IN_RADIAN; Mpr = 306.0253 + 385.81691806 * lun + 0.0107306 * T * T + 0.00001236 * T * T * T; Mpr = Mpr / DEG_IN_RADIAN; F = 21.2964 + 390.67050646 * lun - 0.0016528 * T * T - 0.00000239 * T * T * T; F = F / DEG_IN_RADIAN; if((nph == 0) || (nph == 2)) {/* new or full */ cor = (0.1734 - 0.000393*T) * sin(M) + 0.0021 * sin(2*M) - 0.4068 * sin(Mpr) + 0.0161 * sin(2*Mpr) - 0.0004 * sin(3*Mpr) + 0.0104 * sin(2*F) - 0.0051 * sin(M + Mpr) - 0.0074 * sin(M - Mpr) + 0.0004 * sin(2*F+M) - 0.0004 * sin(2*F-M) - 0.0006 * sin(2*F+Mpr) + 0.0010 * sin(2*F-Mpr) + 0.0005 * sin(M+2*Mpr); jd = jd + cor; } else { cor = (0.1721 - 0.0004*T) * sin(M) + 0.0021 * sin(2 * M) - 0.6280 * sin(Mpr) + 0.0089 * sin(2 * Mpr) - 0.0004 * sin(3 * Mpr) + 0.0079 * sin(2*F) - 0.0119 * sin(M + Mpr) - 0.0047 * sin(M - Mpr) + 0.0003 * sin(2 * F + M) - 0.0004 * sin(2 * F - M) - 0.0006 * sin(2 * F + Mpr) + 0.0021 * sin(2 * F - Mpr) + 0.0003 * sin(M + 2 * Mpr) + 0.0004 * sin(M - 2 * Mpr) - 0.0003 * sin(2*M + Mpr); if(nph == 1) cor = cor + 0.0028 - 0.0004 * cos(M) + 0.0003 * cos(Mpr); if(nph == 3) cor = cor - 0.0028 + 0.0004 * cos(M) - 0.0003 * cos(Mpr); jd = jd + cor; } *jdout = jd; } // Added 2003-02-04 from Skycal 5 for moonrise/moonset static double circulo (double x) { /* assuming x is an angle in degrees, returns modulo 360 degrees. */ int n; n = (int)(x / 360.); return(x - 360. * n); } // Added 2003-02-04 from Skycal 5 for moonrise/moonset static void geocent (double geolong, double geolat, double height, double *x_geo, double *y_geo, double *z_geo) /* computes the geocentric coordinates from the geodetic (standard map-type) longitude, latitude, and height. These are assumed to be in decimal hours, decimal degrees, and meters respectively. Notation generally follows 1992 Astr Almanac, p. K11 */ { double denom, C_geo, S_geo; geolat = geolat / DEG_IN_RADIAN; geolong = geolong / HRS_IN_RADIAN; denom = (1. - FLATTEN) * sin(geolat); denom = cos(geolat) * cos(geolat) + denom*denom; C_geo = 1. / sqrt(denom); S_geo = (1. - FLATTEN) * (1. - FLATTEN) * C_geo; C_geo = C_geo + height / EQUAT_RAD; /* deviation from almanac notation -- include height here. */ S_geo = S_geo + height / EQUAT_RAD; *x_geo = C_geo * cos(geolat) * cos(geolong); *y_geo = C_geo * cos(geolat) * sin(geolong); *z_geo = S_geo * sin(geolat); } // Added 2003-02-04 from Skycal 5 for moonrise/moonset static void eclrot (double jd, double *x unusedParameter, double *y, double *z) /* rotates ecliptic rectangular coords x, y, z to equatorial (all assumed of date.) */ { double incl; double ypr,zpr; double T; T = (jd - J2000) / 36525; /* centuries since J2000 */ incl = (23.439291 + T * (-0.0130042 - 0.00000016 * T))/DEG_IN_RADIAN; /* 1992 Astron Almanac, p. B18, dropping the cubic term, which is 2 milli-arcsec! */ ypr = cos(incl) * *y - sin(incl) * *z; zpr = sin(incl) * *y + cos(incl) * *z; *y = ypr; *z = zpr; /* x remains the same. */ } // Code by James.Ashton@anu.edu.au 2013-05-24 // Use the polynomial approximations from // http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html // good for 1999 BC to 3000 AD, derived by NASA from // Morrison, L. and Stephenson, F. R., "Historical Values of the Earth's // Clock Error delt and the Calculation of Eclipses", J. Hist. Astron., // Vol. 35 Part 3, August 2004, No. 120, pp 327-336 (2004). static double etcorr (double jd) { double y = (jd - J2000) / GREG_DAYS_IN_YEAR + 2000; double t, delt; if (y < -500) { t = (y-1820)/100; delt = -20 + 32 * t * t; } else if (y < 500) { t = y/100; 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; } else if (y < 1600) { t = (y-1000)/100; 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; } else if (y < 1700) { t = y - 1600; delt = 120 - 0.9808 * t - 0.01532 * t * t + t * t * t / 7129; } else if (y < 1800) { t = y - 1700; delt = 8.83 + 0.1603 * t - 0.0059285 * t * t + 0.00013336 * t * t * t - t * t * t * t / 1174000; } else if (y < 1860) { t = y - 1800; 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; } else if (y < 1900) { t = y - 1860; 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; } else if (y < 1920) { t = y - 1900; delt = -2.79 + 1.494119 * t - 0.0598939 * t * t + 0.0061966 * t * t * t - 0.000197 * t * t * t * t; } else if (y < 1941) { t = y - 1920; delt = 21.20 + 0.84493 * t - 0.076100 * t * t + 0.0020936 * t * t * t; } else if (y < 1961) { t = y - 1950; delt = 29.07 + 0.407 * t - t * t/233 + t * t * t / 2547; } else if (y < 1986) { t = y - 1975; delt = 45.45 + 1.067 * t - t * t/260 - t * t * t / 718; } else if (y < 2005) { t = y - 2000; 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; } else if (y < 2050) { t = y - 2000; delt = 62.92 + 0.32217 * t + 0.005589 * t * t; } else if (y < 2150) { t = (y - 1820)/100; delt = -20 + 32 * t * t - 0.5628 * (2150 - y); } else { t = (y - 1820) / 100; delt = -20 + 32 * t * t; } return(jd + delt/SEC_IN_DAY); } // Added 2003-02-04 from Skycal 5 for moonrise/moonset static void accumoon (double jd, double geolat, double lst, double elevsea, double *topora, double *topodec, double *topodist) // double jd,geolat,lst,elevsea; /* jd, dec. degr., dec. hrs., meters */ /* More accurate (but more elaborate and slower) lunar ephemeris, from Jean Meeus' *Astronomical Formulae For Calculators*, pub. Willman-Bell. Includes all the terms given there. */ { /* double *eclatit,*eclongit, *pie,*ra,*dec,*dist; geocent quantities, formerly handed out but not in this version */ double pie, dist; /* horiz parallax */ double Lpr,M,Mpr,D,F,Om,T,Tsq,Tcb; double e,lambda,B,beta,om1,om2; double sinx, x, y, z, l, m, n; double x_geo, y_geo, z_geo; /* geocentric position of *observer* */ jd = etcorr(jd); /* approximate correction to ephemeris time */ T = (jd - 2415020.) / 36525.; /* this based around 1900 ... */ Tsq = T * T; Tcb = Tsq * T; Lpr = 270.434164 + 481267.8831 * T - 0.001133 * Tsq + 0.0000019 * Tcb; M = 358.475833 + 35999.0498*T - 0.000150*Tsq - 0.0000033*Tcb; Mpr = 296.104608 + 477198.8491*T + 0.009192*Tsq + 0.0000144*Tcb; D = 350.737486 + 445267.1142*T - 0.001436 * Tsq + 0.0000019*Tcb; F = 11.250889 + 483202.0251*T -0.003211 * Tsq - 0.0000003*Tcb; Om = 259.183275 - 1934.1420*T + 0.002078*Tsq + 0.0000022*Tcb; Lpr = circulo(Lpr); Mpr = circulo(Mpr); M = circulo(M); D = circulo(D); F = circulo(F); Om = circulo(Om); sinx = sin((51.2 + 20.2 * T)/DEG_IN_RADIAN); Lpr = Lpr + 0.000233 * sinx; M = M - 0.001778 * sinx; Mpr = Mpr + 0.000817 * sinx; D = D + 0.002011 * sinx; sinx = 0.003964 * sin((346.560+132.870*T -0.0091731*Tsq)/DEG_IN_RADIAN); Lpr = Lpr + sinx; Mpr = Mpr + sinx; D = D + sinx; F = F + sinx; sinx = sin(Om/DEG_IN_RADIAN); Lpr = Lpr + 0.001964 * sinx; Mpr = Mpr + 0.002541 * sinx; D = D + 0.001964 * sinx; F = F - 0.024691 * sinx; F = F - 0.004328 * sin((Om + 275.05 -2.30*T)/DEG_IN_RADIAN); e = 1 - 0.002495 * T - 0.00000752 * Tsq; M = M / DEG_IN_RADIAN; /* these will all be arguments ... */ Mpr = Mpr / DEG_IN_RADIAN; D = D / DEG_IN_RADIAN; F = F / DEG_IN_RADIAN; lambda = Lpr + 6.288750 * sin(Mpr) + 1.274018 * sin(2*D - Mpr) + 0.658309 * sin(2*D) + 0.213616 * sin(2*Mpr) - e * 0.185596 * sin(M) - 0.114336 * sin(2*F) + 0.058793 * sin(2*D - 2*Mpr) + e * 0.057212 * sin(2*D - M - Mpr) + 0.053320 * sin(2*D + Mpr) + e * 0.045874 * sin(2*D - M) + e * 0.041024 * sin(Mpr - M) - 0.034718 * sin(D) - e * 0.030465 * sin(M+Mpr) + 0.015326 * sin(2*D - 2*F) - 0.012528 * sin(2*F + Mpr) - 0.010980 * sin(2*F - Mpr) + 0.010674 * sin(4*D - Mpr) + 0.010034 * sin(3*Mpr) + 0.008548 * sin(4*D - 2*Mpr) - e * 0.007910 * sin(M - Mpr + 2*D) - e * 0.006783 * sin(2*D + M) + 0.005162 * sin(Mpr - D); /* And furthermore.....*/ lambda = lambda + e * 0.005000 * sin(M + D) + e * 0.004049 * sin(Mpr - M + 2*D) + 0.003996 * sin(2*Mpr + 2*D) + 0.003862 * sin(4*D) + 0.003665 * sin(2*D - 3*Mpr) + e * 0.002695 * sin(2*Mpr - M) + 0.002602 * sin(Mpr - 2*F - 2*D) + e * 0.002396 * sin(2*D - M - 2*Mpr) - 0.002349 * sin(Mpr + D) + e * e * 0.002249 * sin(2*D - 2*M) - e * 0.002125 * sin(2*Mpr + M) - e * e * 0.002079 * sin(2*M) + e * e * 0.002059 * sin(2*D - Mpr - 2*M) - 0.001773 * sin(Mpr + 2*D - 2*F) - 0.001595 * sin(2*F + 2*D) + e * 0.001220 * sin(4*D - M - Mpr) - 0.001110 * sin(2*Mpr + 2*F) + 0.000892 * sin(Mpr - 3*D) - e * 0.000811 * sin(M + Mpr + 2*D) + e * 0.000761 * sin(4*D - M - 2*Mpr) + e * e * 0.000717 * sin(Mpr - 2*M) + e * e * 0.000704 * sin(Mpr - 2 * M - 2*D) + e * 0.000693 * sin(M - 2*Mpr + 2*D) + e * 0.000598 * sin(2*D - M - 2*F) + 0.000550 * sin(Mpr + 4*D) + 0.000538 * sin(4*Mpr) + e * 0.000521 * sin(4*D - M) + 0.000486 * sin(2*Mpr - D); /* *eclongit = lambda; */ B = 5.128189 * sin(F) + 0.280606 * sin(Mpr + F) + 0.277693 * sin(Mpr - F) + 0.173238 * sin(2*D - F) + 0.055413 * sin(2*D + F - Mpr) + 0.046272 * sin(2*D - F - Mpr) + 0.032573 * sin(2*D + F) + 0.017198 * sin(2*Mpr + F) + 0.009267 * sin(2*D + Mpr - F) + 0.008823 * sin(2*Mpr - F) + e * 0.008247 * sin(2*D - M - F) + 0.004323 * sin(2*D - F - 2*Mpr) + 0.004200 * sin(2*D + F + Mpr) + e * 0.003372 * sin(F - M - 2*D) + 0.002472 * sin(2*D + F - M - Mpr) + e * 0.002222 * sin(2*D + F - M) + e * 0.002072 * sin(2*D - F - M - Mpr) + e * 0.001877 * sin(F - M + Mpr) + 0.001828 * sin(4*D - F - Mpr) - e * 0.001803 * sin(F + M) - 0.001750 * sin(3*F) + e * 0.001570 * sin(Mpr - M - F) - 0.001487 * sin(F + D) - e * 0.001481 * sin(F + M + Mpr) + e * 0.001417 * sin(F - M - Mpr) + e * 0.001350 * sin(F - M) + 0.001330 * sin(F - D) + 0.001106 * sin(F + 3*Mpr) + 0.001020 * sin(4*D - F) + 0.000833 * sin(F + 4*D - Mpr); /* not only that, but */ B = B + 0.000781 * sin(Mpr - 3*F) + 0.000670 * sin(F + 4*D - 2*Mpr) + 0.000606 * sin(2*D - 3*F) + 0.000597 * sin(2*D + 2*Mpr - F) + e * 0.000492 * sin(2*D + Mpr - M - F) + 0.000450 * sin(2*Mpr - F - 2*D) + 0.000439 * sin(3*Mpr - F) + 0.000423 * sin(F + 2*D + 2*Mpr) + 0.000422 * sin(2*D - F - 3*Mpr) - e * 0.000367 * sin(M + F + 2*D - Mpr) - e * 0.000353 * sin(M + F + 2*D) + 0.000331 * sin(F + 4*D) + e * 0.000317 * sin(2*D + F - M + Mpr) + e * e * 0.000306 * sin(2*D - 2*M - F) - 0.000283 * sin(Mpr + 3*F); om1 = 0.0004664 * cos(Om/DEG_IN_RADIAN); om2 = 0.0000754 * cos((Om + 275.05 - 2.30*T)/DEG_IN_RADIAN); beta = B * (1. - om1 - om2); /* *eclatit = beta; */ pie = 0.950724 + 0.051818 * cos(Mpr) + 0.009531 * cos(2*D - Mpr) + 0.007843 * cos(2*D) + 0.002824 * cos(2*Mpr) + 0.000857 * cos(2*D + Mpr) + e * 0.000533 * cos(2*D - M) + e * 0.000401 * cos(2*D - M - Mpr) + e * 0.000320 * cos(Mpr - M) - 0.000271 * cos(D) - e * 0.000264 * cos(M + Mpr) - 0.000198 * cos(2*F - Mpr) + 0.000173 * cos(3*Mpr) + 0.000167 * cos(4*D - Mpr) - e * 0.000111 * cos(M) + 0.000103 * cos(4*D - 2*Mpr) - 0.000084 * cos(2*Mpr - 2*D) - e * 0.000083 * cos(2*D + M) + 0.000079 * cos(2*D + 2*Mpr) + 0.000072 * cos(4*D) + e * 0.000064 * cos(2*D - M + Mpr) - e * 0.000063 * cos(2*D + M - Mpr) + e * 0.000041 * cos(M + D) + e * 0.000035 * cos(2*Mpr - M) - 0.000033 * cos(3*Mpr - 2*D) - 0.000030 * cos(Mpr + D) - 0.000029 * cos(2*F - 2*D) - e * 0.000029 * cos(2*Mpr + M) + e * e * 0.000026 * cos(2*D - 2*M) - 0.000023 * cos(2*F - 2*D + Mpr) + e * 0.000019 * cos(4*D - M - Mpr); beta = beta/DEG_IN_RADIAN; lambda = lambda/DEG_IN_RADIAN; l = cos(lambda) * cos(beta); m = sin(lambda) * cos(beta); n = sin(beta); eclrot(jd,&l,&m,&n); dist = 1/sin((pie)/DEG_IN_RADIAN); x = l * dist; y = m * dist; z = n * dist; /* *ra = atan_circ(l,m) * DEG_IN_RADIAN; *dec = asin(n) * DEG_IN_RADIAN; */ geocent(lst,geolat,elevsea,&x_geo,&y_geo,&z_geo); x = x - x_geo; /* topocentric correction using elliptical earth fig. */ y = y - y_geo; z = z - z_geo; *topodist = sqrt(x*x + y*y + z*z); l = x / (*topodist); m = y / (*topodist); n = z / (*topodist); *topora = atan_circ(l,m) * HRS_IN_RADIAN; *topodec = asin(n) * DEG_IN_RADIAN; } // This began as print_phase in skycalc.c. // No important changes in Skycal V5. static void find_next_moon_phase (double &jd, int &phase) { double newjd, lastnewjd, nextjd; short kount=0; // Originally, there was no problem with getting snagged, but since // I introduced the roundoff error going back and forth with Timestamp, // now it's a problem. // Move ahead by 1 second to avoid snagging. jd += 1.0 / SEC_IN_DAY; // Find current lunation. I doubled the safety margin since it // seemed biased for forwards search. Backwards search has since // been deleted, but little reason to mess with it. int nlast = (int)((jd - 2415020.5) / 29.5307 - 2); flmoon(nlast,0,&lastnewjd); flmoon(++nlast,0,&newjd); while (newjd <= jd) { lastnewjd = newjd; flmoon(++nlast,0,&newjd); require (kount++ < 5); // Original limit was 35 (!) } // We might save some work here by estimating, i.e.: // x = jd - lastnewjd; // noctiles = (int)(x / 3.69134); /* 3.69134 = 1/8 month; truncate. */ // However.... assert (lastnewjd <= jd && newjd > jd); phase = 1; // Lunation is lastnewjd's lunation flmoon (--nlast, phase, &nextjd); // Phase = 1 if (nextjd <= jd) { flmoon (nlast, ++phase, &nextjd); // Phase = 2 if (nextjd <= jd) { flmoon (nlast, ++phase, &nextjd); // Phase = 3 if (nextjd <= jd) { phase = 0; nextjd = newjd; } } } jd = nextjd; } void libxtide::Skycal::findNextMoonPhase (Timestamp t, TideEvent &tideEvent_out) { int phase; double jd (t.jd()); find_next_moon_phase (jd, phase); tideEvent_out.eventTime = jd; switch (phase) { case 0: tideEvent_out.eventType = TideEvent::newmoon; break; case 1: tideEvent_out.eventType = TideEvent::firstquarter; break; case 2: tideEvent_out.eventType = TideEvent::fullmoon; break; case 3: tideEvent_out.eventType = TideEvent::lastquarter; break; default: assert (false); } } // Here's another opportunity for Jeff Dairiki to write a better root // finder :-) // // jd_sun_alt needed good initial guesses to find sunrises and // sunsets. This was not a problem since good guesses were easy to // come by. The original skycalendar did this with estimates based on // the local midnight: // // jd = date_to_jd(date); /* local midnight */ // jdmid = jd + zone(use_dst,stdz,jd,jdbdst,jdedst) / 24.; /* corresponding ut */ // stmid = lst(jdmid,longit); // lpsun(jdmid,&rasun,&decsun); // hasunset = ha_alt(decsun,lat,-(0.83+horiz)); // jdsunset = jdmid + adj_time(rasun+hasunset-stmid)/24.; /* initial guess */ // jdsunset = jd_sun_alt(-(0.83+horiz), jdsunset,lat,longit); /* refinement */ // jdsunrise = jdmid + adj_time(rasun-hasunset-stmid)/24.; // jdsunrise = jd_sun_alt(-(0.83+horiz),jdsunrise,lat,longit); // // While efficient, this is an inconvenient way to go about it when // I'm looking for the next event from time t, and don't even know // when midnight is. So I messed with jd_sun_alt to make it converge // better from bad initial guesses, and substituted three bad guesses // for one good one. Normally, two would suffice, but I wanted to // add a safety margin in case one of them happens to land at a point // that nukes jd_sun_alt. // 2003-02-04 // Expanded to handle moonrise/moonset as well. // Set lunar to true for moonrise/set. static void find_next_rise_or_set (double &jd, double lat, double longit, bool lunar, bool &is_rise) { // Move ahead by precision interval to avoid snagging. jd += libxtide::Global::eventPrecisionJD; double jdorig = jd; double inc = 1.0 / 6.0; // 4 hours // First we want to know what we are looking for. bool looking_for = (altitude (jdorig, lat, longit, lunar) < riseAltitude); // Now give it a decent try. Because jd_alt is so unpredictable, // we can even find things out of order (which is one reason we need // to know what we're looking for). double jdlooper = jdorig; do { jd = jd_alt (riseAltitude, jdlooper, lat, longit, lunar, is_rise); jdlooper += inc; // Loop either on error return (which is a negative number), or if we // found an event in the wrong direction, or the wrong kind of event. } while ((jd < 0.0) || (jd <= jdorig) || (is_rise != looking_for)); } void libxtide::Skycal::findNextRiseOrSet (Timestamp t, const Coordinates &c, RiseSetType riseSetType, TideEvent &tideEvent_out) { assert (!(c.isNull())); bool isRise; double jd = t.jd(); // skycal "longit" is measured in HOURS WEST, not degrees east. // (lat is unchanged) find_next_rise_or_set (jd, c.lat(), -(c.lng())/15.0, (riseSetType==lunar), isRise); tideEvent_out.eventTime = jd; if (isRise) tideEvent_out.eventType = ((riseSetType == lunar) ? TideEvent::moonrise : TideEvent::sunrise); else tideEvent_out.eventType = ((riseSetType == lunar) ? TideEvent::moonset : TideEvent::sunset); } // Simple question deserving a simple answer... const bool libxtide::Skycal::sunIsUp (Timestamp t, const Coordinates &c) { assert (!(c.isNull())); return (altitude (t.jd(), c.lat(), -(c.lng())/15.0, 0) >= riseAltitude); } #ifdef EXPERIMENTAL_MOON_AGE_NOT_PHASE void libxtide::Skycal::findNewMoons (Timestamp t, Timestamp &prev_out, Timestamp &next_out) { double jd(t.jd()), newjd, lastnewjd; short kount=0; // Duplicated from find_next_moon_phase "find current lunation" block. // IDK whether 29.5307 really wants to be 29.530588853? int nlast = (int)((jd - 2415020.5) / 29.5307 - 2); flmoon(nlast,0,&lastnewjd); flmoon(++nlast,0,&newjd); while (newjd <= jd) { lastnewjd = newjd; flmoon(++nlast,0,&newjd); require (kount++ < 5); // Original limit was 35 (!) } // And presto, we're done. prev_out = lastnewjd; next_out = newjd; } #endif