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 +#include +#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;