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 74425ff..e34f8c0 100644
--- a/man/remind.1
+++ b/man/remind.1
@@ -2520,6 +2520,11 @@ For example, the following returns the date and time of the next full moon:
This function is similar to \fBmoondate\fR and \fBmoontime\fR, but returns
a DATETIME result.
.TP
+.B moonillum([d_date [,t_time]]) \fRor\fB moonillum(q_datetime)
+This function returns the percentage of the moon which is illuminated on
+\fIdate\fR and \fItime\fR, which default to \fBtoday()\fR and midnight,
+respectively.
+.TP
.B moonphase([d_date [,t_time]]) \fRor\fB moonphase(q_datetime)
This function returns the phase of the moon on \fIdate\fR and \fItime\fR,
which default to \fBtoday()\fR and midnight, respectively. The returned
diff --git a/src/funcs.c b/src/funcs.c
index fbf6066..3d51171 100644
--- a/src/funcs.c
+++ b/src/funcs.c
@@ -97,6 +97,7 @@ static int FMinute (func_info *);
static int FMinsfromutc (func_info *);
static int FMoondate (func_info *);
static int FMoondatetime (func_info *);
+static int FMoonillum (func_info *);
static int FMoonphase (func_info *);
static int FMoontime (func_info *);
static int FMon (func_info *);
@@ -237,6 +238,7 @@ BuiltinFunc Func[] = {
{ "monnum", 1, 1, 1, FMonnum },
{ "moondate", 1, 3, 0, FMoondate },
{ "moondatetime", 1, 3, 0, FMoondatetime },
+ { "moonillum", 0, 2, 0, FMoonillum },
{ "moonphase", 0, 2, 0, FMoonphase },
{ "moontime", 1, 3, 0, FMoontime },
{ "ndawn", 0, 1, 0, FNDawn},
@@ -2290,6 +2292,46 @@ static int FMoonphase(func_info *info)
/***************************************************************/
/* */
+/* FMoonillum */
+/* */
+/* Illuminated percentage of moon for specified date/time. */
+/* */
+/***************************************************************/
+static int FMoonillum(func_info *info)
+{
+ int date, time;
+
+ switch(Nargs) {
+ case 0:
+ date = JulianToday;
+ time = 0;
+ break;
+ case 1:
+ if (!HASDATE(ARG(0))) return E_BAD_TYPE;
+ date = DATEPART(ARG(0));
+ if (HASTIME(ARG(0))) {
+ time = TIMEPART(ARG(0));
+ } else {
+ time = 0;
+ }
+ break;
+ case 2:
+ if (ARG(0).type == DATETIME_TYPE) return E_2MANY_ARGS;
+ if (ARG(0).type != DATE_TYPE && ARG(1).type != TIME_TYPE) return E_BAD_TYPE;
+ date = ARGV(0);
+ time = ARGV(1);
+ break;
+
+ default: return E_SWERR;
+ }
+
+ RetVal.type = INT_TYPE;
+ RETVAL = MoonIllum(date, time);
+ return OK;
+}
+
+/***************************************************************/
+/* */
/* FMoondate */
/* */
/* Hunt for next occurrence of specified moon phase */
diff --git a/src/moon.c b/src/moon.c
index 35d1d2a..6510ad6 100644
--- a/src/moon.c
+++ b/src/moon.c
@@ -577,3 +577,32 @@ void HuntPhase(int startdate, int starttim, int phas, int *date, int *time)
t1 = h*60 + min;
UTCToLocal(d1, t1, date, time);
}
+
+/***************************************************************/
+/* */
+/* MoonIllum */
+/* */
+/* Interface routine dealing in Remind representations. */
+/* Given a local date and time, returns the percentage */
+/* of the moon illuminated. */
+/* */
+/***************************************************************/
+int MoonIllum(int date, int time)
+{
+ int utcd, utct;
+ int y, m, d;
+ double jd, il;
+
+ /* Convert from local to UTC */
+ LocalToUTC(date, time, &utcd, &utct);
+
+ /* Convert from Remind representation to year/mon/day */
+ 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);
+
+ /* Calculate moon phase */
+ phase(jd, &il, NULL, NULL, NULL, NULL, NULL);
+ return (int) (100.0 * il + 0.5);
+}
diff --git a/src/protos.h b/src/protos.h
index 9179348..a1c7bc4 100644
--- a/src/protos.h
+++ b/src/protos.h
@@ -132,6 +132,7 @@ void FillParagraph (char const *s);
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);
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);
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..c35b4a9 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -29,10 +29,12 @@ MANS= $(srcdir)/../man/rem2ps.1 $(srcdir)/../man/remind.1 \
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 \
- sort.c token.c trigger.c userfns.c utils.c var.c
+ sort.c token.c trigger.c userfns.c utils.c var.c \
+ elp2000-82b/elp2000-82b.c elp2000-82b/arguments.c elp2000-82b/series.c
REMINDHDRS=config.h custom.h dynbuf.h err.h expr.h globals.h lang.h \
- md5.h protos.h rem2ps.h types.h version.h
+ md5.h protos.h rem2ps.h types.h version.h \
+ elp2000-82b/elp2000-82b.h elp2000-82b/arguments.h elp2000-82b/series.h
REMINDOBJS= $(REMINDSRCS:.c=.o)
all: remind rem2ps
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..2d1632f 100644
--- a/src/moon.c
+++ b/src/moon.c
@@ -70,6 +70,8 @@
#include "globals.h"
#include "err.h"
+#include "elp2000-82b/elp2000-82b.h"
+
/* Function prototypes */
static long jdate (int y, int mon, int day);
static double jtime (int y, int mon, int day, int hour, int min, int sec);
@@ -83,7 +85,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 +137,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 +528,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 +560,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 +610,134 @@ 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 l, b, e;
+ double par, a, ut, ut0, cost;
+ spherical_point mp;
+
+ 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 */
+ mp = geocentric_moon_position(ec0);
+ l = mp.longitude / 3600.0;
+ b = mp.latitude / 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.distance ));
+
+ /* 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
