diff options
Diffstat (limited to 'main/moon')
| -rw-r--r-- | main/moon | 601 |
1 files changed, 601 insertions, 0 deletions
diff --git a/main/moon b/main/moon new file mode 100644 index 0000000..0374f01 --- /dev/null +++ b/main/moon @@ -0,0 +1,601 @@ +diff --git a/movement/lib/sunriset/moonriset.c b/movement/lib/sunriset/moonriset.c +new file mode 100644 +index 0000000..ded9f33 +--- /dev/null ++++ b/movement/lib/sunriset/moonriset.c +@@ -0,0 +1,491 @@ ++#include <math.h> ++#include <stdio.h> ++#include "sunriset.h" ++ ++ ++#define E0BASE ((21.448 / 60. + 26.) / 60. + 23.) ++#define MOON_RADIUS 3474.8 ++ ++#define ARRAY_SIZE(a) (sizeof(a)/sizeof(a[0])) ++ ++#define days_since_2000_Jan_0(y,m,d) \ ++ (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L) ++ ++typedef struct ++{ ++ int d; ++ int m; ++ int mprime; ++ int f; ++ int csin; ++ int ccos; ++} Table; ++ ++static const Table table_47a[] = { ++ {0, 0, 1, 0, 6288774, -20905355}, ++ {2, 0, -1, 0, 1274027, -3699111}, ++ {2, 0, 0, 0, 658314, -2955968}, ++ {0, 0, 2, 0, 213618, -569925}, ++ {0, 1, 0, 0, -185116, 48888}, ++ {0, 0, 0, 2, -114332, -3149}, ++ {2, 0, -2, 0, 58793, 246158}, ++ {2, -1, -1, 0, 57066, -152138}, ++ {2, 0, 1, 0, 53322, -170733}, ++ {2, -1, 0, 0, 45758, -204586}, ++ {0, 1, -1, 0, -40923, -129620}, ++ {1, 0, 0, 0, -34720, 108743}, ++ {0, 1, 1, 0, -30383, 104755}, ++ {2, 0, 0, -2, 15327, 10321}, ++ {0, 0, 1, 2, -12528, 0}, ++ {0, 0, 1, -2, 10980, 79661}, ++ {4, 0, -1, 0, 10675, -34782}, ++ {0, 0, 3, 0, 10034, -23210}, ++ {4, 0, -2, 0, 8548, -21636}, ++ {2, 1, -1, 0, -7888, 24208}, ++ {2, 1, 0, 0, -6766, 30824}, ++ {1, 0, -1, 0, -5163, -8379}, ++ {1, 1, 0, 0, 4987, -16675}, ++ {2, -1, 1, 0, 4036, -12831}, ++ {2, 0, 2, 0, 3994, -10445}, ++ {4, 0, 0, 0, 3861, -11650}, ++ {2, 0, -3, 0, 3665, 14403}, ++ {0, 1, -2, 0, -2689, -7003}, ++ {2, 0, -1, 2, -2602, 0}, ++ {2, -1, -2, 0, 2390, 10056}, ++ {1, 0, 1, 0, -2348, 6322}, ++ {2, -2, 0, 0, 2236, -9884}, ++ {0, 1, 2, 0, -2120, 5751}, ++ {0, 2, 0, 0, -2069, 0}, ++ {2, -2, -1, 0, 2048, -4950}, ++ {2, 0, 1, -2, -1773, 4130}, ++ {2, 0, 0, 2, -1595, 0}, ++ {4, -1, -1, 0, 1215, -3958}, ++ {0, 0, 2, 2, -1110, 0}, ++ {3, 0, -1, 0, -892, 3258}, ++ {2, 1, 1, 0, -810, 2616}, ++ {4, -1, -2, 0, 759, -1897}, ++ {0, 2, -1, 0, -713, -2117}, ++ {2, 2, -1, 0, -700, 2354}, ++ {2, 1, -2, 0, 691, 0}, ++ {2, -1, 0, -2, 596, 0}, ++ {4, 0, 1, 0, 549, -1423}, ++ {0, 0, 4, 0, 537, -1117}, ++ {4, -1, 0, 0, 520, -1571}, ++ {1, 0, -2, 0, -487, -1739}, ++ {2, 1, 0, -2, -399, 0}, ++ {0, 0, 2, -2, -381, -4421}, ++ {1, 1, 1, 0, 351, 0}, ++ {3, 0, -2, 0, -340, 0}, ++ {4, 0, -3, 0, 330, 0}, ++ {2, -1, 2, 0, 327, 0}, ++ {0, 2, 1, 0, -323, 1165}, ++ {1, 1, -1, 0, 299, 0}, ++ {2, 0, 3, 0, 294, 0}, ++ {2, 0, -1, -2, 0, 8752} ++}; ++ ++static const Table table_47b[] = { ++ {0, 0, 0, 1, 5128122, 0}, ++ {0, 0, 1, 1, 280602, 0}, ++ {0, 0, 1, -1, 277693, 0}, ++ {2, 0, 0, -1, 173237, 0}, ++ {2, 0, -1, 1, 55413, 0}, ++ {2, 0, -1, -1, 46271, 0}, ++ {2, 0, 0, 1, 32573, 0}, ++ {0, 0, 2, 1, 17198, 0}, ++ {2, 0, 1, -1, 9266, 0}, ++ {0, 0, 2, -1, 8822, 0}, ++ {2, -1, 0, -1, 8216, 0}, ++ {2, 0, -2, -1, 4324, 0}, ++ {2, 0, 1, 1, 4200, 0}, ++ {2, 1, 0, -1, -3359, 0}, ++ {2, -1, -1, 1, 2463, 0}, ++ {2, -1, 0, 1, 2211, 0}, ++ {2, -1, -1, -1, 2065, 0}, ++ {0, 1, -1, -1, -1870, 0}, ++ {4, 0, -1, -1, 1828, 0}, ++ {0, 1, 0, 1, -1794, 0}, ++ {0, 0, 0, 3, -1749, 0}, ++ {0, 1, -1, 1, -1565, 0}, ++ {1, 0, 0, 1, -1491, 0}, ++ {0, 1, 1, 1, -1475, 0}, ++ {0, 1, 1, -1, -1410, 0}, ++ {0, 1, 0, -1, -1344, 0}, ++ {1, 0, 0, -1, -1335, 0}, ++ {0, 0, 3, 1, 1107, 0}, ++ {4, 0, 0, -1, 1021, 0}, ++ {4, 0, -1, 1, 833, 0}, ++ {0, 0, 1, -3, 777, 0}, ++ {4, 0, -2, 1, 671, 0}, ++ {2, 0, 0, -3, 607, 0}, ++ {2, 0, 2, -1, 596, 0}, ++ {2, -1, 1, -1, 491, 0}, ++ {2, 0, -2, 1, -451, 0}, ++ {0, 0, 3, -1, 439, 0}, ++ {2, 0, 2, 1, 422, 0}, ++ {2, 0, -3, -1, 421, 0}, ++ {2, 1, -1, 1, -366, 0}, ++ {2, 1, 0, 1, -351, 0}, ++ {4, 0, 0, 1, 331, 0}, ++ {2, -1, 1, 1, 315, 0}, ++ {2, -2, 0, -1, 302, 0}, ++ {0, 0, 1, 3, -283, 0}, ++ {2, 1, 1, -1, -229, 0}, ++ {1, 1, 0, -1, 223, 0}, ++ {1, 1, 0, 1, 223, 0}, ++ {0, 1, -2, -1, -220, 0}, ++ {2, 1, -1, -1, -220, 0}, ++ {1, 0, 1, 1, -185, 0}, ++ {2, -1, -2, -1, 181, 0}, ++ {0, 1, 2, 1, -177, 0}, ++ {4, 0, -2, -1, 176, 0}, ++ {4, -1, -1, -1, 166, 0}, ++ {1, 0, 1, -1, -164, 0}, ++ {4, 0, 1, -1, 132, 0}, ++ {1, 0, -1, -1, -119, 0}, ++ {4, -1, 0, -1, 115, 0}, ++ {2, -2, 0, 1, 107, 0}, ++}; ++ ++ ++static double ++deg2rad (double d) ++{ ++ return d * (M_PI / 180.); ++} ++ ++static double ++fold_deg2rad (double d) ++{ ++ d /= 360.; ++ d = d - floor (d); ++ return d * 2.0 * M_PI; ++} ++ ++ ++static void ++compute_table (const Table * t, unsigned n, double *ret_sin, double *ret_cos, ++ double D, double M, double Mprime, double F, double E, ++ double E2) ++{ ++ unsigned i; ++ double sigma_sin = 0; ++ double sigma_cos = 0; ++ ++ for (i = 0; i < n; ++i) ++ { ++ double csin = t[i].csin; ++ double ccos = t[i].ccos; ++ double arg = 0; ++ ++ switch (t[i].m) ++ { ++ case -2: ++ case 2: ++ csin *= E2; ++ ++ if (t[i].ccos) ++ ccos *= E2; ++ ++ break; ++ ++ case -1: ++ case 1: ++ csin *= E; ++ ++ if (t[i].ccos) ++ ccos *= E; ++ ++ break; ++ } ++ ++ if (t[i].d) ++ arg += D * (double) t[i].d; ++ ++ if (t[i].m) ++ arg += M * (double) t[i].m; ++ ++ if (t[i].mprime) ++ arg += Mprime * (double) t[i].mprime; ++ ++ if (t[i].f) ++ arg += F * (double) t[i].f; ++ ++ sigma_sin += csin * sin (arg); ++ ++ if (t[i].ccos) ++ sigma_cos += ccos * cos (arg); ++ } ++ ++ *ret_sin = sigma_sin; ++ ++ if (ret_cos) ++ *ret_cos = sigma_cos; ++} ++ ++ ++static int ++moon_above_horizon (double jd, double lst, double lat, int limb, double altit) ++{ ++ double T = (jd) / 36525.; ++ double L, Lprime, D, M, Mprime, F, A1, A2, A3, E, E2; ++ double sigma_l, sigma_r, sigma_b; ++ double lambda, beta, big_delta; ++ double omega, delta_psi, delta_epsilon; ++ double epsilon; ++ double alpha, delta; ++ ++ //printf("j0=%f T=%f lst=%f alit=%f\n",jd+2451545.,T,lst,altit); ++ ++ ++ // Sun's mean longitude ++ L = fold_deg2rad (36000.7698 * T + 280.4665); ++ ++ // Moon's mean longitude. ++ Lprime = fold_deg2rad ((((-(T / 65194000.) + 1. / 538841.) * T - 0.0015786) * T + 481267.88123421) * T + 218.3164477); // Meeus (47.1) ++ ++ // Moon's mean elongation. ++ D = fold_deg2rad ((((T / 113065000. + 1. / 545868.) * T - 0.0018819) * T + 445267.1114034) * T + 297.8501921); //Meeus (47.2) ++ ++ // Sun's mean anomaly. ++ M = fold_deg2rad (((T / 24490000. - 0.0001536) * T + 35999.0502909) * T + 357.5291092); // Meeus (47.3) ++ ++ // Moon's mean anomaly. ++ Mprime = fold_deg2rad ((((-T / 14712000. + 1. / 69699.) * T + 0.0087414) * T + 477198.8675055) * T + 134.9633964); // Meeus (47.4) ++ ++ // Moon's argument of latitude (mean distance ++ // from ascending node). ++ ++ F = fold_deg2rad ((((T / 863310000. - 1. / 3526000.) * T - 0.0036539) * T + 483202.0175233) * T + 93.2720950); // Meeus (47.5) ++ ++ // further arguments ++ A1 = fold_deg2rad (131.849 * T + 119.75); // Venus ++ A2 = fold_deg2rad (479264.290 * T + 53.09); // Jupiter ++ A3 = fold_deg2rad (481266.484 * T + 313.45); // undocumented ++ ++ // Eccentricity correction factor. ++ E = (-0.0000074 * T - 0.002516) * T + 1.; // Meeus (47.6) ++ E2 = E * E; ++ ++ ++ compute_table (table_47a, ARRAY_SIZE (table_47a), &sigma_l, &sigma_r, D, ++ M, Mprime, F, E, E2); ++ compute_table (table_47b, ARRAY_SIZE (table_47b), &sigma_b, (void *) 0, D, ++ M, Mprime, F, E, E2); ++ ++ sigma_l += 3958. * sin (A1) + 1962. * sin (Lprime - F) + 318. * sin (A2); ++ sigma_b += ++ -2235. * sin (Lprime) + 382. * sin (A3) + 175. * sin (A1 - F) + ++ 175. * sin (A1 + F) + 127. * sin (Lprime - Mprime) - ++ 115. * sin (Lprime + Mprime); ++ ++ // Coordinates of Moon (finally!) ++ ++ lambda = deg2rad (sigma_l / 1000000.) + Lprime; ++ beta = deg2rad (sigma_b / 1000000.); ++ big_delta = sigma_r / 1000. + 385000.56; ++ ++#if 0 ++ printf ("JDE=%+f\n", j0); ++ printf ("T =%+f\n", T); ++ printf ("L' =%+f\n", Lprime * (180 / M_PI)); ++ printf ("D =%+f\n", D * (180 / M_PI)); ++ printf ("M =%+f\n", M * (180 / M_PI)); ++ printf ("M' =%+f\n", Mprime * (180 / M_PI)); ++ printf ("F =%+f\n", F * (180 / M_PI)); ++ printf ("A1 =%+f\n", A1 * (180 / M_PI)); ++ printf ("A2 =%+f\n", A2 * (180 / M_PI)); ++ printf ("A3 =%+f\n", A3 * (180 / M_PI)); ++ printf ("E =%+f\n", E); ++ printf ("Sl =%+f\n", sigma_l); ++ printf ("Sb =%+f\n", sigma_b); ++ printf ("Sr =%+f\n", sigma_r); ++ ++ printf ("l=%+f\n", lambda * (180 / M_PI)); ++ printf ("b=%+f\n", beta * (180 / M_PI)); ++ printf ("D=%+f\n", big_delta); ++#endif ++ ++ ++ omega = ++ fold_deg2rad (((T / 450000 + .0020708) * T - 1934.136261) * T + ++ 125.04452); ++ ++ //budget - not doing the whole of chapter 22 ++ delta_psi = ++ deg2rad ((-17.20 * sin (omega) - 1.32 * sin (2 * L) - ++ 0.23 * sin (2 * Lprime) + 0.21 * sin (2 * omega)) / 3600.); ++ delta_epsilon = ++ deg2rad ((9.20 * cos (omega) + 0.57 * cos (2 * L) + ++ 0.10 * cos (2 * Lprime) - 0.09 * cos (2 * omega)) / 3600); ++ ++ lambda += delta_psi; ++ ++ epsilon = ++ deg2rad (((0.001813 * T - 0.00059) * T - 46.8150) * T / 3600 + ++ E0BASE); ++ epsilon += delta_epsilon; ++ ++#if 0 ++ printf ("delta_psi=%+f\n", delta_psi * (180 / M_PI)); ++ printf ("delta_epsilon=%+f\n", delta_epsilon * (180 / M_PI)); ++ printf ("epsilon=%+f\n", epsilon * (180 / M_PI)); ++#endif ++ ++ ++ alpha = ++ atan2 (sin (lambda) * cos (epsilon) - tan (beta) * sin (epsilon), ++ cos (lambda)); ++ delta = ++ asin (sin (beta) * cos (epsilon) + ++ cos (beta) * sin (epsilon) * sin (lambda)); ++ ++#if 0 ++ printf ("alpha=%+f\n", alpha * (180. / M_PI)); ++ printf ("delta=%+f\n", delta * (180. / M_PI)); ++#endif ++ ++ double H = deg2rad (lst) - alpha; ++ ++ lat = deg2rad (lat); ++ ++ double h = ++ asin (sin (lat) * sin (delta) + cos (lat) * cos (delta) * cos (H)); ++ ++// printf("rawh=%f\n",h); ++ ++ if (limb) ++ h += atan2 (MOON_RADIUS / 2, big_delta); ++ ++ h += deg2rad (altit); ++ ++ return h > 0; ++} ++ ++#if 0 ++double ++correct_refraction (double a) ++{ ++ ++ ++ if (a > ((-2.0 / 180.) * (M_PI))) ++ { ++ //double deg = a * (180./M_PI); ++ //double correction = deg2rad( (1.02 / tan (deg2rad (deg + 10.3/(deg + 5.11))) + .0019279)/60.); ++ //a+= deg2rad (correction / 60.); ++ a += 0.0002967 / tan (a + 0.00312536 / (a + 0.08901179)); ++ } ++ return a; ++} ++#endif ++ ++static double ++localSiderealTime (double jd, double longitude) ++{ ++ double lSideTime = (15.0L * (6.697374558L + 0.06570982441908L * jd + ++ remainder (jd, 1) * 24 + 12 + ++ 0.000026 * (jd / 36525) * (jd / 36525)) ++ - longitude) / 360; ++ lSideTime -= floor (lSideTime); ++ lSideTime *= 360; // Convert to degrees. ++//printf("lst %f %f => %f\n",jd+2451545.,longitude,lSideTime); ++ return (lSideTime); ++} ++ ++ ++int ++__moonriset__ (int year, int month, int day, int hour, int min, double lon, ++ double lat, double altit, int upper_limb, double *when) ++{ ++ double jd, jo, jt, o, step; ++ double lst; ++ int up, start_up = -1; ++ ++ jd = days_since_2000_Jan_0 (year, month, day) - 1.5; ++ ++ ++ jo = (double) min; ++ jo = (jo / 60.) + (double) hour; ++ jo = (jo / 24.); ++ ++ jt = jd + jo; ++ ++ ++ step = 0.1; ++ ++ for (o = 0; o < 2.5; o += step) ++ { ++ jt = jd + jo + o; ++ ++ ++ lst = localSiderealTime (jt, lon); ++ up = moon_above_horizon (jt, lst, lat, upper_limb, altit); ++ ++//printf("%f %f %d %d\n",jt,jt-jd,start_up,up); ++ ++ if (start_up == -1) ++ { ++ start_up = up; ++ continue; ++ } ++ ++ if (up != start_up) ++ { ++ if (step < 30. / 86400.) ++ break; ++ o -= step; ++ step *= 0.5; ++ } ++ } ++ ++ if (up == start_up) ++ { ++ return -1; ++ } ++//printf("moon goes to state %d at jd=%f jo=%f o=%f s=%f l=%f\n",up,jd,jo,o,jd+jo+o,jt); ++ *when = 24. * (o + jo); ++ return up; ++} ++ ++ ++#if 0 ++static int ++show_time (char *p, double s) ++{ ++ double h, m; ++ ++ if (s < 0) ++ s += 24.; ++ ++ h = floor (s); ++ s = (s - h) * 60.; ++ m = floor (s); ++ s = (s - m) * 60.; ++ ++ printf ("%s %02d:%02d:%08.5f\n", p, (int) h, (int) m, s); ++} ++ ++ ++int ++main (int argc, char *argv[]) ++{ ++ double rise, set; ++ ++ moonriset (2023, 9, 23, 15, 36, 0.1280, 52.2202, -19.152 / 60., 1, &rise, ++ &set); ++ show_time ("RISE ", rise); ++ show_time ("SET ", set); ++ ++ moonriset (2023, 9, 24, 0, 00, 0.1280, 52.2202, -19.152 / 60., 1, &rise, ++ &set); ++ show_time ("RISE ", rise); ++ show_time ("SET ", set); ++ ++ ++} ++ ++ ++ ++ ++ ++#endif +diff --git a/movement/lib/sunriset/sunriset.h b/movement/lib/sunriset/sunriset.h +index 7b2b6c8..8a79e88 100644 +--- a/movement/lib/sunriset/sunriset.h ++++ b/movement/lib/sunriset/sunriset.h +@@ -22,9 +22,13 @@ double __daylen__( int year, int month, int day, double lon, double lat, + double altit, int upper_limb ); + + int __sunriset__( int year, int month, int day, double lon, double lat, +- double altit, int upper_limb, double *rise, double *set ); ++ double altit, int upper_limb, double *rise, double *set); ++ ++int __moonriset__( int year, int month, int day, int hour, int min, double lon, double lat, ++ double altit, int upper_limb, double *when); + + void sun_RA_dec( double d, double *RA, double *dec, double *r ); ++void moon_RA_dec( double d, double *RA, double *dec, double *r ); + + double revolution( double x ); + +@@ -65,6 +69,12 @@ double GMST0( double d ); + __daylen__( year, month, day, lon, lat, -18.0, 0 ) + + ++ ++/* This macro computes times for moonrise/moonset. */ ++/* It's right according to meeus but it's not matching anytihng else */ ++#define moon_rise_set(year,month, day, hour,min, lon, lat, when) \ ++ __moonriset__( year, month, day, hour,min, lon, lat, -19.152/60.0, 1,when) ++ + /* This macro computes times for sunrise/sunset. */ + /* Sunrise/set is considered to occur when the Sun's upper limb is */ + /* 35 arc minutes below the horizon (this accounts for the refraction */ +diff --git a/movement/make/Makefile b/movement/make/Makefile +index 7110be2..db93b19 100644 +--- a/movement/make/Makefile ++++ b/movement/make/Makefile +@@ -41,6 +41,7 @@ SRCS += \ + ../lib/TOTP/TOTP.c \ + ../lib/base32/base32.c \ + ../lib/sunriset/sunriset.c \ ++ ../lib/sunriset/moonriset.c \ + ../lib/vsop87/vsop87a_milli.c \ + ../lib/astrolib/astrolib.c \ + ../lib/morsecalc/calc.c \ +diff --git a/movement/watch_faces/complication/sunrise_sunset_face.c b/movement/watch_faces/complication/sunrise_sunset_face.c +index 2a929c8..9bd3369 100644 +--- a/movement/watch_faces/complication/sunrise_sunset_face.c ++++ b/movement/watch_faces/complication/sunrise_sunset_face.c +@@ -78,6 +78,46 @@ static void _sunrise_sunset_face_update(movement_settings_t *settings, sunrise_s + // to deal with this, we set aside the offset in hours, and add it back before converting it to a watch_date_time. + double hours_from_utc = ((double)movement_timezone_offsets[settings->bit.time_zone]) / 60.0; + ++ ++ if (state->rise_index == 2) { ++ int wot; ++ watch_clear_display(); ++ watch_start_character_blink('C', 100); ++ ++ wot = moon_rise_set(scratch_time.unit.year + WATCH_RTC_REFERENCE_YEAR, scratch_time.unit.month, scratch_time.unit.day,scratch_time.unit.hour,scratch_time.unit.minute, lon, lat, &rise); ++ watch_stop_blink(); ++ ++ rise += hours_from_utc; ++ if (rise>24.) { ++ scratch_time.unit.day++; ++ rise -= 24.; ++ } ++ minutes = 60.0 * fmod(rise, 1); ++ seconds = 60.0 * fmod(minutes, 1); ++ ++ scratch_time.unit.hour = floor(rise); ++ ++ if (seconds < 30) scratch_time.unit.minute = floor(minutes); ++ else scratch_time.unit.minute = ceil(minutes); ++ ++ if (scratch_time.unit.minute == 60) { ++ scratch_time.unit.minute = 0; ++ scratch_time.unit.hour = (scratch_time.unit.hour + 1) % 24; ++ } ++ ++ watch_set_colon(); ++ ++ if (!settings->bit.clock_mode_24h) { ++ if (watch_utility_convert_to_12_hour(&scratch_time)) watch_set_indicator(WATCH_INDICATOR_PM); ++ else watch_clear_indicator(WATCH_INDICATOR_PM); ++ } ++ ++ sprintf(buf, "LU%2d%2d%02d%s", scratch_time.unit.day, scratch_time.unit.hour, scratch_time.unit.minute,wot ? "RI":"SE"); ++ watch_display_string(buf, 0); ++ return; ++ } ++ ++ + // we loop twice because if it's after sunset today, we need to recalculate to display values for tomorrow. + for(int i = 0; i < 2; i++) { + uint8_t result = sun_rise_set(scratch_time.unit.year + WATCH_RTC_REFERENCE_YEAR, scratch_time.unit.month, scratch_time.unit.day, lon, lat, &rise, &set); +@@ -369,7 +409,7 @@ bool sunrise_sunset_face_loop(movement_event_t event, movement_settings_t *setti + _sunrise_sunset_face_advance_digit(state); + _sunrise_sunset_face_update_settings_display(event, context); + } else { +- state->rise_index = (state->rise_index + 1) % 2; ++ state->rise_index = (state->rise_index + 1) % 3; + _sunrise_sunset_face_update(settings, state); + } + break; |