|
|
|
|
@@ -2,36 +2,35 @@
|
|
|
|
|
//
|
|
|
|
|
// 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.
|
|
|
|
|
// 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.
|
|
|
|
|
// 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 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:
|
|
|
|
|
// 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
|
|
|
|
|
// 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.
|
|
|
|
|
// 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 <math.h>
|
|
|
|
|
#include <stdbool.h>
|
|
|
|
|
#include <stdlib.h>
|
|
|
|
|
#include "moonrise.h"
|
|
|
|
|
#include <math.h>
|
|
|
|
|
|
|
|
|
|
#define K1 15*(M_PI/180)*1.0027379
|
|
|
|
|
#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
|
|
|
|
|
@@ -40,15 +39,15 @@
|
|
|
|
|
//
|
|
|
|
|
// 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, time_t t) {
|
|
|
|
|
MoonRise self;
|
|
|
|
|
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.
|
|
|
|
|
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) ;
|
|
|
|
|
// offsetDays -= (double)MR_WINDOW / (2 * 24) ;
|
|
|
|
|
|
|
|
|
|
// Calculate coordinates at start, middle, and end of search period.
|
|
|
|
|
for (int i = 0; i < 3; i++) {
|
|
|
|
|
@@ -64,142 +63,150 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
|
|
|
|
|
|
|
|
|
|
// Initialize interpolation array.
|
|
|
|
|
skyCoordinates mpWindow[3];
|
|
|
|
|
mpWindow[0].RA = moonPosition[0].RA;
|
|
|
|
|
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);
|
|
|
|
|
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;
|
|
|
|
|
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;
|
|
|
|
|
// 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;
|
|
|
|
|
// 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;
|
|
|
|
|
// 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);
|
|
|
|
|
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));
|
|
|
|
|
// 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;
|
|
|
|
|
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.
|
|
|
|
|
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;
|
|
|
|
|
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];
|
|
|
|
|
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.
|
|
|
|
|
if (d < 0)
|
|
|
|
|
goto noevent; // No event this hour.
|
|
|
|
|
|
|
|
|
|
d = sqrt(d);
|
|
|
|
|
e = (-b + d) / (2 * a);
|
|
|
|
|
if ((e < 0) || (e > 1))
|
|
|
|
|
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).
|
|
|
|
|
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.
|
|
|
|
|
time_t eventTime;
|
|
|
|
|
eventTime = self.queryTime + (time) *60 *60;
|
|
|
|
|
// 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)
|
|
|
|
|
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 ||
|
|
|
|
|
// 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) &&
|
|
|
|
|
llabs(self.riseTime - self.queryTime) > llabs(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;
|
|
|
|
|
}
|
|
|
|
|
(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 ||
|
|
|
|
|
}
|
|
|
|
|
if ((VHz[0] > 0) && (VHz[2] < 0)) {
|
|
|
|
|
if (!self.hasSet ||
|
|
|
|
|
((self.setTime < self.queryTime) == (eventTime < self.queryTime) &&
|
|
|
|
|
llabs(self.setTime - self.queryTime) > llabs(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;
|
|
|
|
|
}
|
|
|
|
|
(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)
|
|
|
|
|
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)
|
|
|
|
|
else if (self.hasRise && !self.hasSet)
|
|
|
|
|
self.isVisible = (self.queryTime > self.riseTime);
|
|
|
|
|
else if (!self.hasRise && self.hasSet)
|
|
|
|
|
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)));
|
|
|
|
|
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.
|
|
|
|
|
mpWindow[0] = mpWindow[2]; // Advance to next interval.
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
return self;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// Moon position using fundamental arguments
|
|
|
|
|
// Moon position using fundamental arguments
|
|
|
|
|
// (Van Flandern & Pulkkinen, 1979)
|
|
|
|
|
skyCoordinates moon(double dayOffset) {
|
|
|
|
|
double l = 0.606434 + 0.03660110129 * dayOffset;
|
|
|
|
|
@@ -275,7 +282,7 @@ double interpolate(double f0, double f1, double f2, double p) {
|
|
|
|
|
|
|
|
|
|
// Determine Julian date from Unix time.
|
|
|
|
|
// Provides marginally accurate results with Arduino 4-byte double.
|
|
|
|
|
double julianDate(time_t t) {
|
|
|
|
|
double julianDate(uint32_t t) {
|
|
|
|
|
return (t / 86400.0L + 2440587.5);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|