// Compute times of moonrise and moonset at a specified latitude and longitude. // // This software minimizes computational work by performing the full calculation // of the lunar position three times, at the beginning, middle, and end of the // period of interest. Three point interpolation is used to predict the // position for each hour, and the arithmetic mean is used to predict the // half-hour positions. // // The full computational burden is negligible on modern computers, but the // algorithm is effective and still useful for small embedded systems. // // This software was originally adapted to javascript by Stephen R. Schmitt // from a BASIC program from the 'Astronomical Computing' column of Sky & // Telescope, April 1989, page 78. // // Subsequently adapted from Stephen R. Schmitt's javascript to c++ for the // Arduino by Cyrus Rahman. // // Subsequently adapted from Cyrus Rahman's Arduino C++ to C for the Sensor // Watch by hueso, this work is subject to Stephen Schmitt's copyright: // // Copyright 2007 Stephen R. Schmitt // Subsequent work Copyright 2020 Cyrus Rahman // You may use or modify this source code in any way you find useful, provided // that you agree that the author(s) have no warranty, obligations or liability. // You must determine the suitability of this source code for your use. // // Redistributions of this source code must retain this copyright notice. #include "moonrise.h" #include #define K1 15 * (M_PI / 180) * 1.0027379 // Determine the nearest moon rise or set event previous, and the nearest // moon rise or set event subsequent, to the specified time in seconds since the // Unix epoch (January 1, 1970) and at the specified latitude and longitude in // degrees. // // We look for events from MR_WINDOW/2 hours in the past to MR_WINDOW/2 hours // in the future. MoonRise MoonRise_calculate(double latitude, double longitude, uint32_t t) { MoonRise self = {}; skyCoordinates moonPosition[3]; double offsetDays; self.queryTime = t; offsetDays = julianDate(t) - 2451545L; // Days since Jan 1, 2000, 1200UTC. // Begin testing (MR_WINDOW / 2) hours before requested time. // offsetDays -= (double)MR_WINDOW / (2 * 24) ; // Calculate coordinates at start, middle, and end of search period. for (int i = 0; i < 3; i++) { moonPosition[i] = moon(offsetDays + i * (double)MR_WINDOW / (2 * 24)); } // If the RA wraps around during this period, unwrap it to keep the // sequence smooth for interpolation. if (moonPosition[1].RA <= moonPosition[0].RA) moonPosition[1].RA += 2 * M_PI; if (moonPosition[2].RA <= moonPosition[1].RA) moonPosition[2].RA += 2 * M_PI; // Initialize interpolation array. skyCoordinates mpWindow[3]; mpWindow[0].RA = moonPosition[0].RA; mpWindow[0].declination = moonPosition[0].declination; mpWindow[0].distance = moonPosition[0].distance; for (int k = 0; k < MR_WINDOW; k++) { // Check each interval of search period float ph = (float)(k + 1) / MR_WINDOW; mpWindow[2].RA = interpolate(moonPosition[0].RA, moonPosition[1].RA, moonPosition[2].RA, ph); mpWindow[2].declination = interpolate(moonPosition[0].declination, moonPosition[1].declination, moonPosition[2].declination, ph); mpWindow[2].distance = moonPosition[2].distance; // Look for moonrise/set events during this interval. { double ha[3], VHz[3]; double lSideTime; // Get (local_sidereal_time - MR_WINDOW / 2) hours in radians. lSideTime = localSiderealTime(offsetDays, longitude) * 2 * M_PI / 360; // Calculate Hour Angle. ha[0] = lSideTime - mpWindow[0].RA + k * K1; ha[2] = lSideTime - mpWindow[2].RA + k * K1 + K1; // Hour Angle and declination at half hour. ha[1] = (ha[2] + ha[0]) / 2; mpWindow[1].declination = (mpWindow[2].declination + mpWindow[0].declination) / 2; double s = sin(M_PI / 180 * latitude); double c = cos(M_PI / 180 * latitude); // refraction + semidiameter at horizon + distance correction double z = cos(M_PI / 180 * (90.567 - 41.685 / mpWindow[0].distance)); VHz[0] = s * sin(mpWindow[0].declination) + c * cos(mpWindow[0].declination) * cos(ha[0]) - z; VHz[2] = s * sin(mpWindow[2].declination) + c * cos(mpWindow[2].declination) * cos(ha[2]) - z; if (signbit(VHz[0]) == signbit(VHz[2])) goto noevent; // No event this hour. VHz[1] = s * sin(mpWindow[1].declination) + c * cos(mpWindow[1].declination) * cos(ha[1]) - z; double a, b, d, e, time; a = 2 * VHz[2] - 4 * VHz[1] + 2 * VHz[0]; b = 4 * VHz[1] - 3 * VHz[0] - VHz[2]; d = b * b - 4 * a * VHz[0]; if (d < 0) goto noevent; // No event this hour. d = sqrt(d); e = (-b + d) / (2 * a); if ((e < 0) || (e > 1)) e = (-b - d) / (2 * a); time = k + e + 1 / 120; // Time since k=0 of event (in hours). // The time we started searching + the time from the start of the search // to the event is the time of the event. uint32_t eventTime; eventTime = self.queryTime + (time) * 60 * 60; double hz, nz, dz, az; hz = ha[0] + e * (ha[2] - ha[0]); // Azimuth of the moon at the event. nz = -cos(mpWindow[1].declination) * sin(hz); dz = c * sin(mpWindow[1].declination) - s * cos(mpWindow[1].declination) * cos(hz); az = atan2(nz, dz) / (M_PI / 180); if (az < 0) az += 360; // If there is no previously recorded event of this type, save this event. // // If this event is previous to queryTime, and is the nearest event to // queryTime of events of its type previous to queryType, save this event, // replacing the previously recorded event of its type. Events subsequent // to queryTime are treated similarly, although since events are tested in // chronological order no replacements will occur as successive events // will be further from queryTime. // // If this event is subsequent to queryTime and there is an event of its // type previous to queryTime, then there is an event of the other type // between the two events of this event's type. If the event of the other // type is previous to queryTime, then it is the nearest event to // queryTime that is previous to queryTime. In this case save the current // event, replacing the previously recorded event of its type. Otherwise // discard the current event. // if ((VHz[0] < 0) && (VHz[2] > 0)) { if (!self.hasRise || ((self.riseTime < self.queryTime) == (eventTime < self.queryTime) && (self.riseTime - self.queryTime) > (eventTime - self.queryTime)) || ((self.riseTime < self.queryTime) != (eventTime < self.queryTime) && (self.hasSet && (self.riseTime < self.queryTime) == (self.setTime < self.queryTime)))) { self.riseTime = eventTime; self.riseAz = az; self.hasRise = true; } } if ((VHz[0] > 0) && (VHz[2] < 0)) { if (!self.hasSet || ((self.setTime < self.queryTime) == (eventTime < self.queryTime) && (self.setTime - self.queryTime) > (eventTime - self.queryTime)) || ((self.setTime < self.queryTime) != (eventTime < self.queryTime) && (self.hasRise && (self.setTime < self.queryTime) == (self.riseTime < self.queryTime)))) { self.setTime = eventTime; self.setAz = az; self.hasSet = true; } } noevent: // There are obscure cases in the polar regions that require extra logic. if (!self.hasRise && !self.hasSet) self.isVisible = !signbit(VHz[2]); else if (self.hasRise && !self.hasSet) self.isVisible = (self.queryTime > self.riseTime); else if (!self.hasRise && self.hasSet) self.isVisible = (self.queryTime < self.setTime); else self.isVisible = ((self.riseTime < self.setTime && self.riseTime < self.queryTime && self.setTime > self.queryTime) || (self.riseTime > self.setTime && (self.riseTime < self.queryTime || self.setTime > self.queryTime))); } if (self.hasSet && self.hasRise) break; mpWindow[0] = mpWindow[2]; // Advance to next interval. } return self; } // Moon position using fundamental arguments // (Van Flandern & Pulkkinen, 1979) skyCoordinates moon(double dayOffset) { double l = 0.606434 + 0.03660110129 * dayOffset; double m = 0.374897 + 0.03629164709 * dayOffset; double f = 0.259091 + 0.03674819520 * dayOffset; double d = 0.827362 + 0.03386319198 * dayOffset; double n = 0.347343 - 0.00014709391 * dayOffset; double g = 0.993126 + 0.00273777850 * dayOffset; l = 2 * M_PI * (l - floor(l)); m = 2 * M_PI * (m - floor(m)); f = 2 * M_PI * (f - floor(f)); d = 2 * M_PI * (d - floor(d)); n = 2 * M_PI * (n - floor(n)); g = 2 * M_PI * (g - floor(g)); double v, u, w; v = 0.39558 * sin(f + n) + 0.08200 * sin(f) + 0.03257 * sin(m - f - n) + 0.01092 * sin(m + f + n) + 0.00666 * sin(m - f) - 0.00644 * sin(m + f - 2*d + n) - 0.00331 * sin(f - 2*d + n) - 0.00304 * sin(f - 2*d) - 0.00240 * sin(m - f - 2*d - n) + 0.00226 * sin(m + f) - 0.00108 * sin(m + f - 2*d) - 0.00079 * sin(f - n) + 0.00078 * sin(f + 2*d + n); u = 1 - 0.10828 * cos(m) - 0.01880 * cos(m - 2*d) - 0.01479 * cos(2*d) + 0.00181 * cos(2*m - 2*d) - 0.00147 * cos(2*m) - 0.00105 * cos(2*d - g) - 0.00075 * cos(m - 2*d + g); w = 0.10478 * sin(m) - 0.04105 * sin(2*f + 2*n) - 0.02130 * sin(m - 2*d) - 0.01779 * sin(2*f + n) + 0.01774 * sin(n) + 0.00987 * sin(2*d) - 0.00338 * sin(m - 2*f - 2*n) - 0.00309 * sin(g) - 0.00190 * sin(2*f) - 0.00144 * sin(m + n) - 0.00144 * sin(m - 2*f - n) - 0.00113 * sin(m + 2*f + 2*n) - 0.00094 * sin(m - 2*d + g) - 0.00092 * sin(2*m - 2*d); double s; skyCoordinates sc; s = w / sqrt(u - v*v); sc.RA = l + atan(s / sqrt(1 - s*s)); // Right ascension s = v / sqrt(u); sc.declination = atan(s / sqrt(1 - s*s)); // Declination sc.distance = 60.40974 * sqrt(u); // Distance return(sc); } // 3-point interpolation double interpolate(double f0, double f1, double f2, double p) { double a = f1 - f0; double b = f2 - f1 - a; return(f0 + p * (2*a + b * (2*p - 1))); } // Determine Julian date from Unix time. // Provides marginally accurate results with Arduino 4-byte double. double julianDate(uint32_t t) { return (t / 86400.0L + 2440587.5); } // Local Sidereal Time // Provides local sidereal time in degrees, requires longitude in degrees // and time in fractional Julian days since Jan 1, 2000, 1200UTC (e.g. the // Julian date - 2451545). // cf. USNO Astronomical Almanac and // https://astronomy.stackexchange.com/questions/24859/local-sidereal-time double localSiderealTime(double offsetDays, double longitude) { double lSideTime = (15.0L * (6.697374558L + 0.06570982441908L * offsetDays + remainder(offsetDays, 1) * 24 + 12 + 0.000026 * (offsetDays / 36525) * (offsetDays / 36525)) + longitude) / 360; lSideTime -= floor(lSideTime); lSideTime *= 360; // Convert to degrees. return(lSideTime); }