Here is a new patch (to replace the original moonrise.patch) which uses the Stellarium implementation of ELP 2000-82b. As an added bonus, the final size of the executable has dropped to 1.1MB.
The file elp82b.c can be found here: https://github.com/Xitian9/Remind/blob/stellarium/src/elp82b.c I've editted it slightly from the version in Stellarium to remove some unneeded code and unit conversions. Best, Stephen On 8 February 2017 at 15:12, Stephen Morgan <[email protected]> wrote: > Hello again, > > Actually, perhaps a better solution is as follows. Stellarium is a big, > active, GPL planetarium software project, and their code includes an > implementation of ELP 2000-82b. It wouldn't necessarily relieve you of any > due diligence you feel necessary, but I would at least have more confidence > in the code being well-maintained and free from bugs. I could rewrite the > function to use their implementation instead; how does that sound? > > Their code can be downloaded from > https://sourceforge.net/p/stellarium/news/2016/12/stellarium-0151/ > and is located in the archive at stellarium-0.15-1/src/core/ > planetsephems/elp82b.{c,h}. > > By the way, as a matter of etiquette, I didn't feel I should be spamming > the list with 3.0MB files, but telling people to download a 150MB archive > just to extract two files inside is also problematic. Let me know if I > should send the file myself when I've finished. > > Best regards, > Stephen Morgan > > On 7 February 2017 at 10:49, Stephen Morgan <[email protected]> > wrote: > >> Hello everyone, >> >> I've written a couple of patches to add some moon-related functionality >> to remind. The first is a simple patch to create the moonillum function, >> which gets the percentage of the moon's surface which is illuminated at a >> given time. >> >> The second patch, which is a bit more troublesome, allows remind to get >> the times of moonrise and moonset on a given date, accurate to within 3 >> minutes of the United States Naval Observatory moonrise/set tables (RMS >> variation 0.7 minutes, at least for the dates I've checked). >> >> The only problem is that the implementation I've used adds a significant >> amount of code, mostly copied verbatim from this implementation of ELP >> 2000-82B (https://github.com/variar/elp2000-82b). This adds 2.7MB of >> code, and increases the size of the resulting binary from ~500kB to ~3.1MB. >> This problem could be solved by using a less accurate algorithm which is >> sufficient for these purposes (using ELP 2000-82B is rather like >> sandblasting a soup cracker), but that would require a bit more effort to >> find and code than I'm perhaps willing to put in. >> >> Alternately, a preprocessor flag could be added exclude the bulky code if >> you don't care about the moonrise/set times. >> >> To save on space in this email I've omitted the ELP 2000-82B code from >> the patch. Just copy it directly into a subdirectory elp2000-82b in the src >> directory. >> >> I hope others find these useful! >> >> Best regards, >> Stephen Morgan >> > >
diff --git a/man/remind.1 b/man/remind.1 index e34f8c0..73b6a07 100644 --- a/man/remind.1 +++ b/man/remind.1 @@ -2531,6 +2531,27 @@ which default to \fBtoday()\fR and midnight, respectively. The returned value is an integer from 0 to 359, representing the phase of the moon in degrees. 0 is a new moon, 180 is a full moon, 90 is first-quarter, etc. .TP +.B moonrise([dq_date]) +Returns a \fBTIME\fR indicating the time of moonrise on the specified +\fIdate\fR (default \fBtoday()\fR.) There may be no moonrise on a particular +day, once a month around first quarter, or because of of high latitudes, in +which case \fBmoonrise()\fR returns the \fBINT\fR 1440. +.TP +.B moonset([dq_date]) +Returns a \fBTIME\fR indicating the time of moonset on the specified +\fIdate\fR (default \fBtoday()\fR.) There may be no moonset on a particular +day, once a month around third quarter, or because of high latitudes, in which +case \fBmoonset()\fR returns the \fBINT\fR 0. +.RS +.PP +The functions \fBmoonrise()\fR and \fBmoonset()\fR are based on +an algorithm in "Explanatory supplement to the astronomical almanac" by the +Nautical Almanac Office, USNO. They require the latitude and longitude to be +specified by setting the appropriate system variables. (See "System +Variables".) The moon functions should be accurate to within about 3 minutes +for latitudes lower than 60 degrees. +.RE +.TP .B ndawn([dq_date]) Returns the time of "nautical dawn" on the specified \fIdate\fR. If \fIdate\fR is omitted, defaults to \fBtoday()\fR. If a \fIdatetime\fR object diff --git a/src/Makefile.in b/src/Makefile.in index 7207a4e..1495ade 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -27,11 +27,11 @@ MANS= $(srcdir)/../man/rem2ps.1 $(srcdir)/../man/remind.1 \ .SUFFIXES: .SUFFIXES: .c .o -REMINDSRCS= calendar.c dynbuf.c dorem.c dosubst.c expr.c files.c funcs.c \ - globals.c hbcal.c init.c main.c md5.c moon.c omit.c queue.c \ +REMINDSRCS= calendar.c dynbuf.c dorem.c dosubst.c elp82b.c expr.c files.c \ + funcs.c globals.c hbcal.c init.c main.c md5.c moon.c omit.c queue.c \ sort.c token.c trigger.c userfns.c utils.c var.c -REMINDHDRS=config.h custom.h dynbuf.h err.h expr.h globals.h lang.h \ +REMINDHDRS=config.h custom.h dynbuf.h elp82b.h err.h expr.h globals.h lang.h \ md5.h protos.h rem2ps.h types.h version.h REMINDOBJS= $(REMINDSRCS:.c=.o) diff --git a/src/elp82b.h b/src/elp82b.h new file mode 100644 index 0000000..58b9492 --- /dev/null +++ b/src/elp82b.h @@ -0,0 +1,79 @@ +/************************************************************************ + +LUNAR SOLUTION ELP2000-82B +by Chapront-Touze M., Chapront J. +ftp://ftp.imcce.fr/pub/ephem/moon/elp82b + +I (Johannes Gajdosik) have just taken the Fortran code and data +obtained from above and used it to create this piece of software. + +I can neither allow nor forbid the usage of ELP2000-82B. +The copyright notice below covers not the works of +Chapront-Touze M. and Chapront J., but just my work, +that is the compilation and rearrangement of +the Fortran code and data obtained from above +into the software supplied in this file. + + +Copyright (c) 2005 Johannes Gajdosik + +Permission is hereby granted, free of charge, to any person obtaining a +copy of this software and associated documentation files (the "Software"), +to deal in the Software without restriction, including without limitation +the rights to use, copy, modify, merge, publish, distribute, sublicense, +and/or sell copies of the Software, and to permit persons to whom the +Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +My implementation of ELP2000-82B has the following modifications compared to +the original Fortran code: +1) fundamentally rearrange the series into optimized instructions + for fast calculation of the results + +****************************************************************/ + + +#ifndef _ELP82B_H_ +#define _ELP82B_H_ + +void GetElp82bSphericalCoor(double t, double r[3]); + + /* Return the spherical coordinates of the earths moon + at time t expressed in julian centuries since J2000. The position is + expressed in arcseconds of longitude, arcseconds of latitude, and + kilometres of distance. The origin of the coordinates is the center of the + earth. The reference frame is "dynamical equinox and ecliptic J2000", + which is the reference frame in VSOP87 and VSOP87A. + + According to vsop87.doc VSOP87 coordinates can be transformed to FK5 by + X cos(psi) -sin(psi) 0 1 0 0 X + Y = sin(psi) cos(psi) 0 * 0 cos(eps) -sin(eps) * Y + Z FK5 0 0 1 0 sin(eps) cos(eps) Z VSOP87 + with psi = -0.0000275 degrees = -0.099 arcsec and + eps = 23.4392803055556 degrees = 23d26m21.4091sec. + + http://ssd.jpl.nasa.gov/horizons_doc.html#frames says: + "J2000" selects an Earth Mean-Equator and dynamical Equinox of + Epoch J2000.0 inertial reference system, where the Epoch of J2000.0 + is the Julian date 2451545.0. "Mean" indicates nutation effects are + ignored in the frame definition. The system is aligned with the + IAU-sponsored J2000 frame of the Radio Source Catalog of the + International Earth Rotational Service (ICRF). + The ICRF is thought to differ from FK5 by at most 0.01 arcsec. + + From this I conclude that in the context of stellarium + ICRF, J2000 and FK5 are the same, while the transformation + ICRF <-> VSOP87 must be done with the matrix given above. + */ + +#endif diff --git a/src/funcs.c b/src/funcs.c index 3d51171..4cf306f 100644 --- a/src/funcs.c +++ b/src/funcs.c @@ -99,6 +99,8 @@ static int FMoondate (func_info *); static int FMoondatetime (func_info *); static int FMoonillum (func_info *); static int FMoonphase (func_info *); +static int FMoonrise (func_info *); +static int FMoonset (func_info *); static int FMoontime (func_info *); static int FMon (func_info *); static int FMonnum (func_info *); @@ -240,6 +242,8 @@ BuiltinFunc Func[] = { { "moondatetime", 1, 3, 0, FMoondatetime }, { "moonillum", 0, 2, 0, FMoonillum }, { "moonphase", 0, 2, 0, FMoonphase }, + { "moonrise", 0, 2, 0, FMoonrise }, + { "moonset", 0, 2, 0, FMoonset }, { "moontime", 1, 3, 0, FMoontime }, { "ndawn", 0, 1, 0, FNDawn}, { "ndusk", 0, 1, 0, FNDusk}, @@ -2332,6 +2336,46 @@ static int FMoonillum(func_info *info) /***************************************************************/ /* */ +/* FMoonriseStuff */ +/* */ +/* Time of moonrise/set on a specified date/time. */ +/* */ +/***************************************************************/ +int FMoonriseStuff(int rise, func_info *info) +{ + int jul = JulianToday; + int r; + + if (Nargs >= 1) { + if (!HASDATE(ARG(0))) return E_BAD_TYPE; + jul = DATEPART(ARG(0)); + } + + r = MoonRise(rise, jul); + if (r == NO_TIME) { + RETVAL = 0; + RetVal.type = INT_TYPE; + } else if (r == -NO_TIME) { + RETVAL = MINUTES_PER_DAY; + RetVal.type = INT_TYPE; + } else { + RETVAL = r; + RetVal.type = TIME_TYPE; + } + return OK; +} + +static int FMoonrise(func_info *info) +{ + return FMoonriseStuff(1, info); +} +static int FMoonset(func_info *info) +{ + return FMoonriseStuff(0, info); +} + +/***************************************************************/ +/* */ /* FMoondate */ /* */ /* Hunt for next occurrence of specified moon phase */ diff --git a/src/moon.c b/src/moon.c index 6510ad6..04e0bb5 100644 --- a/src/moon.c +++ b/src/moon.c @@ -69,6 +69,7 @@ #include "expr.h" #include "globals.h" #include "err.h" +#include "elp82b.h" /* Function prototypes */ static long jdate (int y, int mon, int day); @@ -83,7 +84,7 @@ static double phase (double, double *, double *, double *, double *, double *, d /* Astronomical constants */ -#define epoch 2444238.5 /* 1980 January 0.0 */ +#define epoch 2444238.5 /* 1980 January 0.0 */ /* Constants defining the Sun's apparent orbit */ @@ -135,10 +136,12 @@ static double phase (double, double *, double *, double *, double *, double *, d #define abs(x) ((x) < 0 ? (-(x)) : (x)) /* Absolute val */ #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0))) /* Fix angle */ +#define fixhour(t) ((t) - 24.0 * (floor((t) / 24.0))) /* Fix hour */ #define torad(d) ((d) * (PI / 180.0)) /* Deg->Rad */ #define todeg(d) ((d) * (180.0 / PI)) /* Rad->Deg */ #define dsin(x) (sin(torad((x)))) /* Sin from deg */ #define dcos(x) (cos(torad((x)))) /* Cos from deg */ +#define dtan(x) (tan(torad((x)))) /* Tan from deg */ /***************************************************************/ /* */ @@ -524,7 +527,7 @@ int MoonPhase(int date, int time) FromJulian(utcd, &y, &m, &d); /* Convert to a true Julian date -- sorry for the name clashes! */ - jd = jtime(y, m, d, (utct / 60), (utct % 60), 0); + jd = jtime(y, m, d, (utct / 60), (utct % 60), 0); /* Calculate moon phase */ mp = 360.0 * phase(jd, NULL, NULL, NULL, NULL, NULL, NULL); @@ -556,7 +559,7 @@ void HuntPhase(int startdate, int starttim, int phas, int *date, int *time) /* Convert from Remind representation to year/mon/day */ FromJulian(utcd, &y, &m, &d); /* Convert to a true Julian date -- sorry for the name clashes! */ - jdorig = jtime(y, m, d, (utct / 60), (utct % 60), 0); + jdorig = jtime(y, m, d, (utct / 60), (utct % 60), 0); jd = jdorig - 45.0; nt1 = meanphase(jd, 0.0, &k1); while(1) { @@ -606,3 +609,133 @@ int MoonIllum(int date, int time) phase(jd, &il, NULL, NULL, NULL, NULL, NULL); return (int) (100.0 * il + 0.5); } + +/***************************************************************/ +/* */ +/* MoonRise -- Calculate time of moonrise/set: */ +/* */ +/* Algorithm from Explanatory supplement to the astronomical */ +/* almanac, from The nautical almanac office of the US naval */ +/* observatory, 1992, pp. 486-489. */ +/* */ +/***************************************************************/ +int MoonRise(int rise, int jul) +{ + int mins; + int year, mon, day; + + double latdeg, longdeg; + + double ec, ec0, d0; + double gst, gha; + double dec, ra, cosra; + double mp[3], l, b, e; + double par, a, ut, ut0, cost; + + int ncosfail, ndayfail; + + /* Get offset from UTC */ + if (CalculateUTC) { + if (CalcMinsFromUTC(jul, 12*60, &mins, NULL)) { + Eprint(ErrMsg[E_MKTIME_PROBLEM]); + return NO_TIME; + } + } else mins = MinsFromUTC; + + /* Get latitude and longitude */ + longdeg = (double) LongDeg + (double) LongMin / 60.0 + + (double) LongSec / 3600.0; + + latdeg = (double) LatDeg + (double) LatMin / 60.0 + + (double) LatSec / 3600.0; + + FromJulian(jul, &year, &mon, &day); + +/* Julian centuries since J2000; assumes a base of 1990 */ +#if BASE != 1990 +#error Moon calculations assume a BASE of 1990! +#endif + ec = (jul - 3652.5) / 36525; + e = 23.4392917 - 0.013005 * ec; /* Mean obliquity of the ecliptic. */ + + ut = 12.0; /* First iteration */ + ncosfail = 0; + ndayfail = 0; + + do { + ut0 = ut; + + /* Julian centuries since J2000 */ + ec0 = ec + (ut0 - (double) mins / 60.0) / 876600.0; + d0 = ec0 * 36525; /* Julian days since J2000 */ + + /* Calculate the moon's position */ + GetElp82bSphericalCoor(ec0, mp); + l = mp[0] / 3600.0; + b = mp[1] / 3600.0; + + /* Convert from ecliptic to equatorial co-ordinates */ + dec = todeg(asin( dsin(b) * dcos(e) + dcos(b) * dsin(e) * dsin(l) )); + cosra = dcos(l) * dcos(b) / dcos(dec); + + /* Check the sign of sin(ra), to determine which branch of acos to use */ + if (dcos(b) * dsin(l) * dcos(e) >= dsin(b) * dsin(e)) + ra = todeg(acos(cosra)); + else + ra = -todeg(acos(cosra)); + + /* Greenwich mean sidereal time in hours */ + gst = 18.697374558 + 24.06570982441908 * d0; + /* Greenwich hour angle of the moon in hours */ + gha = gst - ra / 15.0; + + /* Horizontal parallax of the moon.*/ + par = todeg(atan( 6378.137 / mp[2] )); + + /* Apparent altitude of upper limb of moon at moonrise/set */ + a = -34.0 / 60.0 + 0.7275 * par; + + /* Assuming the moon remains stationary, get the local hour angle of rise/set */ + cost = (dsin(a) - dsin(latdeg) * dsin(dec)) / (dcos(latdeg) * dcos(dec)); + + /* Check if the moon never rises/sets */ + if (cost > 1) { + cost = 1; /* Never sets */ + if (ncosfail < 0) ncosfail = 0; + ncosfail++; + } else if (cost < -1) { + cost = -1; /* Never rises */ + if (ncosfail > 0) ncosfail = 0; + ncosfail--; + } else + ncosfail = 0; + + /* Calculate UTC time of moonrise/set. */ + if (rise == 1) + ut = ut0 - gha - ( -longdeg + todeg(acos(cost)) ) / 15; + else + ut = ut0 - gha - ( -longdeg - todeg(acos(cost)) ) / 15; + ut = fixhour(ut); + + /* Check if we've rolled over to the next/previous day */ + if (abs(ut - ut0) > 12) + ndayfail++; + + /* Debugging + printf("Lat %.2f, Lon %.2f, Days %.1f, ut0 %09f, GHA %09f, t %f, ut %09f\n", + latdeg, longdeg, d0, ut0, fixhour(gha), todeg(acos(cost)), ut); + */ + + /* Fail with no rise/set */ + if (abs(ncosfail) > 3 || ndayfail > 1) { + if (rise) return -NO_TIME; + else return NO_TIME; + } + + } while (abs(ut - ut0) >= 0.008); + + mins = (int) floor(ut * 60 + 0.5); + if (mins == 1440) mins = 1439; + + return mins; +} diff --git a/src/protos.h b/src/protos.h index a1c7bc4..9bf6a1a 100644 --- a/src/protos.h +++ b/src/protos.h @@ -133,6 +133,8 @@ void LocalToUTC (int locdate, int loctime, int *utcdate, int *utctime); void UTCToLocal (int utcdate, int utctime, int *locdate, int *loctime); int MoonPhase (int date, int time); int MoonIllum (int date, int time); +int MoonIllum (int date, int time); +int MoonRise (int rise, int jul); void HuntPhase (int startdate, int starttim, int phas, int *date, int *time); int CompareRems (int dat1, int tim1, int prio1, int dat2, int tim2, int prio2, int bydate, int bytime, int byprio, int untimed_first); void SigIntHandler (int d);
_______________________________________________ Remind-fans mailing list [email protected] http://lists.roaringpenguin.com/cgi-bin/mailman/listinfo/remind-fans Remind is at http://www.roaringpenguin.com/products/remind
