diff --git a/config/library.csv b/config/library.csv index d20497f4..8b8db5f8 100644 --- a/config/library.csv +++ b/config/library.csv @@ -9,6 +9,10 @@ "SiDensity", "« → T '(8.35123E20 (cm↑-3))*exp((-7 555.17 K)÷T)' »" "Fanning", "« → Roughness Reynolds 'IFTE(Reynolds≤2100;16/Reynolds;(4.781-((-2۰log((Roughness)/3.7+12/Reynolds)-4.781)^2/((-2۰log((Roughness)/3.7+2.51۰(-2۰log((Roughness)/3.7+12/Reynolds))/Reynolds))-2۰(-2۰log((Roughness)/3.7+12/Reynolds))+4.781)))^-2)' »" +"Astronomy" + + "PeriSel", "=" + "Mathematics" "CollatzBenchmark", "=" diff --git a/doc/calc-help/constants.md b/doc/calc-help/constants.md index 4800ae16..3adbd71b 100644 --- a/doc/calc-help/constants.md +++ b/doc/calc-help/constants.md @@ -1439,14 +1439,35 @@ Moon argument of perigee Measured. Moon's argument of perigee. It is the most recent point in its orbit when it was closest to the Earth. [Reference 25](#reference-25) -### T₀☽ constant +### Porb☽ constant + +Moon orbital period -Moon perigee passage +The Moon's anomalistic month — the mean interval between successive perigee +passages (about 27.55 days). It is distinct from the rotation period `Prot☽` and +from the sidereal month. Its absolute uncertainty (±1.12 days) is the measured 1σ +spread of real perigee-to-perigee intervals over 2000–2050, made large by the +solar perturbation; this variability is why the Moon differs from the planets. +The value is the standard anomalistic month used by `ⓁPeriSel`. [Materials 22](#materials-22) + +### T₀☽ constant -Calculated from measurement (Meeus ch.50, ±1 hour). See MPERC.txt. Moon's -most recent time of perigee passage (closest approach to Earth). The perigee -recurs approximately every 27.55 days. Computed using the Meeus algorithm -(15 terms, ±1 hour accuracy). Value in JDN. [Materials 20](#materials-20) [Materials 22](#materials-22) +Moon time of perigee passage + +Computed, not stored: `T₀☽` evaluates `JDN→(ⓁPeriSel(→Num(JDN(DateTime))))`, +calling the library routine `ⓁPeriSel` on the current date to obtain the nearest +lunar perigee, then converting the Julian Day Number to a date. `ⓁPeriSel` +implements Meeus' Chapter 50 perigee formula in full — all 60 Table 50.A periodic +terms, including the time-dependent coefficients — and matches the PyMeeus +reference implementation to better than 0.001 minute over 1950–2100. The accuracy +of the perigee *instant* against the true Earth–Moon distance minimum is limited by +the Chapter 50 method itself: about 5 minutes typical (1σ), up to ~30 minutes in +rare cases; this is the basis for the ±5 minutes (`300_s`) uncertainty. `ⓁPeriSel` +accepts either a date or a Julian Day Number as input. The perigee recurs about +every 27.55 days (`Porb☽`). Validated against observed perigees: JDN 2457706.5 → +2457706.974818 (2016-Nov-14, k=224) and 2464291.5 → 2464292.421948 (2034-Nov-25, +k=463). Note: Meeus' worked example at p.357 (JDN 2447442.35) is an apogee, not a +perigee, and must not be used to check this value. [Materials 22](#materials-22) ## Mars constants diff --git a/library/PeriSel.48s b/library/PeriSel.48s new file mode 100644 index 00000000..0c41b46e --- /dev/null +++ b/library/PeriSel.48s @@ -0,0 +1,277 @@ +« +@ ============================================================ +@ PeriSel — time of lunar perigee passage (Selene), for DB48x +@ Meeus, "Astronomical Algorithms" 2nd ed., ch. 50 — full Table 50.A +@ Merged implementation: +@ - data-driven structure (matrix x array, explicit 50.A tables) : @evgaster +@ - angle-mode-independent sine via 1_° : J. Wilson +@ +@ Input : a Julian Day Number (JDN, real) OR a date (_date-unit value). +@ A date in DB48x is a `_date`-unit value (YYYYMMDD.HHMMSS_date); a JDN is a +@ plain real. JDN() on a plain-real JDN would ERROR (read as an invalid +@ YYYYMMDD date). So the input preamble tests "does the input carry a unit +@ (date)?": if yes -> JDN converts it; if a plain real -> use as-is. +@ Output: JDN of the nearest perigee (Dynamical Time, TD/TT) +@ Pure function of its input — does NOT read DateTime. +@ +@ Accuracy of the perigee INSTANT vs the true Earth-Moon distance minimum is a +@ limit of the ch.50 method itself: ~5 min typical (median ~3 min), up to ~30 min +@ in rare cases (measured against Meeus ch.47). This is the basis for the ±5 min +@ (300_s) absolute uncertainty carried by the T₀☽ constant. The formula itself +@ matches PyMeeus (full Meeus) to < 0.001 min over 1950-2100. +@ +@ Two changes vs @evgaster's original are tagged "@ MERGED:" below: +@ 1) Each sine argument is scaled by 1_° so sin is mode-independent. +@ 2) The global Degrees-mode save/switch/restore is therefore removed. +@ JDE-row T^3 coefficient = -0.000001098, the value published in Meeus eq.50.1 +@ (confirmed in the book; matches @evgaster's original matrix). NOTE: PyMeeus +@ uses +0.000001098 instead; the resulting difference is < 0.1 s over 1950-2100 +@ (at most 2*0.000001098*T^3 d, i.e. < 0.2 s even a century from epoch 2000), far +@ below the ~5 min method accuracy, so it does not affect any perigee instant. +@ +@ Epoch (k = 0 mean perigee): JDE 2451534.6698 ; anomalistic month 27.55454989 d +@ +@ VALIDATION (full formula vs PyMeeus and real perigees): +@ 2457706.5 PeriSel -> 2457706.974818 2016-Nov-14 k = 224 +@ 2464291.5 PeriSel -> 2464292.421948 2034-Nov-25 k = 463 +@ ============================================================ + + @ Step 0: input preamble — accept a date OR a JDN (PR1701-conformance §4). + @ A DB48x date is a `_date`-unit value; a JDN is a plain real. JDN() on a + @ plain real ERRORS, so guard: if the input carries a unit (date), convert + @ it with JDN; a plain real passes through unchanged. ( date|JDN -- JDN ) + Duplicate TypeName "unit" == « JDN » IFT + + @ Step 1: k = nearest integer anomalistic-month count since the k=0 epoch. + @ ( JDN -- k ) + 2451534.6698 - 27.55454989 ÷ 0 Round + + @ Step 2: T = k / 1325.55 (Julian centuries scale used by Table 50.A). + @ ( k -- k T ) + Duplicate 1325.55 ÷ + + → k T + « + + @ Step 3: D, M, F (deg) and JDE (mean perigee) via matrix x array. + @ Each matrix row dotted with [ 1 , k , T^2 , T^3 , T^4 ]. + @ Rows: D ; M ; F ; JDE. + @ JDE-row T^3 coefficient -0.000001098 = the published Meeus eq.50.1 value + @ (PyMeeus uses +; the difference is < 0.1 s over 1950-2100 — see header). + [ + [ 171.9179 335.9106046 -0.0100383 -0.00001156 0.000000055 ] + [ 347.3477 27.1577721 -0.0008130 -0.0000010 0 ] + [ 316.6109 364.5287911 -0.0125053 -0.0000148 0 ] + [ 2451534.6698 27.55454989 -0.0006691 -0.000001098 0.0000000052 ] + ] + + @ The multiplier array [ 1 , k , T^2 , T^3 , T^4 ]. + 1 k T 2 ^ T 3 ^ T 4 ^ 5 ToArray + + @ Matrix x array -> the four results D, M, F, JDE. + @ ( Matrix Array -- D M F JDE ) + * FromArray Drop + + → D M F JDE + « + + @ Reduce the angles mod 360 to keep sin accurate, then store as globals. + D 360 mod 'D' Store + M 360 mod 'M' Store + F 360 mod 'F' Store + + + @ Step 4: Table 50.A "Argument of sine", as { Dfactor Mfactor Ffactor }. + @ Leading 0's kept; trailing 0's added programmatically below. + { + { 2 } + { 4 } + { 6 } + { 8 } + { 2 -1 } + { 0 1 } + { 10 } + { 4 -1 } + { 6 -1 } + { 12 } + { 1 } + { 8 -1 } + { 14 } + { 0 0 2 } + { 3 } + { 10 -1 } + { 16 } + { 12 -1 } + { 5 } + { 2 0 2 } + { 18 } + { 14 -1 } + { 7 } + { 2 1 } + { 20 } + { 1 1 } + { 16 -1 } + { 4 1 } + { 9 } + { 4 0 2 } + { 2 -2 } + { 4 -2 } + { 6 -2 } + { 22 } + { 18 -1 } + { 6 1 } + { 11 } + { 8 1 } + { 4 0 -2 } + { 6 0 2 } + { 3 1 } + { 5 1 } + { 13 } + { 20 -1 } + { 3 2 } + { 4 -2 2 } + { 1 2 } + { 22 -1 } + { 0 0 4 } + { 6 0 -2 } + { 2 1 -2 } + { 0 2 } + { 0 -1 2 } + { 2 0 4 } + { 0 -2 2 } + { 2 2 -2 } + { 24 } + { 4 0 -4 } + { 2 2 } + { 1 -1 } + } + + @ Pad each argument list to 3 elements with trailing 0's. + « + Duplicate Size + Case + 1 When Drop { 0 0 } + End + 2 When Drop { 0 } + End + 3 When Drop End + End + » + Map + + + @ For each padded { dD dM dF }, evaluate the Table 50.A sine term + @ directly: sin((dD·D + dM·M + dF·F) · 1_°). D, M and F are read from + @ the enclosing locals at run time. + @ NOTE: this replaces @evgaster's ToProgram-built mapper. A program + @ assembled at run time with ToProgram cannot be evaluated when PeriSel + @ is invoked inside an algebraic (e.g. the T₀☽ constant: it stays the + @ symbolic 'PeriSel(x)' and ⒸT₀☽ →Num errors). Expressing the mapper as + @ a fixed program avoids ToProgram while giving identical results. + @ MERGED: the argument is scaled by 1_° so sin is angle-mode-independent + @ (RAD/DEG/GRAD all agree; no global Degrees-mode save/switch/restore). + « FromList Drop → dD dM dF + « dD D * dM M * + dF F * + 1_° * sin » » + Map + + + @ Step 5: Table 50.A "Coefficient", as { coefficient Tcoefficient }. + @ Trailing 0 (no T-term) added programmatically below. + { + { -1.6769 } + { +0.4589 } + { -0.1856 } + { +0.0883 } + { -0.0773 +0.00019 } + { +0.0502 -0.00013 } + { -0.0460 } + { +0.0422 -0.00011 } + { -0.0256 } + { +0.0253 } + { +0.0237 } + { +0.0162 } + { -0.0145 } + { +0.0129 } + { -0.0112 } + { -0.0104 } + { +0.0086 } + { +0.0069 } + { +0.0066 } + { -0.0053 } + { -0.0052 } + { -0.0046 } + { -0.0041 } + { +0.0040 } + { +0.0032 } + { -0.0032 } + { +0.0031 } + { -0.0029 } + { +0.0027 } + { +0.0027 } + { -0.0027 } + { +0.0024 } + { -0.0021 } + { -0.0021 } + { -0.0021 } + { +0.0019 } + { -0.0018 } + { -0.0014 } + { -0.0014 } + { -0.0014 } + { +0.0014 } + { -0.0014 } + { +0.0013 } + { +0.0013 } + { +0.0011 } + { -0.0011 } + { -0.0010 } + { -0.0009 } + { -0.0008 } + { +0.0008 } + { +0.0008 } + { +0.0007 } + { +0.0007 } + { +0.0007 } + { -0.0006 } + { -0.0006 } + { +0.0006 } + { +0.0005 } + { +0.0005 } + { -0.0004 } + } + + @ Add a 0 T-coefficient where it is missing. + « + If + Duplicate Size 1 = + Then + 0 + + End + » + Map + + + @ Form each coefficient = coefficient + Tcoefficient * T. + « + FromList Drop + T * + + + » + Map + + + @ Step 6: sum JDE + Sigma( coefficient * sine ) on the fly. + → s c + « + JDE + 1 s Size + For i + c i Get + s i Get + * + + + Next + » + + + » + » +» diff --git a/src/constants.cc b/src/constants.cc index a3d5d122..14ad10e0 100644 --- a/src/constants.cc +++ b/src/constants.cc @@ -1471,6 +1471,10 @@ static const cstring basic_constants[] = " 0.1_s " " 'ROUND(ⓈProt☽/ⒸProt☽;-2)' " " 2.361E6_s ]", + // *Moon orbital period - anomalistic month (perigee to perigee) [M22] + "Porb☽", "[ 2380713.11_s " + " 97165_s " + " 0 ]", // *Moon axial tilt - Measurement [22] "ϵ☽", "[ 1.5424_° " " 0.0001_° " @@ -1501,9 +1505,10 @@ static const cstring basic_constants[] = " 0.01_° " " 'ROUND(Ⓢω☽/Ⓒω☽;-2)' " " 318.2_° ]", - // *Moon perigee passage - Calculation from measurement [M20] [M22] - "T₀☽", "[ 20260517.5833_date " - " 0_date " + // *Moon time of perigee passage - PeriSel (full Meeus ch.50); unc = ch.50 + // method accuracy vs true perigee ~5 min (1σ), up to ~30 min rare [M22] + "T₀☽", "[ 'JDN→(ⓁPeriSel(→Num(JDN(DateTime))))' " + " 300_s " " 0 ]", "Astronomy/Mars", nullptr,