This commit is contained in:
hueso 2025-05-16 18:41:17 -03:00
parent 48399312be
commit 13863d32ca
3 changed files with 151 additions and 151 deletions

View File

@ -2,34 +2,35 @@
// //
// This software minimizes computational work by performing the full calculation // This software minimizes computational work by performing the full calculation
// of the lunar position three times, at the beginning, middle, and end of the // 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 // period of interest. Three point interpolation is used to predict the
// for each hour, and the arithmetic mean is used to predict the half-hour positions. // 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 // The full computational burden is negligible on modern computers, but the
// algorithm is effective and still useful for small embedded systems. // algorithm is effective and still useful for small embedded systems.
// //
// This software was originally adapted to javascript by Stephen R. Schmitt // This software was originally adapted to javascript by Stephen R. Schmitt
// from a BASIC program from the 'Astronomical Computing' column of Sky & Telescope, // from a BASIC program from the 'Astronomical Computing' column of Sky &
// April 1989, page 78. // Telescope, April 1989, page 78.
// //
// Subsequently adapted from Stephen R. Schmitt's javascript to c++ for the Arduino // Subsequently adapted from Stephen R. Schmitt's javascript to c++ for the
// by Cyrus Rahman. // Arduino by Cyrus Rahman.
// //
// Subsequently adapted from Cyrus Rahman's Arduino C++ to C for the Sensor Watch // Subsequently adapted from Cyrus Rahman's Arduino C++ to C for the Sensor
// by hueso, this work is subject to Stephen Schmitt's copyright: // 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 // Subsequent work Copyright 2020 Cyrus Rahman
// You may use or modify this source code in any way you find useful, provided // 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 // that you agree that the author(s) have no warranty, obligations or liability.
// must determine the suitability of this source code for your use. // You must determine the suitability of this source code for your use.
// //
// Redistributions of this source code must retain this copyright notice. // Redistributions of this source code must retain this copyright notice.
#include "moonrise.h"
#include <math.h> #include <math.h>
#include <stdbool.h> #include <stdbool.h>
#include <stdlib.h> #include <stdlib.h>
#include "moonrise.h"
#define K1 15 * (M_PI / 180) * 1.0027379 #define K1 15 * (M_PI / 180) * 1.0027379
@ -40,8 +41,8 @@
// //
// We look for events from MR_WINDOW/2 hours in the past to MR_WINDOW/2 hours // We look for events from MR_WINDOW/2 hours in the past to MR_WINDOW/2 hours
// in the future. // in the future.
MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) { MoonRise MoonRise_calculate(double latitude, double longitude, uint32_t t) {
MoonRise self; MoonRise self = {};
skyCoordinates moonPosition[3]; skyCoordinates moonPosition[3];
double offsetDays; double offsetDays;
@ -71,11 +72,10 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
for (int k = 0; k < MR_WINDOW; k++) { // Check each interval of search period for (int k = 0; k < MR_WINDOW; k++) { // Check each interval of search period
float ph = (float)(k + 1) / MR_WINDOW; float ph = (float)(k + 1) / MR_WINDOW;
mpWindow[2].RA = interpolate(moonPosition[0].RA, mpWindow[2].RA = interpolate(moonPosition[0].RA, moonPosition[1].RA,
moonPosition[1].RA,
moonPosition[2].RA, ph); moonPosition[2].RA, ph);
mpWindow[2].declination = interpolate(moonPosition[0].declination, mpWindow[2].declination =
moonPosition[1].declination, interpolate(moonPosition[0].declination, moonPosition[1].declination,
moonPosition[2].declination, ph); moonPosition[2].declination, ph);
mpWindow[2].distance = moonPosition[2].distance; mpWindow[2].distance = moonPosition[2].distance;
@ -93,7 +93,8 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
// Hour Angle and declination at half hour. // Hour Angle and declination at half hour.
ha[1] = (ha[2] + ha[0]) / 2; ha[1] = (ha[2] + ha[0]) / 2;
mpWindow[1].declination = (mpWindow[2].declination + mpWindow[0].declination)/2; mpWindow[1].declination =
(mpWindow[2].declination + mpWindow[0].declination) / 2;
double s = sin(M_PI / 180 * latitude); double s = sin(M_PI / 180 * latitude);
double c = cos(M_PI / 180 * latitude); double c = cos(M_PI / 180 * latitude);
@ -101,13 +102,16 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
// refraction + semidiameter at horizon + distance correction // refraction + semidiameter at horizon + distance correction
double z = cos(M_PI / 180 * (90.567 - 41.685 / mpWindow[0].distance)); 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[0] = s * sin(mpWindow[0].declination) +
VHz[2] = s * sin(mpWindow[2].declination) + c * cos(mpWindow[2].declination) * cos(ha[2]) - z; 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])) if (signbit(VHz[0]) == signbit(VHz[2]))
goto noevent; // No event this hour. 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; double a, b, d, e, time;
a = 2 * VHz[2] - 4 * VHz[1] + 2 * VHz[0]; a = 2 * VHz[2] - 4 * VHz[1] + 2 * VHz[0];
@ -123,43 +127,44 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
e = (-b - d) / (2 * a); 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 // The time we started searching + the time from the start of the search
// event is the time of the event. // to the event is the time of the event.
time_t eventTime; uint32_t eventTime;
eventTime = self.queryTime + (time) * 60 * 60; eventTime = self.queryTime + (time) * 60 * 60;
double hz, nz, dz, az; double hz, nz, dz, az;
hz = ha[0] + e * (ha[2] - ha[0]); // Azimuth of the moon at the event. hz = ha[0] + e * (ha[2] - ha[0]); // Azimuth of the moon at the event.
nz = -cos(mpWindow[1].declination) * sin(hz); nz = -cos(mpWindow[1].declination) * sin(hz);
dz = c * sin(mpWindow[1].declination) - s * cos(mpWindow[1].declination) * cos(hz); dz = c * sin(mpWindow[1].declination) -
s * cos(mpWindow[1].declination) * cos(hz);
az = atan2(nz, dz) / (M_PI / 180); az = atan2(nz, dz) / (M_PI / 180);
if (az < 0) if (az < 0)
az += 360; az += 360;
// If there is no previously recorded event of this type, save this event. // 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 // If this event is previous to queryTime, and is the nearest event to
// of events of its type previous to queryType, save this event, replacing the // queryTime of events of its type previous to queryType, save this event,
// previously recorded event of its type. Events subsequent to queryTime are // replacing the previously recorded event of its type. Events subsequent
// treated similarly, although since events are tested in chronological order // to queryTime are treated similarly, although since events are tested in
// no replacements will occur as successive events will be further from // chronological order no replacements will occur as successive events
// queryTime. // will be further from queryTime.
// //
// If this event is subsequent to queryTime and there is an event of its type // If this event is subsequent to queryTime and there is an event of its
// previous to queryTime, then there is an event of the other type between the // type previous to queryTime, then there is an event of the other type
// two events of this event's type. If the event of the other type is // between the two events of this event's type. If the event of the other
// previous to queryTime, then it is the nearest event to queryTime that is // type is previous to queryTime, then it is the nearest event to
// previous to queryTime. In this case save the current event, replacing // queryTime that is previous to queryTime. In this case save the current
// the previously recorded event of its type. Otherwise discard the current // event, replacing the previously recorded event of its type. Otherwise
// event. // discard the current event.
// //
if ((VHz[0] < 0) && (VHz[2] > 0)) { if ((VHz[0] < 0) && (VHz[2] > 0)) {
if (!self.hasRise || if (!self.hasRise ||
((self.riseTime < self.queryTime) == (eventTime < self.queryTime) && ((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.riseTime < self.queryTime) != (eventTime < self.queryTime) &&
(self.hasSet && (self.hasSet && (self.riseTime < self.queryTime) ==
(self.riseTime < self.queryTime) == (self.setTime < self.queryTime)))) { (self.setTime < self.queryTime)))) {
self.riseTime = eventTime; self.riseTime = eventTime;
self.riseAz = az; self.riseAz = az;
self.hasRise = true; self.hasRise = true;
@ -168,10 +173,10 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
if ((VHz[0] > 0) && (VHz[2] < 0)) { if ((VHz[0] > 0) && (VHz[2] < 0)) {
if (!self.hasSet || if (!self.hasSet ||
((self.setTime < self.queryTime) == (eventTime < self.queryTime) && ((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.setTime < self.queryTime) != (eventTime < self.queryTime) &&
(self.hasRise && (self.hasRise && (self.setTime < self.queryTime) ==
(self.setTime < self.queryTime) == (self.riseTime < self.queryTime)))) { (self.riseTime < self.queryTime)))) {
self.setTime = eventTime; self.setTime = eventTime;
self.setAz = az; self.setAz = az;
self.hasSet = true; self.hasSet = true;
@ -187,18 +192,22 @@ MoonRise MoonRise_calculate(double latitude, double longitude, time_t t) {
else if (!self.hasRise && self.hasSet) else if (!self.hasRise && self.hasSet)
self.isVisible = (self.queryTime < self.setTime); self.isVisible = (self.queryTime < self.setTime);
else else
self.isVisible = ((self.riseTime < self.setTime && self.riseTime < self.queryTime && self.setTime > self.queryTime) || 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) ||
(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; return self;
} }
// Moon position using fundamental arguments // Moon position using fundamental arguments
// (Van Flandern & Pulkkinen, 1979) // (Van Flandern & Pulkkinen, 1979)
skyCoordinates moon(double dayOffset) { skyCoordinates moon(double dayOffset) {
@ -275,7 +284,7 @@ double interpolate(double f0, double f1, double f2, double p) {
// Determine Julian date from Unix time. // Determine Julian date from Unix time.
// Provides marginally accurate results with Arduino 4-byte double. // 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); return (t / 86400.0L + 2440587.5);
} }

View File

@ -1,8 +1,7 @@
#ifndef MoonRise_h #ifndef MoonRise_h
#define MoonRise_h #define MoonRise_h
#include <time.h> #include <stdint.h>
// Size of event search window in hours. // Size of event search window in hours.
// Events further away from the search time than MR_WINDOW/2 will not be // Events further away from the search time than MR_WINDOW/2 will not be
// found. At higher latitudes the moon rise/set intervals become larger, so if // found. At higher latitudes the moon rise/set intervals become larger, so if
@ -10,7 +9,7 @@
// windows will increase interpolation error. Useful values are probably from // windows will increase interpolation error. Useful values are probably from
// 12 - 48 but will depend upon your application. // 12 - 48 but will depend upon your application.
#define MR_WINDOW 72 // Even integer #define MR_WINDOW 48 // Even integer
typedef struct { typedef struct {
double RA; // Right ascension double RA; // Right ascension
@ -19,9 +18,9 @@ typedef struct {
} skyCoordinates; } skyCoordinates;
typedef struct { typedef struct {
time_t queryTime; uint32_t queryTime;
time_t riseTime; uint32_t riseTime;
time_t setTime; uint32_t setTime;
float riseAz; float riseAz;
float setAz; float setAz;
bool hasRise; bool hasRise;
@ -29,14 +28,14 @@ typedef struct {
bool isVisible; bool isVisible;
} MoonRise; } MoonRise;
MoonRise MoonRise_calculate(double latitude, double longitude, time_t t); MoonRise MoonRise_calculate(double latitude, double longitude, uint32_t t);
// private: // private:
void testMoonRiseSet(MoonRise *self, int i, double offsetDays, double latitude, double longitude, void testMoonRiseSet(MoonRise *self, int i, double offsetDays, double latitude,
skyCoordinates *mp); double longitude, skyCoordinates *mp);
skyCoordinates moon(double dayOffset); skyCoordinates moon(double dayOffset);
double interpolate(double f0, double f1, double f2, double p); double interpolate(double f0, double f1, double f2, double p);
double julianDate(time_t t); double julianDate(uint32_t t);
double localSiderealTime(double offsetDays, double longitude); double localSiderealTime(double offsetDays, double longitude);
#endif #endif

View File

@ -52,7 +52,7 @@ static void _moonrise_face_update(movement_settings_t *settings, moonrise_state_
if (movement_location.reg == 0) { if (movement_location.reg == 0) {
watch_clear_colon(); watch_clear_colon();
watch_display_string("MR no Loc", 0); watch_display_string("Mz no Loc", 0);
return; return;
} }
@ -60,15 +60,10 @@ static void _moonrise_face_update(movement_settings_t *settings, moonrise_state_
watch_date_time scratch_time; // scratchpad, contains different values at different times watch_date_time scratch_time; // scratchpad, contains different values at different times
scratch_time.reg = date_time.reg; scratch_time.reg = date_time.reg;
// Weird quirky unsigned things were happening when I tried to cast these directly to doubles below. double lat = (double)movement_location.bit.latitude / 100.0;
// it looks redundant, but extracting them to local int16's seemed to fix it. double lon = (double)movement_location.bit.longitude / 100.0;
int16_t lat_centi = (int16_t)movement_location.bit.latitude;
int16_t lon_centi = (int16_t)movement_location.bit.longitude;
double lat = (double)lat_centi / 100.0; uint32_t t = watch_utility_date_time_to_unix_time(date_time, movement_timezone_offsets[settings->bit.time_zone] * 60);
double lon = (double)lon_centi / 100.0;
time_t t = watch_utility_date_time_to_unix_time(scratch_time, movement_timezone_offsets[settings->bit.time_zone] * 60);
MoonRise mr = MoonRise_calculate(lat, lon, t); MoonRise mr = MoonRise_calculate(lat, lon, t);
if(mr.isVisible) if(mr.isVisible)
@ -76,11 +71,12 @@ static void _moonrise_face_update(movement_settings_t *settings, moonrise_state_
else else
watch_clear_indicator(WATCH_INDICATOR_LAP); watch_clear_indicator(WATCH_INDICATOR_LAP);
if (!mr.hasRise && !mr.hasSet) { if ( (state->rise_index == 0 && !mr.hasRise) ||
(state->rise_index == 1 && !mr.hasSet) ) {
watch_clear_colon(); watch_clear_colon();
watch_clear_indicator(WATCH_INDICATOR_PM); watch_clear_indicator(WATCH_INDICATOR_PM);
watch_clear_indicator(WATCH_INDICATOR_24H); watch_clear_indicator(WATCH_INDICATOR_24H);
snprintf(buf, sizeof(buf), "MR%2d none ", scratch_time.unit.day); snprintf(buf, sizeof(buf), "%s%2d none ", state->rise_index ? "M_" : "M~", scratch_time.unit.day);
watch_display_string(buf, 0); watch_display_string(buf, 0);
return; return;
} }
@ -95,7 +91,6 @@ static void _moonrise_face_update(movement_settings_t *settings, moonrise_state_
state->rise_set_expires.reg = scratch_time.reg; state->rise_set_expires.reg = scratch_time.reg;
bool set_leading_zero = false; bool set_leading_zero = false;
if (!settings->bit.clock_mode_24h) if (!settings->bit.clock_mode_24h)
if (watch_utility_convert_to_12_hour(&scratch_time)) if (watch_utility_convert_to_12_hour(&scratch_time))
@ -105,12 +100,9 @@ static void _moonrise_face_update(movement_settings_t *settings, moonrise_state_
else if (settings->bit.clock_24h_leading_zero && scratch_time.unit.hour < 10) { else if (settings->bit.clock_24h_leading_zero && scratch_time.unit.hour < 10) {
set_leading_zero = true; set_leading_zero = true;
} }
snprintf(buf, sizeof(buf), "M %2d%2d%02d%2s", scratch_time.unit.day, scratch_time.unit.hour, scratch_time.unit.minute,longLatPresets[state->longLatToUse].name); snprintf(buf, sizeof(buf), "%s%2d%2d%02d%2s", state->rise_index ? "M_" : "M~", scratch_time.unit.day, scratch_time.unit.hour, scratch_time.unit.minute,longLatPresets[state->longLatToUse].name);
watch_display_string(buf, 0); watch_display_string(buf, 0);
if(state->rise_index == 0)
watch_set_pixel(0,11);
else
watch_set_pixel(2,11);
if (set_leading_zero) if (set_leading_zero)
watch_display_string("0", 4); watch_display_string("0", 4);
return; return;