OcalEN Methodology and Specifications

1. Overview and scope

OcalEN Methodology and Specifications documents the server-side, web-based lunar ephemeris calculator behind OcalEN, exposing a complete and reproducible computational pipeline through a public web interface. Given a date, time, and observer location, it returns the geocentric and topocentric position of the Moon (right ascension, declination, ecliptic longitude and latitude, altitude and azimuth), phase, illuminated fraction, age, magnitude, angular diameter, optical and physical libration, selenographic colongitude, sub-Earth and sub-solar coordinates, ICRS state vectors, and a full uncertainty budget. Auxiliary outputs include eclipse imminence parameters (Besselian elements), lunar occultation candidates against three star catalogues (Hipparcos, Tycho-2, Gaia DR3), apsides, super-Moon and micro-Moon flags, and Lunar Laser Ranging (LLR) ping distances against the five operational retroreflectors.

Scope and intended use. OcalEN is built for educational exposure of frame mechanics, time-scale chains and post-Newtonian corrections; for sophisticated amateur observation planning (astrophotography, occultation prediction, eclipse logistics); for the teaching of computational astronomy; and for parametric exploration of relativistic and gravitational effects on the Earth-Moon system. The methodology is open and reproducible from public papers and IAU/IERS conventions; every call emits a deterministic Reproducibility ID, a versioned JSON schema, an uncertainty bookkeeping that follows JCGM 100:2008 (GUM), and a CCSDS 502.0-B-3 OEM export.

What OcalEN is not. It is not real-time IERS-tracked: Earth Orientation Parameters (DUT1, polar motion, Delta-AT) come from a periodically refreshed local snapshot of IERS Bulletin A, with a documented maximum age of 365 days outside which the engine refuses to evaluate. The full list of explicit limitations is given in Section 12.7.

The algorithmic foundation is layered. The default mode combines a compact ELP-2000/82B truncation for the Moon (Chapront-Touze and Chapront 1988) with a compact VSOP87D truncation for the Sun and inner planets (Bretagnon and Francou 1988); time scales follow the standard chain (UTC, UT1, TAI, TT, TDB) using a periodically refreshed IERS Bulletin A data file for DUT1 and leap seconds. Frame transformations use Frame Bias x Precession (IAU 2006 P03; Capitaine, Wallace and Chapront 2003) x Nutation (IAU 2000A; full series in high-precision mode, otherwise a compact subset retaining the dominant signal power). In high-precision mode the pipeline upgrades to the CIO-based GCRS, CIRS, ITRS sequence and adds Schwarzschild light deflection by the Sun, Earth self-deflection, and Lense-Thirring frame dragging. A JPL DE440 path (Park et al. 2021) is provided as a higher-precision option over its native coverage interval (1550-2650).

Position accuracy in the default mode is dominated by the ELP-2000/82B truncation residual, at the few-arcsecond level over 1900-2100. The high-precision mode (DE440 plus full IAU 2000A nutation) has a design target of sub-arcsecond apparent direction; the current measured residual against the limited Horizons smoke fixtures is at the few-arcsecond / few-kilometre level for typical mid-latitude epochs (Section 12.8 reports the validation status in detail). When the IERS UT1-UTC snapshot is one year old, that age contributes up to about 5 mas to the timing-driven component of the residual. Outside the IERS coverage window the dominant residual becomes the Espenak-Meeus extrapolated Delta-T model, which reaches one second of timing uncertainty at the boundaries of its calibration interval. Detailed propagation and a per-quantity budget are presented in Section 11; the public tolerance thresholds the engine targets are documented in Section 11.10, and the validation status against JPL Horizons and other baselines is documented in Section 12.8.

The document is organised so that each section aligns with a stage of the computational graph and ends with the bibliographic anchors used by the implementation. The section numbering is contiguous and free of placeholders: every numbered section is directly implemented in the engine, and every implemented behaviour is documented under its corresponding number. Sections 4 to 9 form the algorithmic core (ephemerides, time, frames, topocentrics, relativity, observability). Section 10 enumerates the JSON contract emitted by the calculator. Section 11 presents the GUM-conformant uncertainty budget and the CCSDS 502.0-B-3 export. Section 12 documents reproducibility primitives and the regression baseline against JPL Horizons; the limitations subsection sec-limitations is embedded in 12.7.

The intended use cases divide into four groups. Pedagogical exposure of frame mechanics, time-scale chains, and post-Newtonian corrections, in a setting where the underlying formulas and references are visible alongside the numerical output. Observational planning at the precision typical of mid-aperture amateur telescopes and small institutional setups, where sub-arcsecond pointing is helpful but sub-mas tracking is not required. Methodology cross-checks using the deterministic Reproducibility ID, the explicit declaration of model versions and the EOP snapshot timestamp. Citation in academic work via the deterministic Reproducibility ID and the BibTeX block emitted by the tool.

Localisation status. EN canonicals for the lunar calendar (year and month modes), the daily moon-phase page, and the live moon panel are not yet deployed: only the PT canonical paths exist for those surfaces. Tools available in both PT and EN today are the popular calculator (/calculadora-lunar ↔ /lunar-calculator), the scientific calculator (/calculadora-lunar-cientifica ↔ /scientific-lunar-calculator), and the astrophotography calculator (/calculadora-lunar-astrofotografia ↔ /lunar-astrophotography-calculator); these expose reciprocal hreflang in the sitemap. The methodology paper itself (this document) is published in EN.

2. Computational architecture

The pipeline is organised as a directed acyclic graph of well-typed stages. Each stage consumes a typed input (epoch, observer geodetic position, model selection flags, sandbox overrides) and emits a typed intermediate object that is consumed downstream. A single engine entry point performs the full traversal; the same entry point is invoked from the user-facing form, from the JSON API, and from the regression test suite, guaranteeing a single source of truth across the three entry points. Intermediate quantities are carried in immutable, well-typed value objects so that post-processing layers (apparent place, refraction, observability synthesis) cannot silently mutate upstream state.

At a high level, the traversal covers: input validation and time-scale resolution; ephemeris evaluation in the appropriate barycentric or geocentric frame; frame transformation to GCRS and to the requested intermediate frame; light-time iteration with optional Shapiro coupling; the three-component aberration; topocentric reduction with WGS84 station coordinates and atmospheric refraction; observability synthesis (phase, libration, eclipse, occultation, LLR ping, apsides); and uncertainty quantification with optional CCSDS export. Each stage is purely functional and free of global state; intermediate caches are content-addressable and explicit; failure modes surface as typed exceptions with sufficient context for the regression suite to attach a fixture.

The same canonical phase-instant root finder used inside the scientific calculator (Section 9, principal-phase computation via bisection on the Sun-Moon elongation crossings of 0°, 90°, 180°, 270°) also drives the public calendário lunar listings rendered across the rest of the site, ensuring minute-level accuracy and a single source of truth for principal-phase timestamps across all surfaces. The historical mean-cycle approximation, which assumed a constant 29.530589-day spacing between the four principal phases of a single synodic cycle and could drift by up to roughly fifteen hours at the quarters, is no longer in the active code path.

2.1 Modules

The implementation is organised as a flat set of PHP modules, each independently unit-testable and dependency-injectable; the engine wires them together at construction time and exposes a single high-level entry point. The functional roles of the active modules at the version corresponding to this document are listed below; concrete class names, file paths, and internal call sequences are intentionally not enumerated here, as they form part of the implementation detail rather than the public methodology.

• Engine orchestrator. Top-level coordinator that hosts the stage graph, the per-call cache, and the provenance recorder.
• Lunar and solar theories. Compact analytical models for the Moon (ELP-2000/82B truncated) and the Sun (VSOP87D truncated to a compact term subset for the heliocentric position of the Earth).
• Numerical-ephemeris reader. Reader for the high-precision JPL DE440 path, plus a reader for the corresponding lunar physical-orientation product, both consumed through a uniform interface.
• Time-scale module. Pure-functional resolution of TT, UT1, UTC, TAI, TDB, with TDB-TT from a reduced Fairhead-Bretagnon series and Delta-T extended outside the IERS window via Espenak-Meeus 2006.
• Frame-chain stage. Implementation of the canonical IAU 2006 / IERS 2010 transformations (frame bias, precession, nutation, polar motion, Earth Rotation Angle, CIO locator) consistent with the SOFA reference functions cited in Sections 6 and 7.
• Earth-orientation layer. Consumer for the IERS Bulletin A data product, exposing DUT1, polar motion (x_p, y_p), celestial pole offsets (dX, dY) and leap seconds, with an extrapolation policy when the requested epoch lies outside the data window.
• Relativistic stack. Optional post-Newtonian corrections (Shapiro, de Sitter, Lense-Thirring, gravitational redshift, Wigner rotation) plus a residual diagnostic block; off by default and decorated with a sandbox flag in the output when active.
• Uncertainty-budget module. Implementation of the JCGM 100/101/102 framework: sensitivity Jacobian, Welch-Satterthwaite degrees of freedom, Monte Carlo sampler, and CCSDS 502.0-B-3 OEM export.
• Observability layer. Phase, libration, eclipse Besselian elements, lunar occultations against three star catalogues (Hipparcos, Tycho-2, Gaia DR3), apsides and super/micro Moon flagging, and LLR retroreflector ping path.
• Data-refresh layer. Maintenance-key-protected refresh of the IERS Bulletin A data product and the auxiliary Earth and lunar physical-orientation files, behind a small administrative endpoint.

2.2 Caching

The engine separates three levels of caching. Per-call memoisation reuses previously computed nutation matrices and IERS data interpolations within a single evaluation; the cache key is the request epoch in TT, with model selection flags participating in equality checks. Per-process caching exploits the PHP opcache for the constant analytical-series tables, avoiding a multi-megabyte constant reload on every request when the high-precision mode is active. Per-deployment caching takes the form of a periodically refreshed local copy of the IERS Bulletin A data product; its age is exposed in the provenance block and a warning is emitted when it exceeds the configured threshold (default 30 days, hard limit 365 days).

The numerical-ephemeris reader maintains its own LRU cache of recently decoded Chebyshev coefficient records to amortise disk I/O across closely spaced epoch queries. The default cache holds 64 records, which corresponds to roughly 5.6 years of continuous coverage at the standard 32-day record interval used by the bundled high-precision path; the size is configurable so that batches spanning longer windows (multi-year monthly tables, multi-decade occultation surveys) avoid eviction thrashing.

The binary orientation reader maintains an analogous record-level LRU cache to amortise disk I/O in tight loops (libration heatmap, libration compass tile). The default cache holds 16 records, configurable via constructor parameter; this is sufficient for the sub-day Chebyshev sub-intervals used by the IAU lunar and high-precision Earth orientation products without thrashing across a one-month grid. BE-IEEE kernels are explicitly rejected with an actionable error pointing to NAIF toxfr conversion, defensively consistent with the numerical-ephemeris reader.

2.3 Provenance and reproducibility primitives

Every evaluation produces a provenance sub-object that records the engine identifier, the engine version, the active model selection (precession, nutation, frame mode, refraction, advanced GR toggles), the IERS data identifier and age, the active sandbox overrides (if any), and the deterministic Reproducibility ID introduced in Section 12.2. Two requests that produce identical Reproducibility IDs return byte-for-byte identical numerical outputs, modulo non-determinism in concurrent IERS data refreshes which is ruled out by the same-version constraint. The provenance is the canonical means by which a result is cited in academic work and the canonical artefact compared against in the regression suite.

2.4 Stage graph

[Inputs] -> validate -> resolve_time_scales
                       -> ephemeris (BCRS)
                       -> frame_chain (GCRS, CIRS, ITRS)
                       -> light_time_iteration (+ Shapiro coupling)
                       -> aberration (annual + diurnal + galactic)
                       -> topocentric_reduction (WGS84 + refraction)
                       -> observability_synthesis (phase, eclipse, occultation, ...)
                       -> uncertainty_quantification (+ CCSDS export)
                       -> [Outputs]

Each arrow is a pure function; each node has a unique identifier consumed by the regression-fixture matcher. The graph is non-recursive (no cycle), so the cost of a single call is the sum of stage costs; in practice the dominant terms are the full IAU 2000A nutation series in the high-precision mode and the Chebyshev evaluation of the DE440 numerical ephemeris when that mode is selected, both well below 100 ms per request on a commodity shared-hosting environment.

3. Inputs accepted

This section describes the inputs the form accepts, grouped by physical meaning rather than by URL parameter. The visible form labels are the canonical identifiers; each input has a documented default chosen so that an untouched form returns results consistent with the current IERS bulletins and the IAU 2006 / 2000A resolutions. New inputs are added without breaking existing permalinks (Section 12.4); valid ranges for date and time are bounded by the underlying ephemeris kernel coverage and by the time-scale tabulations adopted in Section 5. Detailed enumeration of the option values exposed by each selector is reserved to the form itself; here we describe what each input is, what physical quantity it controls, and the standard that defines it.

3.1 Date and time

The required temporal input is a calendar date plus a time of day in UTC. The user can equivalently supply a Julian Date, a Modified Julian Date, or a Julian Date already expressed in TT, UT1 or TAI for the advanced time-scale stack of 3.4; when more than one time-scale field is filled, the engine resolves them by precedence rules and stamps the active branch in the provenance block. The calendar selector toggles between the Gregorian and Julian conventions for dates that span the 1582 reform; ISO 8601 governs the surface formats. The display time zone follows the IANA zoneinfo identifier supplied by the form and affects only the human-readable layer of the result, not the numerical computation. The valid epoch range is bounded by the underlying ephemeris kernel coverage (Section 4.1) and by the Espenak-Meeus 2006 Delta-T extrapolation window adopted in Section 5.4.

3.2 Observer location and atmosphere

The observer panel collects the geodetic latitude, geodetic longitude, ellipsoidal height, and IANA time zone for the station, plus a binary observer-mode selector that switches between a topocentric reduction (which includes the parallactic shift induced by the observer's position on the rotating Earth) and a geocentric reduction (which suppresses it; see Meeus 1998 Ch. 47). Latitude and longitude follow the WGS84 conventions for sign and origin (positive east of Greenwich, positive north of the equator); altitude is the height above the WGS84 ellipsoid in metres. The form populates the location panel from the visitor's IP-derived geolocation, falling back to a reference site at Brasilia when geolocation is unavailable.

The atmosphere subsection collects the meteorological state of the column above the observer: pressure, temperature, relative humidity, the reference wavelength used for chromatic refraction, the temperature lapse rate, the integrated ozone column, the precipitable water vapour column, and the aerosol optical depth at the 550 nm reference wavelength. Defaults follow ICAO Standard Atmosphere conventions; the chromatic-refraction reference wavelength corresponds to the photopic peak of human vision and matches the value adopted by the Astronomical Almanac for visual observation. The CO₂ mixing ratio used by the Ciddor (1996) refractive-index formula and the cloud-cover fraction reported in the METAR convention are exposed alongside the basic state and feed the higher-precision refraction path of 3.3.

3.3 Engine and frame model selection

This subsection describes the model selectors that determine which physical theory drives each layer of the apparent-place reduction. The default selection follows the IAU 2006 / IAU 2000A resolutions and the IERS Conventions 2010 and yields apparent topocentric coordinates suitable for visual and CCD observation; alternate choices are supported for cross-validation against pre-2003 publications and against the current institutional ephemeris services. The chosen path is stamped in the provenance block (Section 10.6) so the result is reproducible bit-for-bit at any future date.

Time scale engine selection. The user can request the lite analytical engine (a compact-truncation evaluation of the published ELP-2000/82B and VSOP87D series), the high-precision DE440 path (Park et al. 2021), or let the engine auto-select between the two based on the requested epoch and the kernel availability detected at boot time. The auto path falls back to the analytical layer outside the DE440 native coverage and stamps the active branch in the provenance block.

Apparent versus geometric reduction. The apparent path applies the full chain of stellar aberration (Section 7.3), nutation (Section 6.3), and light-time correction (Section 7.2); the geometric path returns the purely geometric direction at the requested epoch and matches the corresponding JPL Horizons quantity for cross-validation.

Nutation model selection. The default follows the IAU 2000A series (Mathews et al. 2002), which retains the dominant signal power across the luni-solar and planetary contributions. The user can substitute the abridged IAU 2000B subset of McCarthy and Luzum (2003) for faster evaluation at the milliarcsecond level, or the equinox-based IAU 1980 series (Wahr 1981) for compatibility with pre-2003 publications.

Precession model selection. The default follows the IAU 2006 P03 series (Capitaine, Wallace and Chapront 2003), which is the IAU-recommended replacement for the IAU 1976 series (Lieske 1977); both are exposed for cross-validation. The frame-realisation selector chooses between the GCRS aligned with the ICRS through the IAU 2006 frame bias and the legacy ICRS realisation that predates the IAU 2006 resolutions, again as a reproducibility marker rather than a numerical-result-changing toggle at the lunar distance.

Frame origin selection. The frame-origin selector chooses between the equinox-based pipeline (the classical nutation-precession-bias matrix chain) and the CIO-based pipeline (the X, Y, s parametrisation introduced by the IAU 2006 resolutions), both yielding numerically equivalent apparent coordinates at the milliarcsecond level when the underlying nutation and precession models are matched. The CIO-based pipeline is the IAU-recommended path; the equinox-based pipeline is offered for compatibility with publications that adopt the older convention.

Aberration model selection. The user can apply annual aberration alone (the dominant 20.5 arcsec term induced by Earth's heliocentric motion) or annual plus diurnal aberration (which adds the small correction induced by the observer's daily rotation around the geocentre, of magnitude 0.32 cos(latitude) arcsec). Galactic aberration is exposed as a separate frame-correction flag because it operates on a different physical scale and is disabled by default per the IAU 2006 Resolution B1 prescription.

Atmosphere and refraction model selection. The user can request no refraction, the Bennett (1982) closed-form approximation suitable for general visual observation, the Saemundsson (1986) variant which reduces residual error at low elevation, or the Mendes-Pavlis (2004) FCULa mapping function adopted by the IERS Combination Centre for Lunar Laser Ranging and VLBI data reduction. The chosen model is fed by the meteorological state of 3.2.

Relativistic corrections selection. The default applies the Newtonian light-time path; the user can extend the reduction with the Shapiro time delay (Section 8.1, the logarithmic correction induced by the Earth and Sun gravitational potentials) or the Lense-Thirring frame dragging (Section 8.3, the gravitomagnetic Earth-spin contribution that drives the lunar nodal precession). A master toggle activates the full Section 8 stack including the de Sitter geodetic rotation, the gravitational redshift on the lunar coordinate time, and the Wigner rotation on the relativistic Doppler.

Observability filter set. Ten parameters covering elevation cutoff, airmass, solar elongation bounds, hour-angle cutoff, daylight skip, magnitude limit, time-digit resolution, extra-precision flag, and CSV output formatting are exposed in a dedicated fieldset; see 3.9 for the full description.

Occultation catalog selection. The user selects the stellar catalogue consulted by the occultation-search engine (Section 9.3): Hipparcos as the default bright-star catalogue, Tycho-2 for fainter visual magnitudes (with a manual install step because of the catalogue size), or the Gaia DR3 subset bounded by Gmag for the highest astrometric precision available (also a manual install).

Uncertainty propagation options. Two selectors govern the uncertainty budget of Section 11. The Monte Carlo trial-count selector activates the JCGM 101:2008 cross-validation path with 10000, 100000, or 1000000 trials in addition to the GUM-linear default; trial counts above 100000 add about one to ten seconds of latency. The coverage-factor selector chooses between the one-sigma, two-sigma, and three-sigma intervals reported next to the apparent place; the JCGM 100:2008 Section 6.3 default is the two-sigma value, which corresponds to approximately 95 percent coverage in the Gaussian limit.

3.4 Earth orientation and precision-time stack

The Earth Orientation Parameters describe the deviation of Earth's rotation from a uniform inertial reference; the precision-time stack defines the conversion path between the user-facing UTC instant and the dynamical time argument consumed by the ephemeris. The advanced fieldset exposes optional overrides for the canonical EOP quantities (UT1 minus UTC, polar motion, length of day, the celestial pole offsets dX and dY, and the integer leap-second offset between TAI and UTC) and for the dynamical time-scale stack (Delta-T, the TT to TDB periodic series, and the multi-time-scale output set TT, TAI, UTC, UT1, TDB, TCG, TCB, TCL, GPS as defined by IAU 2000 Resolution B1.7). When the user supplies any of these values explicitly, the engine bypasses the corresponding IERS Bulletin A snapshot for that single call and stamps the manual override into the provenance block. This is the recommended path for epochs far from the present: a researcher reproducing a 1950 occultation from the original Hipparcos circular can supply the contemporary IERS values and obtain a result that matches the original publication.

The free oscillation of Earth's fluid core relative to the mantle, known as the Free Core Nutation, is exposed as an amplitude-and-phase pair following the Lambert (2007) empirical model; this oscillation is not predicted by precession-nutation models and must be supplied empirically when sub-milliarcsecond precision is required. The TT to TDB conversion accounts for the difference between Terrestrial Time (geocentric) and Barycentric Dynamical Time (solar-system barycentric); the user can choose between the Fairhead-Bretagnon 1990 series at the nanosecond level, the same series augmented by Irwin and Fukushima (1999) for sub-nanosecond precision, and the JPL TE405 series derived directly from DE405. The planetary ephemeris selector exposes DE440, DE441, INPOP21a (IMCCE) and EPM2021 (IAA RAS) as declarative reproducibility markers; the native engine path runs DE440 (Park et al. 2021) and the alternate identifiers are recorded in the JSON snapshot without altering the underlying numerical computation, since those datasets are not bundled locally.

Defaults follow the current IERS Bulletin A snapshot (a value near zero for UT1 minus UTC and for polar motion at the present epoch, the integer leap-second count valid for 2026), the canonical IERS Conventions 2010 prescriptions for the celestial pole offsets and the length-of-day correction, and the Lambert (2007) reference values for the Free Core Nutation amplitude and phase. Valid ranges for each quantity are bounded by the documented physical envelope of the corresponding IERS product (typically below one second for UT1 minus UTC, below one arcsecond for polar motion, below one milliarcsecond for the celestial pole offsets); the form rejects values outside this envelope rather than silently clamping them.

3.5 Station displacement layer (IERS 2010 Chapter 7)

This subsection describes the optional station-displacement layer that the engine applies to the WGS84 station position before the topocentric reduction. The IERS Conventions 2010 Chapter 7 specifies six categories of displacement that should be modelled for high-precision geodetic and astrometric work: solid-Earth tides driven by the lunisolar tidal potential, ocean tide loading driven by the redistribution of ocean water under the same potential, atmospheric pressure loading driven by inverse-barometer response, hydrological loading driven by continental water storage variations from GRACE, the pole tide driven by the centrifugal redistribution of polar motion, and the deflection of the vertical decomposed into a north-south component xi and an east-west component eta. The dominant contribution at typical mid-latitude sites is the solid-Earth tide, which reaches approximately 30 cm in peak-to-peak vertical displacement at the M2 maximum (Petit and Luzum 2010 Eq. 7.5); ocean tide loading from the FES2014b global hydrodynamic model (Lyard et al. 2021) reaches roughly 5 cm at coastal stations and is well below 1 cm in continental interiors. When the user activates a category, the engine convolves the corresponding Green's function or empirical model with the WGS84 station position and stamps the active layer in the provenance block; when the category is left at its default off state, the residual is folded into the uncertainty budget of Section 11.

3.6 Modern reference systems

This subsection describes the selectors that lock the celestial-to-terrestrial-to-selenocentric chain to the current canonical realisations. The celestial frame selector chooses between the ICRF3 realisation (Charlot et al. 2020), the legacy ICRF2 realisation, and the GCRS aligned with the ICRS through the IAU 2006 frame bias; the terrestrial frame selector chooses between the ITRF 2020 solution, the ITRF 2014 solution, and the WGS84 G2139 epoch realisation; the time-scale realisation selector chooses between the BIPM-published TT realisation derived from primary frequency standards and the TAI-derived realisation. The selenocentric frame selector chooses between the Mean Earth / Polar Axis realisation distributed with DE440 and the dynamical Principal Axis realisation recommended by the IAU 2009 Working Group on Cartographic Coordinates and Rotational Elements (Archinal et al. 2018); the two realisations differ by a small rotation whose accumulated surface offset is approximately 100 m, documented in JPL Interoffice Memorandum 343R-08-003 (Williams et al. 2014).

These selectors are useful for cross-validation against external products published in a specific realisation: the JPL Horizons output prior to 2018 is referenced to the legacy ICRF2 realisation, while post-2018 output is referenced to ICRF3. The numerical impact of switching between ICRF2 and ICRF3 at the lunar distance is well inside the analytical-engine residual; the choice is recorded in the provenance block as a reproducibility marker rather than as a result-changing toggle. The reporting frame selector, the geodetic earth model selector (WGS84, GRS80, or the IERS 2010 reference ellipsoid), and the geoid model selector (EGM2008, EGM96, or the no-geoid option that reports purely ellipsoidal heights) sit alongside the frame selectors and are consumed by the corresponding output formatting layer.

3.7 Sandbox constant overrides

This subsection describes the sandbox tier exposed for parametric sensitivity studies, partial-derivative cross-checks against published Jacobians, and pedagogical exploration of how individual fundamental constants enter the lunar ephemeris. The user can rewrite, on a per-call basis, the heliocentric, geocentric and selenocentric gravitational parameters fixed by IAU 2015 Resolution B3 and by the DE440 fit (Park et al. 2021); the speed of light defined by the 2019 redefinition of the SI base units; the Newtonian gravitational constant from the CODATA 2018 recommended values; and the Earth zonal harmonics J₂, J₃, J₄ from EGM2008 (Pavlis et al. 2012). Every overridden constant is flagged as sandbox in the JSON output, and the affected result lines are decorated so the caller cannot mistake them for a baseline result.

The same fieldset exposes the parameters of the theoretical-lab tier: the Nordtvedt parameter eta (a strong-equivalence-principle test via Lunar Laser Ranging following Nordtvedt 1968), a cosmological G-dot over G drift (the Williams, Turyshev and Boggs 2012 LLR upper limit is at the 10^-13 per year level), a non-zero graviton mass through a Yukawa-screened gravitational potential, and a local cosmological-constant term in the Schwarzschild-de Sitter metric. These inputs are offered as a transparent vehicle for exploring how individual parameters enter the apparent place; production work on the calculator uses the IAU 2015, CODATA 2018, DE440 and EGM2008 nominal values and never reads the override fields. Results carrying any non-default override are flagged in the JSON snapshot and the reproducibility hash so downstream consumers can recognise the sandbox provenance and decline to mix it with observational comparisons.

3.8 Edge-of-physics groups (Group 1, 2, 3)

The form exposes three groups of toggles that collectively wire the relativistic and edge-of-physics layer documented in Section 8. Group 1 (Standard GR Confirmed) activates Shapiro time delay, geodetic (de Sitter) precession, galactic tide residual vector, polar motion W matrix, and lunar Yarkovsky thermal recoil. Group 2 (Lunar Clock and Photonics) activates Lunar Coordinate Time (TCL) export, gravitational redshift between observer and lunar surface, and the laser-ping path through the LLR retroreflectors. Group 3 (Theoretical Lab) activates the sandbox toggles of 3.7 (Nordtvedt eta, G-dot/G, graviton mass, SME coefficients, cosmological constant). Each toggle defaults off; when enabled, the affected output line is decorated with a sandbox flag so it cannot be confused with a baseline result.

3.9 Observability filter set

This subsection describes the observability filter set, a standard set of constraints defined by the JPL Horizons User Manual that the user can apply to any computed result. The filter set covers ten axes. The elevation cutoff rejects sub-horizon epochs. The airmass maximum rejects extreme-airmass epochs near the horizon where atmospheric refraction becomes numerically unstable. The solar elongation lower and upper bounds filter epochs by Sun-target separation. The local hour angle cutoff rejects equatorial-mount pier-flip regions. The daylight skip flag rejects epochs when the Sun is above the horizon. The apparent visual magnitude maximum filters by brightness. Three further axes control output formatting (extra-precision flag, time-digit resolution, CSV formatting toggle) and have no effect on the underlying numerical computation. Default values are chosen so that an untouched form returns every computed result as passing; the user tightens individual axes to match the constraints of a specific observation programme. The filters do not modify the underlying ephemeris computation; they evaluate the computed result against the user-defined criteria and report a pass or fail status for each axis along with an aggregate result.

Airmass is computed via the Pickering 2002 approximation: X = 1 / sin(h + 244 / (165 + 47 * h^1.1)), with h in degrees of apparent altitude. Solar elongation reads from the engine's phase block. The local hour angle is taken from the topocentric reduction and normalised to the [-180, +180] degree interval before the comparison. The Sun altitude used by the daylight filter is computed from the engine's solar right ascension and declination and from local sidereal time using the standard horizontal-coordinate formula sin(h) = sin(phi) sin(dec) + cos(phi) cos(dec) cos(H), with the observer geographic latitude phi taken from the form.

3.10 Public input contract (provenance)

Every input parameter accepted by the calculator (the full set of approximately 130 axes documented in subsections 3.1 to 3.9) is reflected verbatim in the JSON output payload under the scientific.provenance.inputs path. Each entry exposes the same four fields, regardless of whether the value was supplied by the user or fell back to its documented default:

• value. The actual numerical or categorical value used by the engine for this call. After clamping, enum coercion, regex validation, or default substitution have all been applied; this is the value that drives every downstream computation.
• default_used. A Boolean flag. True when the form did not supply the parameter (or supplied an empty string), so the documented default was substituted. False when the user actively supplied the value present in the value field. The flag is independent of whether the supplied value was equal to the default; it only records the route taken through the input layer.
• unit. A short human-readable unit string (degrees, hPa, ISO 8601 date, electron volts, dimensionless, enum, boolean, and so on). This is the canonical unit; any client that re-supplies the input must use the same unit, since the engine performs no unit conversion at the input layer.
• sanitization. A short label describing the rule applied at the input layer. Common patterns are clamp [a, b] for numerical clamps with explicit bounds, enum {x, y, z} for categorical inputs with the accepted vocabulary, regex YYYY-MM-DD for ISO date validation, and is_finite for plain numerical inputs without bounds. Boolean checkbox-style inputs use enum {0, 1, on} or document the absence-implies-true pattern explicitly.

Reproducibility contract. The combination of value, unit, and the engine version stamped into scientific.provenance is sufficient to reproduce a result bit-for-bit. A third party recovers the call by re-supplying every value to the corresponding query string parameter, paying attention only to the unit string (no implicit assumptions about preferred units, default coordinate frames, or default model selectors). The default_used flag is informational; it lets a citation indicate when a value was the documented default versus when it was a user-specified override, which matters for academic transparency.

The contract definition is centralised in engine/InputContract.php and is the single source of truth for the unit and sanitization labels emitted by the JSON payload. Adding a new input to the form requires updating that file's definitions() method and adding the corresponding key to the $_inputContractValues array assembled by the page entry point. Removing an input has the symmetric requirement; a regression test verifies that the inventory of $_GET reads in the page entry point matches the keys in the contract definitions, so a form change cannot ship without the JSON contract being updated in the same release.

4. Theory of ephemerides

The ephemeris layer produces the position and velocity of the Moon, the Sun, and the inner planets at the requested epoch in the appropriate barycentric or geocentric frame. Three theories are used in combination, with explicit precedence rules so that the path actually taken is always recoverable from the provenance block of the JSON output.

4.1 Moon (ELP-2000/82B truncated + DE440 + lunar PCK)

The default Moon position uses the semi-analytical lunar theory ELP-2000/82B of Chapront-Touze and Chapront (1988). The default mode uses a compact truncation of the published series, expressed in the fundamental Delaunay arguments D (Moon's mean elongation from the Sun), M (Sun's mean anomaly), M' (Moon's mean anomaly), and F (Moon's argument of latitude). The truncation is chosen so that the residual against the full series stays below the apparent-direction precision targeted in the default mode, and so that the discarded short-period terms remain inside the documented uncertainty budget of Section 11. The series structure is:

lambda(t) = L0(t) + sum_i  A_i^L * sin(D_i * D + M_i * M + Mp_i * Mp + F_i * F)
beta(t)   = sum_i  A_i^B * sin(D_i * D + M_i * M + Mp_i * Mp + F_i * F)
r(t)      = a + sum_i  A_i^R * cos(D_i * D + M_i * M + Mp_i * Mp + F_i * F)

where  D, M, Mp, F  are linear in T = (JD_TT - 2451545.0) / 36525
and the integer multipliers (D_i, M_i, Mp_i, F_i) come from the
ELP-2000/82B published table (Chapront-Touze and Chapront 1988).

The expected accuracy of the truncated form is 3 to 5 arcsec in apparent angular position over the 1900-2100 window (Meeus 1998 Ch. 47), dominated by the discarded short-period terms below the truncation threshold. Outside this window the residual grows because Delta-T uncertainty propagates into the mean motion arguments through L0(t); at the Late Antique boundary (year 500 CE) the residual reaches approximately 30 arcsec, and at the Babylonian boundary (year -500) it reaches several arcminutes. The Espenak-Meeus 2006 Delta-T extrapolation is the recommended path for these epochs and is selected automatically when the requested date falls outside 1972 to present.

When the high-precision DE440 numerical ephemeris (Park, Folkner, Williams and Boggs 2021) is available locally, the engine swaps the analytical layer for a Chebyshev evaluation of the published coefficients. The reader follows the standard SPICE Type-2 (Chebyshev polynomial) interpolation conventions:

position(t) = sum_{k=0..N}  C_k(component) * T_k(tau),
              tau = 2 * (t - T_mid) / INTLEN
T_k = Chebyshev polynomial of the first kind, degree k.

Coverage spans the native 1550-2650 interval; over this window the 1-sigma position error for the Moon is sub-arcsecond at J2000 and grows toward the kernel boundaries (Park et al. 2021 Table 5). Activation is automatic when the corresponding data is present and the requested epoch falls inside the coverage; otherwise the default mode is used and the chosen engine identifier is logged in the output payload.

For lunar physical libration and orientation (selenographic frame), the engine consumes the corresponding NAIF lunar binary PCK orientation product, which carries the integrated principal-axis frame from the same DE440 dynamical model. The free-libration amplitudes of Eckhardt (1981) are recovered to sub-arcsecond level when the binary orientation product is available; otherwise the Eckhardt 1981 closed-form reduction is used and the libration residual is folded into the uncertainty budget.

The lunar phase angle and the apparent magnitude follow Meeus 1998 Ch. 48 once the Moon and Sun positions are available; the magnitude formula uses the V-band photometric correction with the Hapke phase function approximation (Hapke 1984), which is accurate to 0.05 magnitudes over the full lunation cycle.

4.2 Sun (VSOP87D)

Heliocentric ecliptic spherical coordinates of date for the Sun (the negative of Earth's heliocentric position) are computed from VSOP87D (Bretagnon and Francou 1988) in its "D" variant, which delivers the heliocentric ecliptic spherical coordinates of date directly without the dynamical-frame rotation required by the "A" variant. The default mode retains a compact term subset for the Earth's heliocentric position, with a sub-arcsecond residual at J2000 in heliocentric longitude. The series structure is:

L(t) = L0(T) + L1(T) * tau + L2(T) * tau^2 + ... + L5(T) * tau^5
B(t) = B0(T) + B1(T) * tau + ... + B5(T) * tau^5
R(t) = R0(T) + R1(T) * tau + ... + R5(T) * tau^5
tau = (JD_TT - 2451545.0) / 365250.0
each L_k, B_k, R_k = sum_i  A_ki * cos(B_ki + C_ki * tau)

The Sun's apparent geocentric position is then obtained by reversing the heliocentric vector and applying aberration (Section 7.3) and light-time correction (Section 7.2). For very long baselines (more than about 4000 years from J2000), the truncated VSOP87D residual grows to several arcseconds in solar position; the DE440 path is the recommended option for those epochs when available. The full VSOP87D series recovers the residual to the sub-mas level; the truncated form used here is a deliberate trade-off favouring server cost over photometric precision, since the Sun's position primarily affects phase and elongation and is not a sub-mas observable in this calculator.

4.3 Planets (VSOP87D + DE440)

Inner-planet positions for solar-system perturbations on the apparent place (Schwarzschild deflection by the Sun on photons emitted near Mercury or Venus, planetary aberration of the inner system, Mars and Venus in occultation searches against the lunar limb) come from VSOP87D in the default mode and from DE440 in the high-precision mode. The two paths are interchangeable to within the documented residuals (sub-arcsecond for the inner planets in the truncated VSOP87D, sub-mas for DE440); the choice is logged in the output provenance.

The outer planets (Jupiter, Saturn, Uranus, Neptune) and the four largest asteroids (Ceres, Pallas, Vesta, Hygiea) are not directly modelled in the default path; their gravitational perturbations on the Moon's orbit are absorbed into the ELP-2000/82B fitted constants. When DE440 is active, all eight planets and the largest asteroids are available for the relativistic perturbation analysis of Section 8 (residual perturbation vector); the corresponding effect on the lunar apparent place is below 10 microarcsec and is reported as a diagnostic.

5. Time scales

Five time scales are explicitly maintained throughout the pipeline (UTC, UT1, TAI, TT, TDB), connected by the standard chain documented in IERS Conventions 2010 Chapter 5 and IAU 2000 Resolution B1.7. Each time scale is exposed in the output payload as both Julian Date and ISO 8601 string, and the conversion machinery is purely functional: a single closure produced by TimeScales::closuresFor(JD_UTC) is reused throughout an evaluate() call so that the values stay self-consistent.

5.1 TT, UT1, UTC, TAI, TDB

The scalar relationships, with all values in SI seconds, are:

TT  =  TAI + 32.184 s
TAI =  UTC + Delta-AT          (Delta-AT = total leap seconds, IERS Bulletin C)
UT1 =  UTC + DUT1              (DUT1 = UT1 - UTC, |DUT1| < 0.9 s by definition)
TDB =  TT + Delta-TDB(TT)      (Fairhead-Bretagnon FB91 reduced; see 5.3)
GPS =  TAI - 19.0 s            (constant offset from 1980-01-06 Sunday epoch)

The constant 32.184 s is the historical anchor between Ephemeris Time and TAI established by IAU 1976 Resolution 5 to maintain continuity across the 1977 January 1.0 redefinition. The value of Delta-AT is updated by IERS via Bulletin C and currently sits at 37 s (no leap second has been declared since 2017 January 1.0); the snapshot used by this calculator is the value at the IERS Bulletin A snapshot date and is exposed in the provenance block as delta_at.

The relationship UT1 = UTC + DUT1 is the canonical link between atomic time and Earth-rotation time; the magnitude bound |DUT1| < 0.9 s is enforced by the leap-second insertion mechanism of IAU 2003 Resolution B1, and DUT1 is published daily by IERS in Bulletin A. For epochs before 1972 January 1.0, the UTC scale is not formally defined; the engine uses the convention DUT1 = 0 and Delta-AT = 10 in that regime, which is consistent with the contemporary practice of the US Naval Observatory and matches the Espenak-Meeus 2006 historical Delta-T extrapolation.

The engine resolves the chain in closures so that no global state is mutated. Each time-scale value is exposed in the output payload (UTC, UT1, TAI, TT, TDB, JD, MJD) and can be selected for downstream output through the output_time_scales parameter, which is a multi-valued enum accepting any subset of {TT, TAI, UTC, UT1, TDB, TCG, TCB, TCL, GPS}.

5.2 Earth Rotation Angle (ERA, IAU 2006) and CIO-based GST

The Earth Rotation Angle (ERA) is defined as the geocentric angle between the Celestial Intermediate Origin (CIO) and the Terrestrial Intermediate Origin (TIO), measured along the Celestial Intermediate Pole equator. It replaces the classical Greenwich Sidereal Time in the IAU 2006 framework as the canonical fundamental quantity for Earth rotation, linking the rotational state of the Earth with the celestial reference frame without ambiguity from precession or nutation. The defining expression is a strictly linear function of UT1 (Capitaine and Wallace 2006):

Theta(UT1) = 2 pi * (0.7790572732640 + 1.00273781191135448 * Tu)
Tu = JD(UT1) - 2451545.0

The constants in this expression encode two physical facts: 0.7790572732640 is the value of the angle (in fraction of a turn) at the J2000.0 epoch, and 1.00273781191135448 is the ratio of the mean sidereal day to the mean solar day. The expression is exact in the sense that it defines the relation between UT1 and rotation, modulo the specification of UT1 itself which is anchored to the IAU/IERS conventions. The implementation matches the SOFA reference function iauEra00 bit-for-bit on the canonical test vectors documented in the SOFA cookbook (IAU SOFA Board 2021).

Greenwich Sidereal Time in the IAU 2006 framework is recovered from ERA via the equation of the origins (EO):

GST(IAU 2006) = ERA - EO(t)
EO(t) = polynomial in T (centuries since J2000.0)

The polynomial in EO is dominated by precession; the leading term in arcseconds per century is approximately -4612.156534 (Capitaine, Wallace and Chapront 2003 Eq. 28). The full series, accurate to several microarcseconds over 1900-2100, is:

EO(arcsec) = -0.014506
           - 4612.156534 * T
           -    1.3915817 * T^2
           +    0.00000044 * T^3
           +    0.000029956 * T^4
           +    0.0000000368 * T^5
T = (JD_TT - 2451545.0) / 36525

For higher precision the nutation contribution must be added explicitly. This implementation uses the truncated polynomial form for the secular component and relies on the NPB chain to absorb the nutation residual; the resulting GST agrees with the SOFA reference function iauGst06a at the 5 microarcsec level over the full IAU 2006 calibration interval.

5.3 TDB-TT Fairhead-Bretagnon (FB91 reduced)

Barycentric Dynamical Time (TDB) differs from Terrestrial Time (TT) by a periodic series with no secular drift, dominated by the eccentric Earth orbit (annual term) and the Moon (monthly term). The TDB-TT difference is computed via the Fairhead-Bretagnon series (Fairhead and Bretagnon 1990). The implementation uses the reduced 41-term form catalogued in Hilton and Hohenkerk (2002), which retains all amplitudes above 0.1 microsecond:

TDB - TT = sum_{i=1..41}  A_i * sin(omega_i * T + phi_i)      [microseconds]
T = (JD_TT - 2451545.0) / 365250.0
omega_i, phi_i, A_i = constants from the FB91 reduced table.

The dominant term has amplitude 1.6568 ms with a one-year period (corresponding to Earth's heliocentric mean motion); the secondary term has amplitude 22.418 microseconds with a half-year period; the lunar contribution at one month period reaches 13.840 microseconds. The full 41-term sum reproduces the JPL TE405 reference at 0.1 microsecond over the 1900-2100 window (Hilton and Hohenkerk 2002 Table 2). The optional tt_tdb_model = te405 selector swaps the analytical reduced form for a direct TE405 numerical lookup, which is recommended for sub-nanosecond timing applications.

5.4 Delta-T (Espenak-Meeus extrapolation)

The quantity Delta-T = TT - UT1 measures the departure of Earth's rotational time from atomic time. Inside the IERS coverage window (1972 to current), Delta-T is recovered analytically from the chain Delta-T = 32.184 s + Delta-AT - DUT1; both Delta-AT and DUT1 come from the IERS Bulletin A snapshot. Outside the window, the Espenak-Meeus 2006 polynomial (NASA TP-2006-214141) provides the modelled Delta-T for historical reconstruction (back to year -1999) and forward extrapolation. The polynomial is piecewise on epoch:

y = year + (month - 0.5) / 12

case y < -500:        Delta-T = -20 + 32 * u^2,                u = (y - 1820) / 100
case -500 <= y < 500: Delta-T = polynomial of degree 6 in u    u = y / 100
case 500 <= y < 1600: Delta-T = polynomial of degree 6 in u    u = (y - 1000) / 100
case 1600 <= y < 1700, 1700 <= y < 1800, ..., 2005 <= y < 2050: piecewise polynomial
case y >= 2050:       Delta-T = -20 + 32 * u^2 - 0.5628 * (2150 - y)

The polynomial reproduces the IERS series to better than 1 second near the present and is the de facto educational standard. For historical eclipse work it provides a useful first approximation; for sub-second timing of pre-1900 events the Stephenson, Morrison and Hohenkerk (2016) tabulation of long-term Earth rotation is the recommended supplement, but it is not currently bundled with this calculator.

Boundary smoothing 2045-2055. The piecewise model has a documented sub-second discontinuity at year=2050 between the 2005-2050 quadratic branch (about 92.13 s at year=2049) and the 2050-2150 parabolic branch (about 93.0 s at year=2050). The engine applies a 10-year linear blend across 2045-2055 to remove the cosmetic jump. The blend amplitude is well below the documented Delta-T uncertainty in this range (about plus or minus 10 seconds post-2050, since the actual values depend on future Earth rotation which has not yet happened), so the smoothing has no science impact and is purely presentational.

Estimate vs canonical Delta-T. The time-scale bundle returns two distinct fields. The field delta_t_seconds is the Espenak-Meeus polynomial estimate (informational, model-based). The field delta_t_actual_seconds is the real Delta-T derived from the canonical chain (jd_tt - jd_ut1) * 86400, with TT obtained from the exact relation TT = TAI + 32.184 s and UT1 anchored on the IERS DUT1 snapshot. Inside the IERS window the two fields can differ by several seconds because the polynomial is a forward fit while the canonical value reflects the measured Delta-AT and DUT1 at the instant. Callers building UQ budgets or comparing against external tools should consume delta_t_actual_seconds; the polynomial estimate is retained as a baseline for historical and far-future epochs where IERS data does not exist.

5.5 IERS Bulletin A snapshots

DUT1, x_p, y_p, dX, dY and Delta-AT (leap seconds) values come from a local copy of the standard USNO daily IERS Bulletin A product, refreshed by an automated daily refresh procedure (Section 12.6); the maximum age between refreshes is documented in the provenance block. The data format is the standard IERS daily product, with each record carrying a Modified Julian Date (MJD) anchor, the predicted DUT1 in seconds, the polar motion (x_p, y_p) in arcseconds, and the celestial pole offsets (dX, dY) in milliarcseconds. The reader locates the bracketing pair on MJD and linearly interpolates each scalar quantity.

For epochs outside the data window, the engine extrapolates linearly from the most recent two records and stamps an explicit warning in the JSON output (extrapolation = true with the days-outside count). Linear extrapolation is reasonable for a few weeks past the data end (DUT1 changes slowly), but degrades quickly beyond a year because the leap-second timing introduces a discontinuity. The hard limit is 365 days outside the most recent data; beyond that the engine refuses to evaluate and returns an error suggesting a data refresh.

6. Frame transformations

The position vector emitted by the ephemeris module is in the GCRS (Geocentric Celestial Reference System), aligned with the BCRS at the geocenter (IERS Conventions 2010 Section 5.3). To translate into the requested frame the engine applies one of two pipelines depending on the frame_mode selection: the equinox-based chain (B, P, N) producing the true equator and equinox of date, or the CIO-based chain (Q, ERA, W) producing the International Terrestrial Reference System through the Celestial and Terrestrial Intermediate Reference Systems.

6.1 Frame bias B (ICRS to J2000)

The frame-bias matrix B aligns the ICRS axes with the J2000 dynamical mean equator and equinox. The rotation has fixed components of the order of 17 milliarcseconds (IERS Conventions 2010 Eq. 5.20):

B = R3(-dr * cos(eps0)) * R2(-dpsi * sin(eps0)) * R1(deps)
where:
  dr   = -0.0146 arcsec (right-ascension offset of the dynamical equinox)
  dpsi = -0.041775 arcsec (frame-bias longitude correction)
  deps = -0.0068192 arcsec (frame-bias obliquity correction)
  eps0 = 23.439279444444 deg (mean obliquity at J2000)

B is computed once per process and cached as a 3x3 matrix in the SOFA polyfill module. The numerical impact of skipping B is below 17 mas in any direction and is well inside the lite-engine residual; nonetheless B is always applied because the cost is one constant matrix multiplication per call.

6.2 Precession (IAU 2006 P03 / IAU 1976)

The default precession model is the IAU 2006 P03 series of Capitaine, Wallace and Chapront (2003), adopted by the IAU in 2006 in Resolution B1. It expresses the precession angles (zeta_A, z_A, theta_A; or the equivalent X, Y components in CIO-based form) as polynomials in the time argument T = (JD_TT - 2451545.0) / 36525, evaluated to fifth order:

zeta_A  =  2.650545      + 2306.083227 T + 0.2988499 T^2
           + 0.01801828 T^3 - 0.000005971 T^4 - 0.0000003173 T^5     [arcsec]
z_A     = -2.650545      + 2306.077181 T + 1.0927348 T^2
           + 0.01826837 T^3 - 0.000028596 T^4 - 0.0000002904 T^5     [arcsec]
theta_A =                  + 2004.191903 T - 0.4294934 T^2
           - 0.04182264 T^3 - 0.000007089 T^4 - 0.0000001274 T^5     [arcsec]

The precession matrix in the equinox-based path is then P(t) = R3(-z_A) * R2(theta_A) * R3(-zeta_A). The legacy IAU 1976 model of Lieske (1977) remains selectable via precession_model = iau1976 for comparison with historical literature; the residual difference between IAU 2006 and IAU 1976 is below 50 milliarcseconds at J2050 and grows to about 1 arcsec at J2200, dominated by the second-order term theta_A.

6.3 Nutation (IAU 2000A / 2000B / IAU 1980)

Three nutation models are implemented and selectable through nutation_model:

• IAU 2000A: full 1365-term luni-solar series plus 687-term planetary series of Mathews, Herring and Buffett (2002). The implementation is a port of the SOFA function iauNut00a with the canonical truncation guards. The series is parameterised in the Delaunay arguments l, l', F, D, Omega plus the eight planetary mean longitudes (Mercury through Neptune) and the general precession in longitude. Sub-microarcsecond accuracy at the J2000 epoch.
• IAU 2000B: 77-term reduced series of McCarthy and Luzum (2003) that retains 1 milliarcsecond agreement with IAU 2000A; convenient when speed matters and sub-mas pointing is not required.
• IAU 1980: legacy series of Wahr (1980), 106 terms; retained for cross-validation with pre-2003 publications. Residual against IAU 2000A reaches 18 mas in nutation longitude due to the geophysical model improvements introduced by Mathews, Herring and Buffett.

The nutation matrix in the equinox-based path is N(t) = R1(-eps_t) * R3(-dpsi) * R1(eps_0 + deps), where eps_0 is the mean obliquity at J2000 (23.4392911 deg, IAU 2006), eps_t is the mean obliquity at the requested epoch, and (dpsi, deps) are the nutation in longitude and obliquity from the chosen series. The CIO-based path bypasses the equinox altogether and uses the (X, Y, s) parametrisation introduced in Section 6.4.

6.4 NPB chain and CIO-based pipeline

Two chains are available depending on frame_mode:

# Equinox-based chain
r_TOD  =  N(t) * P(t) * B * r_GCRS
where TOD is the true equator and equinox of date,
N is the nutation matrix (Section 6.3),
P is the precession matrix (Section 6.2),
B is the frame-bias matrix (Section 6.1).

# CIO-based chain (default in cio mode)
r_CIRS = Q(X, Y, s) * r_GCRS
r_TIRS = R3(-ERA(UT1)) * r_CIRS
r_ITRS = W(x_p, y_p, s_prime) * r_TIRS
r_topo = R(phi, lambda) * r_ITRS - r_observer

The CIO-based chain is the IAU-recommended path for high-accuracy work and avoids the equinox altogether (Capitaine and Wallace 2006). The matrix Q(X, Y, s) is built from the celestial pole coordinates X and Y (computed from the IAU 2006 P03 series with IAU 2000A nutation contributions) and the small CIO locator s, which evaluates to about -94 microarcseconds at the J2000 epoch:

Q = R3(-E) * R2(-d) * R3(E) * R3(s)
sin(d) = sqrt(X^2 + Y^2)
tan(E) = Y / X
s = - X * Y / 2 + polynomial(t)
The polynomial is given in IERS Conventions 2010 Eq. 5.13 and matches
the SOFA reference function iauS06 to better than 1 microarcsec
over 1900-2100.

The polar motion matrix W uses x_p and y_p from the IERS Bulletin A snapshot and the small TIO locator s' (about -47 microarcseconds per Julian century), evaluated as in iauPom00 and iauSp00:

W = R3(-s_prime) * R2(x_p) * R1(y_p)
s_prime(t) = -47e-6 * T  arcseconds   (IERS Conv 2010 Eq. 5.13)
T = (JD_TT - 2451545.0) / 36525.

Polar motion is set to zero outside the IERS coverage window; the resulting residual contributes less than 30 m to topocentric position and less than 0.01 arcsec to apparent direction. The CIO-based chain is also the canonical path for the modular composition of the de Sitter rotation when Group 1 of the relativistic stack is active (Section 8.2).

6.5 NAIF binary Earth orientation product

When the daily-rate NAIF binary Earth orientation product is locally available, the engine reads the daily-rate Earth orientation (precession-nutation residuals, polar motion, and Delta UT1) directly from this product rather than from the IERS Bulletin A data alone. This recovers the JPL operational Earth-orientation product at the daily cadence, with an interpolation residual below 0.05 milliarcseconds in pole position. The selection is opt-in (it is not the default to keep the deployment lean) and logged in the provenance. The reader follows the standard NAIF DAF binary conventions, including automatic endianness detection.

7. Topocentric corrections

The topocentric layer transforms the geocentric apparent place produced by the frame chain into the apparent place seen by an observer at a specific WGS84 location. Four corrections are applied: the geometric subtraction of the observer's geocentric position vector, the iterative light-time correction (with optional Shapiro coupling when the relativistic stack is active), the three-component aberration (annual, diurnal, galactic), and atmospheric refraction in one of four model variants.

7.1 WGS84 geodetic position

The observer position is given in WGS84 geodetic coordinates (phi, lambda, h). It is converted to Earth-centred Earth-fixed (ECEF) Cartesian using the standard ellipsoid parameters (semi-major axis a = 6378137.0 m, inverse flattening 1/f = 298.257223563), via:

N(phi) = a / sqrt(1 - e^2 * sin^2(phi))    (prime vertical radius)
e^2 = 2 f - f^2                            (first eccentricity squared)
X = (N + h) * cos(phi) * cos(lambda)
Y = (N + h) * cos(phi) * sin(lambda)
Z = (N * (1 - e^2) + h) * sin(phi)

The geocentric latitude correction (the difference between geodetic and geocentric latitude) is propagated into the topocentric parallax. For the Moon at perigee (about 356 000 km), the parallax can reach 1.0 deg, hence topocentric and geocentric directions differ by up to one lunar diameter. The maximum diurnal parallax of the Moon, computed from the equatorial horizontal parallax pi = arcsin(R_earth / r_moon), reaches 61.5 arcmin near perigee at full Moon (Meeus 1998 Ch. 47).

When the station displacement layer of Section 3.5 is active, the WGS84 station position is corrected for solid-Earth tides (vertical and horizontal), ocean tide loading, atmospheric pressure loading, hydrological loading, deflection of the vertical (xi, eta), and pole tide. The dominant solid-Earth tide reaches 30 cm in vertical displacement at the M2 maximum (Petit and Luzum 2010 Eq. 7.5); for sub-mas pointing of fixed installations the vertical correction is essential, while for educational use it is below the residual budget and can be safely set to zero.

7.2 Light-time iteration

The position of the Moon at the observation instant is the position of the body at the slightly earlier instant t - tau, where tau is the light-travel time from body to observer. The engine solves the implicit equation

tau = | r_body(t - tau) - r_obs(t) | / c

by simple fixed-point iteration. For the Earth-Moon vector, tau is approximately 1.28 s and convergence reaches 10^-6 s in three iterations. The corresponding apparent-motion correction is approximately 0.6 arcsec per second of light time, hence approximately 0.8 arcsec for the Moon. The iteration is convergent in the strong sense (contraction mapping) because the Moon's apparent angular velocity is below 0.55 arcsec per second and the second derivative of the position with respect to time is small on the 1-second light-time scale.

When the relativistic stack is active (the advanced_gr toggle), the engine extends the light-time fixed-point with the gravitational time delay first derived by Shapiro (1964). The Newtonian iteration is augmented in its final step by the logarithmic Earth-Shapiro term:

tau = r/c + (2 G M_earth / c^3) * ln[(r_obs + r + R_AB) / (r_obs + r - R_AB)]

with r_obs the observer geocentric distance, r the Earth-Moon distance, and R_AB their separation. For ground-based observers the contribution is sub-nanosecond (a few hundred picoseconds, c * tau of the order of a few centimetres) because both endpoints lie shallow in the Earth's gravity well. The term is evaluated only once, on the final iteration of the convergence loop, so the additional cost is a single logarithm per evaluate() call. The dominant Shapiro contribution along the line of sight to the Moon comes from the Sun's potential (Section 8.1), not the Earth's, and is applied to the displayed apparent place by a separate post-hoc layer with a double-counting guard.

7.3 Stellar aberration (annual + diurnal + galactic)

Three aberration components are modelled, all referred to the same observer instant.

Annual aberration. Computed from the heliocentric velocity of the geocentre (or topocentric site) using the Smith 1980 / Klioner 2003 formulation:

u_corrected = (|u| * (u + v/c)) / |u + v/c|
where u is the unit direction vector to the source,
v is the velocity of the observer in the BCRS, taken from VSOP87 or DE440,
and c is the speed of light.

The maximum amplitude is approximately 20.5 arcsec (the so-called constant of aberration k = 20.49552 arcsec; Meeus 1998 Ch. 23). The lite engine implementation is first-order in v/c; the second-order term reaches 1.6e-7 arcsec and is included only in the SOFA polyfill path. For the Moon, which is itself moving with non-negligible velocity in the BCRS, the relevant velocity is the relative velocity of the observer (geocentric site rotated to BCRS) and the apparent place at the retarded epoch, which is what the iterative light-time loop produces.

Lite versus Full precision tier. The lite engine path uses the simplified Meeus 23.1 form Delta_lambda = -k * cos(lambda_sun - lambda) / cos(beta), which omits the e * k eccentricity correction (peak amplitude approximately 0.34 arcsec) and the Meeus 23.2 latitude term Delta_beta. The full DE440 path applies the complete Meeus 23.1 plus 23.2 with eccentricity e and longitude of perihelion pi_sun derived from Meeus 25.4 / 25.5. The JSON output exposes the boolean apparent_corrections.annual_aberration_eccentricity_term_omitted: true in Lite mode, false in Full mode. Consumers that need to distinguish precision tier without inspecting the engine selection can branch on this flag.

Diurnal aberration. The Earth-rotation contribution is approximately 0.32 arcsec * cos(phi) at the equator, decreasing with latitude. It is added at the topocentric stage and its sign is consistent with the rising/setting geometry. The constant k_d = 0.3200 * cos(phi) arcsec follows from the topocentric site velocity v = omega_earth * R_earth * cos(phi), with omega_earth the sidereal rotation rate and R_earth the equatorial radius.

Galactic aberration. Following the IAU 2000 resolution and Klioner et al. (2003), the secular acceleration of the solar system barycentre toward the Galactic centre induces a slow drift of celestial positions. The default model uses the dipole vector A_G = 5.05 microarcsec/yr aimed toward (l, b) = (264.0 deg, 48.3 deg). Both the DC component (constant offset relative to the fiducial epoch) and the AC modulation (annual variation due to Earth's heliocentric motion) are implemented. The cumulative effect over the lifetime of the calculator is below 0.05 mas; the term is enabled when frame_corrections includes galactic_aberration and is off by default.

7.4 Atmospheric refraction (4 models)

Four refraction models are selectable through refraction_model:

1. None: refraction is set to zero; useful for purely geometric reductions and for cross-validation against the JPL Horizons "geometric" output (qty=geometric).
2. Bennett 1982 (default): the marine-navigation form

   R = cot(h + 7.31 / (h + 4.4))      [arcmin, h in degrees apparent]
   R_corrected = R * (P / 1010) * (283 / (T + 273))      [P in hPa, T in deg C]

   Chromatic refraction is approximated by a Cauchy refractivity factor (1 + 0.000293 / lambda^2) relative to lambda = 590 nm. Standard atmosphere (1013.25 hPa, 10 deg C) is the implicit reference.
3. Saemundsson 1986: an alternative form well-suited to low altitudes (h < 15 deg); agrees with Bennett to within 0.01 arcmin at the horizon. The Saemundsson form takes the apparent altitude as input, rather than the true altitude, and is therefore the more natural choice when the calculator is used to plan observations from the apparent altitude on the eyepiece.
4. Mendes-Pavlis 2004: the FCULa/FCULb closed-form mapping function used for satellite laser ranging (SLR) analysis. Couples the zenith hydrostatic delay (Saastamoinen) with a wavelength- and elevation-dependent mapping derived from radiosonde profiles, with a 1-sigma residual below 0.05 mm at 10 deg elevation. This is the recommended path when the result is consumed by a Lunar Laser Ranging analysis or by a satellite-tracking pipeline. The hydrostatic and wet components are exposed separately in the JSON output.

For very low altitudes (h < 5 deg) all four models become noticeably inaccurate because the column-integrated refractivity becomes a strong function of the local atmospheric profile (temperature inversions, mirages, anomalous propagation). The calculator emits a warning in the JSON output when h < 5 deg and the refraction is applied; the warning is suppressed when refraction_model = none.

8. Relativity

The relativistic stack is opt-in via the advanced_gr toggle and the three group toggles introduced in Section 3.8. It groups five primary effects whose status ranges from confirmed by Lunar Laser Ranging (Shapiro, de Sitter, gravitational redshift) to demonstrated by Gravity Probe B (Lense-Thirring) to of educational interest in the photonic regime (Wigner rotation), plus the Group 3 sandbox parameters (Nordtvedt eta, G-dot/G drift, graviton mass, SME coefficients, cosmological constant). Each effect is exposed in its own JSON block; the underlying primary ephemeris is unchanged unless explicitly composed (the de Sitter rotation is an example, see 8.2). Some legacy diagnostics from earlier engine generations (galactic tide, Yarkovsky residual, lunar core-mantle friction) remain available as informational outputs but are not part of the apparent-place pipeline.

8.1 Shapiro time delay

The Shapiro time delay (Shapiro 1964) is the additional propagation time experienced by a light signal traversing a gravitational well, beyond the flat-space geometric path. The calculator computes both the Earth and the Sun contributions:

Delta_T = (2 G M / c^3) * ln[(r_A + r_B + R_AB) / (r_A + r_B - R_AB)]
where:
  G = 6.67430e-11 m^3/(kg s^2) (CODATA 2018)
  M = mass of the gravitating body (Earth, Moon, or Sun)
  c = 299792458 m/s (defined SI)
  r_A, r_B = distances from the gravitating body to the two endpoints
  R_AB = separation between the two endpoints

For an Earth, Moon round trip, the Earth and Moon potentials together contribute approximately 25 ns (Murphy 2013). The corresponding distance increment is c * Delta_T, of order metres at typical elongations. Near solar conjunction, the Sun term reaches hundreds of microseconds (the photons traverse the Sun's deep gravitational well) and is the dominant relativistic light-time term in the LLR error budget. The Shapiro delay was first measured in 1964 by radar reflection from Mercury and Venus (Shapiro et al. 1968); subsequent confirmations reached 0.001 percent precision via Cassini spacecraft tracking near solar conjunction (Bertotti, Iess and Tortora 2003) and continue to be a benchmark for tests of the post-Newtonian parameter gamma in the Will (2018) parametrised post-Newtonian (PPN) formalism.

The Earth term is folded into the light-time iteration of Section 7.2 by extending the Newtonian fixed point with the logarithmic correction; the Sun term is added as a residual to the displayed apparent place by the relativistic post-processing layer with a double-counting guard that subtracts the Earth term already absorbed in the loop. The implementation matches the SOFA reference function for light deflection on the canonical test vectors documented in the SOFA cookbook.

8.2 De Sitter precession (geodetic)

The Earth-Moon system traverses spacetime curved by the Sun's gravitational field. This motion induces a slow precession of the lunar orbital plane: the geodetic (de Sitter) precession, predicted by de Sitter (1916) and measured by Williams, Turyshev and Boggs (2004) at 19.2 +/- 0.13 mas/yr from 35 years of Lunar Laser Ranging data, in agreement with general relativity at the 0.7 percent level. The rate formula (Brumberg 1991) is

dOmega/dt = (3/2) * (GM_sun) / (c^2 * a_earth) * n_earth
where:
  GM_sun = 1.32712440041e20 m^3/s^2 (IAU 2015 nominal)
  c = 299792458 m/s (defined SI)
  a_earth = 1.495978707e11 m (1 AU exact, IAU 2012)
  n_earth = 1.99098659e-7 rad/s (Earth's mean motion).
Numerical magnitude: dOmega/dt = 1.92e-8 rad/yr = 19.2 mas/yr.

The cumulative rotation since J2000.0 is exposed as a diagnostic. When the relativistic stack is active, it can be composed onto the precession-nutation chain via

R_total = R_nutation * R_precession * R_bias * R_deSitter
n_hat = [0, -sin(eps_0), cos(eps_0)]   (axis perpendicular to ecliptic in J2000 equatorial frame)
eps_0 = 23.4392911 deg
theta = (3/2) * (GM_sun / (c^2 * a_earth)) * n_earth * Delta_t
R_deSitter = I + theta * [n_hat]_x   (small-angle form, for spans below 10 millennia)

For epochs within 2000 to 2200 the small-angle form is sufficient (theta < 4 arcsec); the full Rodrigues form is reserved for spans beyond ten millennia. The composition is opt-in to avoid double counting against the underlying ephemeris series: ELP-2000/82B was fitted to numerically integrated ephemerides that already include geodetic precession at leading order, and DE440 inherits the Einstein-Infeld-Hoffmann (EIH) post-Newtonian equations of motion documented by Park et al. (2021), which likewise carry de Sitter precession in their integration. With the toggle off, the baseline ephemeris is unchanged; with the toggle on, the geodetic rotation appears as a visible delta of about 0.5 arcsec at present, growing roughly linearly with elapsed time from J2000.0.

8.3 Lense-Thirring frame dragging

The rotation of the Earth produces a magnetic-like gravitational field in the linearised Einstein equations (gravito-electromagnetism). The Moon, moving through this field, experiences a Lorentz-like force whose dynamical signature is a secular contribution to the lunar nodal regression of approximately 10 microarcseconds per 30 days (Iorio 2002; Murphy, Nordtvedt and Turyshev 2007). The acceleration formula is:

a_LT = (2 / c^2) * (v x grad phi_LT)
grad phi_LT = (G / r^3) * [J - 3 (J . r_hat) r_hat]
J = I_earth * omega_earth = 5.860e33 kg m^2/s  (Iorio 2002 canonical)
B_g = (G / (c r^3)) * [3 r_hat (J . r_hat) - J]   (gravito-magnetic field)

The integrated effect on the lunar position is well below the LLR detection floor today; the term is reported as a diagnostic and is not added to the displayed apparent place. The Gravity Probe B mission (Everitt et al. 2011) confirmed the Lense-Thirring effect on a gyroscope at the 19 percent level; LAGEOS satellite laser ranging (Ciufolini and Pavlis 2004) reached 10 percent for the LAGEOS nodal precession; the lunar measurement remains an open challenge for next-generation LLR campaigns.

Source of the Earth angular momentum constant. The value J = 5.860e33 kg m^2/s used above is the figure adopted by Iorio 2002 (Class. Quantum Grav. 19, 5473) and Murphy 2013 (Reports on Progress in Physics 76, 076901) and is retained verbatim so that the reported nodal rate of approximately 0.001 mas/yr matches the canonical published prediction. More recent geodetic measurements (LAGEOS-1/2 SLR, GRACE/GRACE-FO; Chao 2005) yield a higher IAU 2015 nominal value of approximately 7.08e33 kg m^2/s based on a refined moment of inertia I_E = 8.034e37 kg m^2 and the IAU 2010 sidereal rotation rate omega_E = 7.292115e-5 rad/s. The discrepancy is approximately 17 percent and scales the predicted nodal rate proportionally to about 0.0012 mas/yr. The engine API exposes $sEarth as an override parameter on nodalRateRadPerSec() and nodalRateMasPerYear() for users wishing to substitute the LAGEOS-based value or any alternate angular momentum.

8.4 Gravitational redshift and Lunar Coordinate Time (LTC)

Photons climbing out of a gravitational well lose energy. For light emitted at the lunar surface and observed at the Earth surface, the deeper Earth potential dominates and produces a net blueshift of approximately 7e-10 in fractional frequency (Pound and Rebka 1960; Vessot et al. 1980):

Delta_nu / nu = (Phi_emit - Phi_obs) / c^2
where Phi is the Newtonian gravitational potential at each endpoint
(positive for a deeper well, sign convention of Misner Thorne Wheeler 1973).

This is the same physics that defines Lunar Coordinate Time (LTC, Ashby and Patla 2024; Kopeikin and Kaplan 2024): a clock at rest on the lunar surface ticks faster than a clock on the Earth's geoid by approximately 56 microseconds/day, with a peak-to-peak modulation of 0.22 microseconds/day driven by the lunar orbital eccentricity. The calculator reports the canonical Ashby-Patla form:

total_us_per_day(f) = 56.0199 - 0.10843 * cos(f)
where f is the lunar true anomaly

This is the IAU 2018 Resolution B2 standard for lunar-surface time-keeping (Lunar Coordinate Time); the calculator exposes both the redshift and the LTC daily rate in the JSON output, and the user can compare two observers (own location and the ISS, or own location and the lunar surface) to see how the daily LTC gain changes with altitude and latitude. The full radial-Doppler form including the Lorentz factor and the Schwarzschild metric correction at both endpoints is:

1 + z = gamma * (1 + v_rad / c) * sqrt((1 - 2 GM_obs / (r_obs c^2)) / (1 - 2 GM_emit / (r_emit c^2)))

8.5 Wigner rotation and relativistic Doppler

For light propagating through a sequence of non-collinear Lorentz boosts (a photon emitted at the Moon, received at a moving observer on Earth), the polarization vector undergoes a small rotation called the Wigner rotation. The rotation angle is computed from the chain of velocities:

theta_W = arctan((gamma_1 + gamma_2 + gamma_12 + 1)^-1 * sin(alpha) * |v1| * |v2| / c^2)
small-angle approximation:
theta_W ~ |v1 x v2| / (2 c^2)

where alpha is the angle between v1 (Earth orbital velocity) and v2 (Moon orbital velocity). For typical Earth-Moon geometries, theta_W is at the picoradian level. The same module also exposes the complete relativistic Doppler shift z_total combining radial and transverse contributions:

1 + z_total = gamma * (1 + v_radial / c)
gamma = 1 / sqrt(1 - v^2 / c^2)
v_reflector = v_CoM + omega_Moon x r_reflector_ICRS
(v_reflector formula is used when targeting a specific LLR retroreflector)

For LLR retroreflectors, the velocity at the reflector point is the lunar centre-of-mass velocity plus the cross product of the lunar angular velocity and the reflector body-fixed position rotated to ICRS via the integrated lunar PCK. The transverse contribution is below 0.1 microarcsec and the radial contribution dominates the residual; the net Doppler shift is exposed in the LLR retroreflector card alongside the geometric and Shapiro paths.

8.6 Theoretical lab (sandbox toggles)

The form exposes a small set of sandbox toggles that explore physics still under investigation. Each toggle defaults off; when enabled, the affected output line is decorated with a sandbox flag so it cannot be confused with a baseline result.

Nordtvedt eta. The Nordtvedt parameter (Nordtvedt 1968; Hofmann and Mueller 2018) tests the strong equivalence principle via Lunar Laser Ranging. In the parametrised post-Newtonian formalism:

eta = 4*beta - gamma - 3 - (10/3)*xi - alpha_1 - (2/3)*alpha_2 - (2/3)*zeta_1 - (1/3)*zeta_2
Delta_r(t) = A_eta * cos(D - D_sun)
A_eta ~ 13.1 m * eta

The current LLR upper limit is |eta| < 4.4e-4 (Hofmann and Mueller 2018), well below the calculator's lite residual; the sandbox value is reported alongside the LLR limit for context.

G-dot/G cosmological drift. Williams, Turyshev and Boggs (2012) place an LLR upper limit of about 7e-13 per year on the cosmological drift of the gravitational constant. The calculator exposes the corresponding orbital effect:

a_dot/a = -G_dot/G    (semi-major axis drift)

Yukawa-screened gravity (graviton mass). A non-zero graviton mass would screen the gravitational potential at large distances:

V(r) = -GM/r * exp(-r/lambda_g)
where lambda_g = h / (m_g * c) is the Compton wavelength of the graviton.

The current LIGO upper limit on m_g is approximately 1.27e-23 eV/c^2, corresponding to lambda_g > 1.6e16 m (Abbott et al. 2017); the lunar orbit at 3.84e8 m is far inside this scale, so the Yukawa screening is unobservable in this calculator and the value reported is illustrative.

Standard Model Extension (SME) coefficients. Kostelecky and Russell (2011) parametrise Lorentz violation in the gravity sector via a set of coefficients s_bar^ij. The relevant coupling for a lunar acceleration anomaly is:

Delta a^i_SME ~ -s_bar^ij * grad_j phi_Newton
|s_bar_eff| = sqrt((s_bar^XX)^2 + (s_bar^YY)^2 + (s_bar^ZZ)^2)

Current LLR limits on the s_bar coefficients are at the 1e-9 level (Battat, Chandler and Stubbs 2007); the calculator reports the user-supplied value alongside the limit.

Cosmological-constant term. The Schwarzschild-de Sitter metric introduces a centrifugal-like term in the Newtonian limit:

a_Lambda = (1/3) * Lambda * c^2 * r
Delta_r_Lambda ~ (1/2) * a_Lambda * t^2 ~ 4.4e-10 m ~ 0.44 nm
(integrated over the lunar lifetime, at the cosmologically observed Lambda)

This is unobservably small at the lunar scale, but the diagnostic is included for pedagogical completeness.

9. Observability

The observability synthesis layer takes the apparent place of the Moon and produces the observable quantities of interest to the user: phase and illumination, libration, eclipse imminence, occultation candidates, apsides, super and micro Moon flags, and Lunar Laser Ranging (LLR) ping distance against five operational retroreflectors. Each block is exposed in its own JSON sub-object and rendered by a dedicated view in the user interface.

9.1 Phase, illumination, libration

The phase angle phi is the Sun-Moon-observer angle and is computed from the geocentric position vectors of the Sun and the Moon evaluated at the same retarded epoch. The illuminated fraction follows directly:

cos(phi) = (-r_moon . r_sun_from_moon) / (|r_moon| * |r_sun_from_moon|)
k = (1 + cos(phi)) / 2     (illuminated fraction, 0 to 1)
age_days = (JD_now - JD_last_new_moon)

The age of the lunation (days since the last new Moon) is computed by bisection on phi = 180 deg over a 30-day window backward from the requested epoch. The principal phases (new, first quarter, full, last quarter) are listed for the next four lunations by zero-crossing of phi - n*90 deg, with n = 0, 1, 2, 3.

The optical libration in longitude (l) and latitude (b) is computed from the Eckhardt 1981 closed-form reduction in lite mode and from the integrated lunar PCK in the high-accuracy path. Selenographic colongitude, sub-Earth point, and sub-solar point are derived from the lunar principal-axis frame:

l = arctan2(-x_E, z_E) - 180 deg     (libration in longitude)
b = arctan2(y_E, sqrt(x_E^2 + z_E^2)) (libration in latitude)
where (x_E, y_E, z_E) is the Earth direction in the lunar PA frame.
Selenographic colongitude C is the lunar terminator longitude where
the Sun is rising; advances by approximately 12.19 deg/day on average.

Free libration amplitudes (Eckhardt 1981) reach 0.04 deg in longitude and 0.02 deg in latitude on the 81-day, 24.2-yr, and 75-yr periods; these are recovered to the 0.01 arcsec level when the binary PCK is loaded, otherwise they are folded into the residual.

9.2 Eclipse imminence (Besselian elements)

For the requested epoch, the engine evaluates the Besselian elements (x, y, d, f1, f2, mu) of the next solar or lunar eclipse. The Besselian framework follows Mucke and Meeus (1992) and the formulation in Meeus (1998 Ch. 54):

x, y = projection of the Moon's shadow axis on the fundamental plane
d    = declination of the shadow axis
f1   = umbra cone half-angle
f2   = penumbra cone half-angle
mu   = ephemeris hour angle of the shadow axis
Type and gamma:
  gamma = closest geocentric approach of the shadow axis (in Earth radii)
  |gamma| < 0.9972 -> central eclipse (total or annular)
  0.9972 <= |gamma| < 1.5433 -> partial eclipse
Magnitude: ratio of the obscured to total angular diameter at greatest eclipse.

For the 30 representative eclipses currently in the regression fixture set, contact times agree with the NASA Five-Millennium Canon (Espenak and Meeus 2006) to within the documented design target of one minute on principal contact times (Section 11.10). The fixture set is small and broader Five-Millennium Canon coverage is planned (Section 12.8).

The central duration (umbral or annular phase at the centre of the path) is obtained by central differencing of the Besselian shadow-axis components (x, y) over a +/- 5 minute window around the instant of greatest eclipse, projecting the resulting shadow speed onto the Earth surface and dividing the umbral diameter by the relative velocity (shadow minus rotation at the subpoint latitude). Both bracket samples (tMax minus 5 min and tMax plus 5 min) are computed with the same ephemeris engine selected for the central instant (DE440 when the SPK kernel is available, ELP2000-82B otherwise), so the derivative does not mix engines and the resulting duration is internally consistent to sub-arcsecond level on the shadow speed.

9.3 Lunar occultations

For the requested epoch and observer, the engine performs a candidate event search over three star catalogues:

• Hipparcos (ESA SP-1200, 1997): 118 218 stars to V approximately 12.4. Default catalogue for naked-eye occultations. Astrometric precision approximately 1 mas at the J1991.25 epoch.
• Tycho-2 (Hog et al. 2000): 2.5 million stars to V approximately 11.5, with proper motions (catalogue completeness above 99 percent at V = 11.5). Loaded lazily by declination band so that an occultation search consumes only the required subset of the catalogue.
• Gaia DR3 (Gaia Collaboration 2023): the deep optional path; about 1.8 billion sources, sub-milliarcsecond astrometry. The bundled subset is the Gmag < 13 cut, distributed in 18 declination bands of 10 deg each, total approximately 290 MB. Used for sub-millisecond timing of grazing events and for the deepest-occultation searches.

The implementation projects the lunar limb (corrected for libration and topocentric parallax) onto the local meridian and tests catalogue stars within a 3-degree window over a +/- 2 hour search interval centred on the requested epoch. The two-stage refinement (Meeus 1998 Ch. 34) first locates the bracketing pair via 1-minute time stepping, then bisects to 1-second precision; predicted contact times are reported with a 1-sigma timing uncertainty propagated from the EOP and ephemeris budgets. The 10 next occultations are listed in the JSON output and rendered in a dedicated UI tile.

9.4 Apsides, super/micro Moon

Perigee and apogee in the requested month are located by sign change of the radial velocity computed from the active engine. The radial velocity v_r = d|r|/dt is approximated by a centred finite difference on the geocentric distance over a +/- 4 hour window; the zero-crossing is bracketed and refined by inverse-quadratic interpolation. The calculator labels a perigee as super Moon when it falls within 90 percent of the long-term distance distribution (about 360 000 km) and the lunar phase is within 24 hours of full Moon. Symmetrically, the micro Moon tag applies near apogee with a new-Moon condition.

The 24-month horizon for super and micro Moon detection is configurable via special_months (range 1-60); the default of 24 captures the dominant 411-day perigee-syzygy beat cycle. The apparent angular diameter at perigee reaches 33.5 arcmin against 29.4 arcmin at apogee, a 14 percent variation that is visible to the naked eye and is the source of the popular super-Moon terminology.

9.5 LLR retroreflectors

Five corner-cube retroreflector arrays were placed on the lunar surface between 1969 and 1973: Apollo 11 (Mare Tranquillitatis), Apollo 14 (Fra Mauro), Apollo 15 (Hadley Rille), Lunokhod 1 (Mare Imbrium), and Lunokhod 2 (Le Monnier crater). They remain the only physical anchors for centimetre-grade Earth-Moon distance measurement (Murphy 2013). The user can select a target reflector instead of the lunar centre of mass; the topocentric laser-ping distance is then computed from the body-fixed selenographic coordinates of the array, rotated through the integrated lunar libration:

R_laser = R_Moon_CoM + M_Lib * r_reflector_PA - R_Topocentric
r_PA = R * [cos(phi_seleno) cos(lambda_seleno), cos(phi_seleno) sin(lambda_seleno), sin(phi_seleno)]
M_Lib = lunar PA-to-ICRS rotation from the binary PCK or the Eckhardt 1981 reduction.

The card reports the geometric distance, the Shapiro extra delay (Earth and Sun terms separately), and the equivalent path-length bias from the GRAIL-era core-mantle friction phase delay (Williams, Konopliv et al. 2014). The Lunokhod 1 array, lost from 1971 to 2010 and recovered by Murphy et al. (2010), is included with its updated body-fixed position; the array efficiency (range bias due to thermal degradation of the corner-cube optical surfaces) is reported when known.

The list of LLR retroreflectors physically deployed on the lunar surface (five corner-cube arrays: Apollo 11, Apollo 14, Apollo 15, Lunokhod 1, Lunokhod 2) is exposed to consumers through a dedicated accessor that excludes pseudo-targets such as the sub-Earth point and the lunar poles; those entries are geometric markers used for distance and visibility studies that need a canonical surface anchor, not retroreflector hardware, and a separate accessor exposes every catalog row for callers that need the surface markers as well.

10. Outputs

The engine emits a single typed JSON payload with versioned schema (schema_version field; see Section 12.3). The web interface renders this payload through dedicated views (Bento cards, scientific telemetry bar, libration compass, illumination ring, eclipse timeline, occultations, retroreflectors, sandbox panels). The payload is also available verbatim through the JSON download button and the Copy JSON action; for academic citation, the BibTeX block is embedded in the same payload at the citation path.

10.1 Coordinates with precision

Right ascension and declination are reported in the requested frame (ICRS, J2000.0, mean of date, true of date) with sub-arcsecond precision in lite mode and sub-milliarcsecond precision when DE440 + IAU 2000A are active. Ecliptic coordinates (lambda, beta) follow the same precision band as RA/Dec. Each output value carries a 1-sigma uncertainty propagated from Section 11; the central value and the 1-sigma figure share the same format string so that the displayed precision matches the actual numerical content.

{
  "ra_deg": 152.123456,        "ra_sigma_arcsec": 0.05,
  "dec_deg": -8.234567,        "dec_sigma_arcsec": 0.05,
  "lambda_deg": 154.789012,    "lambda_sigma_arcsec": 0.05,
  "beta_deg": -3.456789,       "beta_sigma_arcsec": 0.05
}

10.2 Topocentric Az/Alt and rise/set/transit

Topocentric altitude and azimuth are computed at the requested instant via the standard formula:

sin(h) = sin(phi) sin(dec) + cos(phi) cos(dec) cos(H)
sin(A) = -sin(H) cos(dec) / cos(h)
H = LST - RA   (local hour angle)

Rise, transit, and set times are computed by bisection of altitude = h_apparent (default -0 deg 34 arcmin 0 arcsec, the standard lunar refraction-corrected horizon), with a residual below 1 second over the 1900-2100 window. The bisection bracket is refined by quadratic interpolation when the lunar declination crosses the latitude limit, which avoids the spurious solutions that appear at high latitudes near the polar night/day boundaries.

10.3 ICRS state vectors

The Cartesian state vector (x, y, z, vx, vy, vz) is exposed in the ICRS at the requested epoch in km and km/s. This is the primary input expected by the CCSDS 502.0-B-3 OEM exporter (Section 11.8). Velocities are computed from the analytical theory by differentiating the series with respect to time when ELP/VSOP is the active engine, and from the Chebyshev derivative when DE440 is active.

10.4 Eclipse parameters

Gamma, magnitude, type (partial / annular / total / hybrid), contact times (P1, U1, U2, U3, U4 where applicable), and a path summary for total/annular events. The path summary is a coarse 10-point sampling of the umbra ground track, sufficient for a first-look planning map but not for the centimetre-grade contact times required for academic eclipse reports.

10.5 Monthly ephemeris table (30 days)

A 30-day strip surrounding the requested epoch with daily rows: UTC date, distance, RA/Dec, illumination, phase angle, libration. Useful for programming a month of observations from a single page; the format is parallel to the daily ephemerides published by the Astronomical Almanac and the Connaissance des Temps.

10.6 Time scales and provenance

The five time scales (UTC, UT1, TAI, TT, TDB), plus optional TCG, TCB, TCL, GPS, are reported in both Julian Date and ISO 8601 form. The provenance sub-object records the engine identifier, the engine version, the active model selection, the IERS snapshot identifier and age, the active sandbox overrides, and the deterministic Reproducibility ID introduced in Section 12.2.

10.7 Optional outputs (toggleable cards)

The form exposes a set of independent display-quantity toggles that control which result cards are rendered. The numerical engine always computes the full payload; the toggles only control rendering. This means a permalink with all toggles off still produces a complete JSON payload, and a third party can compute against the canonical contract without depending on the user's display preferences. The toggles are exposed in the URL for reproducibility (a permalink that opens the page with a specific set of cards visible is itself a citable artefact).

10.8 Result UI: Panel, Telemetry, Bento, Apparent corrections

The user-facing rendering layer consumes the JSON payload through dedicated presentational views. Each view is a thin layer over the engine output; none of them runs additional astronomical computation. The list below pairs each visible UI element with the JSON path it reads (the canonical contract column), the unit it renders in, and the methodology subsection where the underlying physics is documented.

10.8.1 Panel (HERO block, standalone section below Sandbox)

The Spotlight is a single-glance summary suitable for citation in a slide deck or a cross-engine comparison. After the 2026-05-09 layout reorder, it lives as a standalone section directly below the Sandbox (Constant Overrides) block, no longer at the very top of the result panel. The block reads from the lunar position, phase, topocentric and libration sub-objects of the JSON payload and renders the following quantities:

Limitations. The Spotlight does not currently expose 1-sigma uncertainties next to each quantity; the per-quantity expanded uncertainties live in the uncertainty sub-object of the JSON payload, and a future iteration will add inline error bars. The 3-decimal distance format is a display convenience: the true precision of the displayed value is governed by the engine residual (Section 11 and 12.8), which is at the kilometre level in default mode and is targeted at the sub-kilometre level in the high-precision DE440 path subject to ongoing validation.

10.8.2 Telemetry bar (live, 4 cells)

A horizontal bar with four live-updated cells driven by a small front-end animation loop. The bar is purely cosmetic for the time scales (UTC, TT, TDB, JD); the engine's authoritative time scale is always the one stamped into the JSON payload. The four default cells are:

• Cell 1: UTC clock. Reads the UTC time scale from the payload and increments at requestAnimationFrame; resyncs on every result refresh. Format: ISO 8601.
• Cell 2: TT clock. Reads the Terrestrial Time scale from the payload and applies the TT = TAI + 32.184 s + Delta-AT identity. Format: ISO 8601.
• Cell 3: Lunar Coordinate Time. Reads the lunar coordinate time field from the relativistic stack of Section 8.4 when the Group 2 toggles are active; shows the canonical Ashby-Patla daily rate alongside the cumulative offset. Pauses when Group 2 is off.
• Cell 4: Snapshot age badge. Reads the IERS snapshot age field from the provenance block; goes amber above 30 days and red above 90.

Pause behaviour. The bar pauses on the page-visibility event (when the tab is backgrounded) and on the before-unload event to release the animation callback; it resumes when the tab returns to the foreground.

10.8.3 Apparent corrections card

A dedicated card that itemises every correction applied between the geometric ICRS direction and the apparent topocentric direction reported in the Spotlight. The values come from the apparent-corrections sub-object of the JSON payload and are shown line by line so that the user can audit the path. The four lines are:

Interpretation. The card is intended as a sanity check before a result is cited or compared to an external tool: if the user expected geometric coordinates and the card shows a non-zero refraction line, the form has the refraction toggle in the wrong state. Each line is rendered grey when the corresponding model is set to off.

10.8.4 Provenance block (orientation-product status, data hashes, ephemeris versions)

The provenance block is the audit trail of a single computation. It is rendered as a collapsible card under the result panel and is the canonical citation target for academic work. The block exposes the following fields:

• Engine identifier and version. A short identifier denoting the default-mode or high-precision branch, plus a semantic-version string for the engine release that produced the result.
• Active model selection. The full set of model selectors active for this call (nutation, precession, frame origin, refraction, aberration, advanced relativity, plus the Group 1, 2, 3 toggles). Each field is the literal value emitted by the form.
• IERS data identifier. The IERS snapshot issue date, the integer Delta-AT value at the requested epoch, the UT1 minus UTC value in seconds, the snapshot age in days, and a Boolean flag indicating whether the result was extrapolated past the snapshot end.
• Orientation-product status. A compound status field indicating whether the lunar physical-orientation product and the Earth high-precision orientation product were available for this call, plus content-hash fingerprints of each product when present.
• Sandbox overrides. An array of constants overridden in this call (Section 3.7), with the value used. The array is empty when no override is active.
• Reproducibility ID. The deterministic, human-readable identifier described in Section 12.2.

10.8.5 Bento tile (G heatmap) and standalone sections (J occultations, K uncertainty budget)

Layout reorder of 2026-05-09: the result panel now opens with a single Bento tile (G, monthly libration heatmap) directly below the Telemetry bar, followed by the Sandbox section (Constant Overrides), then three standalone sections containing what were previously the Spotlight HERO and Bento tiles J and K. The Spotlight, occultations and uncertainty budget were promoted out of the dashboard grid because they answer different questions: the Spotlight is the single-glance summary, the occultations table is a forecasting tool over the next 30 days, and the uncertainty budget is a traceable accounting of the result. Folding all three into a single grid above Sandbox flattened the visual hierarchy and made it harder to scan the result for the most critical numbers first. The older tiles (E illumination ring, F next event) were removed in 86.7.3.0.0 to reduce redundancy with the Spotlight; the stochastic gravitational-wave background jitter waveform tile was removed in 86.8.x because the displayed values (RMS amplitude, characteristic strain h_c) were literature constants from NANOGrav 15-yr (Agazie et al. 2023), not a live measurement, and the visualisation risked being read as real-time detection. Sandbox, libration compass and eclipse timeline tiles are not currently rendered in the Bento and are documented in their own subsections (3.7, 3.8 and 9.x) when reactivated.

• Tile G (Bento, above Sandbox): Monthly libration heatmap. A 30-day by 24-hour grid of libration in longitude, libration in latitude, and topocentric distance rendered with a third-party charting library; falls back to a plain HTML table when the chart library fails to load. Source: the libration-grid event series of Section 9.1.
• Section J (standalone, below Sandbox): Next 10 lunar occultations. Reads the occultation event list and renders the predicted contact times (D = disappearance, R = reappearance) with star identifier, apparent visual magnitude, and 1-sigma timing uncertainty propagated from EOP and ephemeris. Catalogue selectable between Hipparcos, Tycho-2 and Gaia DR3. Source: Section 9.3.
• Section K (standalone, below Sandbox): Uncertainty Budget. The Section 11 budget rendered as a stacked bar (per-source contribution to the root-sum-square total) plus a summary row with the combined standard uncertainty u_c, the expanded uncertainty U at the requested coverage factor, the effective degrees of freedom from Welch-Satterthwaite, and the Lense-Thirring nodal-rate diagnostic when Group 2 is active. Optional Monte Carlo cross-validation panel (10000, 100000 or 1000000 trials) and Allan deviation table (with fixed-layout column widths from 2026-05-09 to keep the tau and sigma_y columns aligned across rows of varied exponent length); CCSDS 502.0-B-3 OEM v3.0 export available as a data URL.

10.8.6 Reproducibility ID (visible in result footer)

The Reproducibility ID is rendered verbatim as a copy-button-flanked text field in the result footer. The ID is the canonical citation handle: a publication that quotes the ID and the engine version can be reproduced bit-for-bit at any future date. See Section 12.2 for the construction rule.

10.8.7 CCSDS OEM export

A button labelled Download CCSDS OEM in the result footer triggers the CCSDS 502.0-B-3 v3.0 export (Section 11.8). The output is a UTF-8 text block with the canonical headers, a single-epoch META block, the 6-vector state line, and the 21-element lower-triangular covariance. The Reproducibility ID and the engine version are embedded as COMMENT lines so the source pipeline is recoverable from the file alone.

10.8.8 JSON output schema

The full JSON payload is exposed verbatim through a Download JSON button and a Copy JSON action in the result footer. The schema carries an explicit version field (currently version 3); semver discipline is described in Section 12.3. Backwards-incompatible changes emit a deprecation banner for at least one release cycle.

11. Uncertainty quantification (GUM/JCGM/CCSDS)

Uncertainty quantification follows the Joint Committee for Guides in Metrology framework. The calculator implements the canonical GUM (JCGM 100:2008), the Monte Carlo supplement (JCGM 101:2008), the multivariate extension (JCGM 102:2011), and the CCSDS 502.0-B-3 export. Allan deviation diagnostics follow NIST SP 1065 (Riley 2008). For a topocentric apparent direction the angular uncertainty is dominated by ephemeris truncation in the lite mode, and by the IERS UT1 snapshot age in the DE440 + SOFA mode; the per-source breakdown follows.

For the geocentric distance the budget is dominated by the ELP-2000/82B truncation residual (kilometre level) plus a small DE440 numerical-coefficient quantisation when the high-precision path is active (sub-decametre level). For phase, the relevant quantity is the Moon-Sun elongation, whose uncertainty is at most about 6 arcsec in the default mode.

11.1 Type A vs Type B inputs

Each input is classified into Type A (statistical, derived from a sample of repeated observations or from a published statistical analysis) or Type B (evaluated from a-priori information: instrument specifications, published standard uncertainties, expert judgement). Examples: Bulletin A DUT1 is Type A from the IERS VLBI/GNSS combination (about 1000 daily samples per quarter, with a documented combined formal uncertainty); the IAU 2000A nutation truncation residual is Type B from the published series amplitude with a rectangular distribution over the residual interval. The Type A and Type B classification governs the choice of distribution in the Monte Carlo path and the choice of degrees of freedom in the Welch-Satterthwaite formula.

11.2 Sensitivity coefficients (Jacobian)

For each output quantity y_k = f_k(x_1, ..., x_n) the partial derivatives c_ki = df_k/dx_i are computed numerically by central difference (with adaptive step), assembled into a Jacobian J = [c_ki]:

c_ki = (f_k(x_1, ..., x_i + h_i, ..., x_n) - f_k(x_1, ..., x_i - h_i, ..., x_n)) / (2 h_i)
h_i = max(|x_i| * eps_machine^(1/3), h_min_i)
where eps_machine = 2.22e-16 for IEEE 754 double precision
and h_min_i is a per-input floor protecting against degenerate cases (x_i = 0).

The numeric step is chosen by the standard finite-difference rule (h ~ x * eps_machine^(1/3) balances truncation against round-off), with safeguards against degenerate inputs (x_i = 0 forces the floor). The Jacobian is computed once per output quantity and cached; the cost is 2*N_inputs additional engine evaluations, dominated by the IAU 2000A nutation series and amounting to about 36 ms additional latency for a typical request.

11.3 Combined u_c via RSS with correlation matrix

For each output y_k, the combined standard uncertainty u_c(y_k) is computed in the GUM form including correlations:

u_c^2(y_k) = sum_i  c_ki^2 u^2(x_i)
           + 2 sum_{i<j}  c_ki c_kj u(x_i) u(x_j) r(x_i, x_j)
where:
  c_ki = sensitivity coefficient (Section 11.2)
  u(x_i) = standard uncertainty of input x_i
  r(x_i, x_j) = correlation coefficient (in [-1, +1])

The full input correlation matrix is supplied: DUT1 and Delta-AT are uncorrelated (different physical processes); x_p and y_p are correlated through the IERS combined product (typical r about 0.4 from the LS combination); refraction inputs (P, T, RH) are correlated via meteorological covariance (typical r about 0.6-0.8 between P and T at synoptic scales, weaker between RH and the others).

11.4 Welch-Satterthwaite degrees of freedom

The effective degrees of freedom of u_c(y_k) follow the Welch-Satterthwaite formula:

nu_eff = u_c^4(y_k) / sum_i ((c_ki u(x_i))^4 / nu_i)

nu_i is taken from the published Type-A degrees of freedom (Bulletin A combination provides nu_DUT1 of order 200 from the daily samples per analysis cycle) or set to infinity for Type-B inputs whose standard uncertainty is given without further information. The Welch-Satterthwaite formula is approximate but is the canonical GUM choice; for cases where it breaks down (small nu_eff with strongly non-Gaussian inputs), the Monte Carlo path of Section 11.6 is the recommended fallback.

11.5 Coverage factor and expanded uncertainty

The expanded uncertainty U is reported at three coverage levels:

• k = 1 (approximate 68 percent coverage when nu_eff is large);
• k = 2 (approximate 95 percent coverage); default reported next to the apparent place;
• k = 3 (approximate 99.7 percent coverage); used when a wider safety margin is preferred.

For nu_eff small, the coverage factor is replaced by the Student's t value t_p(nu_eff) per JCGM 100 G.2:

k = t_{0.95}(nu_eff)        (95 percent two-sided coverage)
For nu_eff = 5,    k = 2.571
For nu_eff = 10,   k = 2.228
For nu_eff = 30,   k = 2.042
For nu_eff -> inf, k -> 1.960 (Gaussian limit)

The internal lookup table holds anchor rows at nu = 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100 (NIST SP 1065 / JCGM 100 Tab G.2 values) and linearly interpolates between bracketing keys. Adding the intermediate rows (4, 25, 40, 75, 100) keeps the maximum interpolation error below 0.005 across nu in [3, infinity), which is the regime the Welch-Satterthwaite formula tends to populate. Above nu = 1e6 the asymptote k -> 2.0 is returned by engineering convention (CCSDS); consumers that require the strict normal limit k = 1.960 should consult Section 12.7.

11.6 Monte Carlo (JCGM 101:2008)

For non-linear cases or when the linearised Jacobian is unrepresentative, the engine runs an optional Monte Carlo simulation with M = 10^6 trials. Inputs are drawn from their Type-A or Type-B distributions (Gaussian, Student's t, rectangular, triangular, U-shape as appropriate), propagated through the full evaluate() call, and the output empirical distribution is summarised by mean, standard deviation, 95 percent coverage interval (shortest), and 95 percent probabilistic-symmetric interval. The result is consistent with the GUM linear approximation in the regime where the latter is valid (relative deviation below 1 percent); when the deviation exceeds 1 percent, the calculator emits a warning and the Monte Carlo result is the recommended one.

The Monte Carlo path uses Box-Muller for Gaussian sampling, inverse-CDF for triangular and U-shape, and standard rejection for the Student's t when nu_eff is small. The sample size M = 10^6 is the JCGM 101 recommendation for sub-percent precision on the 95 percent interval; smaller trial counts are available through the form's Monte Carlo selector for educational use.

11.7 Output covariance (JCGM 102:2011)

The full output covariance matrix Sigma_y = J Sigma_x J^T is computed and exposed in the JSON payload. Off-diagonal terms (e.g., the correlation between RA and Dec, between altitude and azimuth, between the geocentric distance and the topocentric parallax) are essential for downstream consumers that combine multiple outputs (occultation timing, mosaic planning, multi-station coincidence). The covariance is exported in the OEM block (Section 11.8) when the user requests CCSDS export.

11.8 CCSDS 502.0-B-3 OEM export

The output state vector (Section 10.3) and its full covariance can be exported in the Orbit Ephemeris Message format defined by CCSDS 502.0-B-3 (CCSDS 2023). This is the standard interchange format used by NASA, ESA, JAXA, and major commercial operators for spacecraft and natural-body ephemerides. The OEM is emitted as a UTF-8 text block with the canonical headers:

CCSDS_OEM_VERS = 3.0
CREATION_DATE = 2026-05-09T15:30:00
ORIGINATOR = OcalEN

META_START
OBJECT_NAME = Moon
OBJECT_ID = 301 (NAIF body ID)
CENTER_NAME = EARTH
REF_FRAME = ICRF
TIME_SYSTEM = TDB
START_TIME = ... STOP_TIME = ...
META_STOP

(state vector lines: epoch, x, y, z, vx, vy, vz)

COVARIANCE_START
EPOCH = ... REF_FRAME = ICRF
(21 elements of the lower-triangular 6x6 covariance, in km, km/s units)
COVARIANCE_STOP

The OEM is downloadable from the result panel and is suitable for direct ingest into NASA's ODM software, ESA's ESOC tools, or any tool implementing the CCSDS standard. The reproducibility ID and the engine version are embedded in the OEM as COMMENT lines, so the source pipeline is recoverable from the file alone.

11.9 Allan deviation (NIST SP 1065)

For users who consume a time series of OcalEN outputs (e.g., a daily LTC drift measurement, or a multi-decade occultation residual stream), the calculator can compute the overlapping Allan deviation sigma_y(tau) following NIST SP 1065 (Riley 2008):

sigma_y^2(tau) = 1 / (2 (M - 2 m + 1) tau^2) sum_{i=0..M-2m} (y_{i+m} - y_i)^2
where m = tau / tau_0 (averaging factor)
and y_i are the consecutive frequency or phase samples.

The result characterises the noise type (white phase, white frequency, flicker, random walk) by the slope of the log-log Allan plot and is rendered in the diagnostics panel. The averaging windows are configurable; the default is a logarithmically spaced grid from tau_0 (the sample interval) to the half-length of the input series.

11.10 Numerical tolerance thresholds (public targets)

The table below declares the tolerance targets that the engine aims for in the canonical DE440 + IAU 2006 + IAU 2000A configuration. They are targets, not certified values: validation against these thresholds is in progress and the current sample set is small. The validation status, the broader baseline plan and the current measured residual are reported in Section 12.8.

These targets are the same ones consumed by the validation status of Section 12.8: a fixture that exceeds the target tolerance is reported as a regression and surfaced in the build log. The targets are conservative for the high-precision configuration; the default-mode residual is dominated by the analytical-truncation contributions of Section 11 Table 11.1 and is not meant to meet these thresholds.

12. Reproducibility and validation

This section documents the four primitives by which an OcalEN computation is reproducible: the canonical permalink, the deterministic Reproducibility ID, the versioned JSON schema, and the regression baseline against JPL Horizons. The limitations subsection (sec-limitations) is embedded in 12.7 as the canonical declaration of what the calculator does not do.

12.1 Permalink format

Every result panel exposes a permalink button. The permalink is a canonical query string with parameters in alphabetical order, with the SHA-256 hash of the canonicalised payload appended as the trailing parameter h=. Two clients that supply the same set of inputs will receive identical permalinks; a third party can verify that a permalink has not been modified by recomputing the hash.

https://www.ocalendario.com.br/calculadora-lunar-cientifica
?date=2026-05-09&time=22:30:00&lat=-23.55&lon=-46.63&alt=760
&engine=de440&nutation_model=iau2000a&refraction_model=bennett
&h=4f7c...e2a1

Live URL state synchronisation is implemented: every change to any form parameter triggers a 300 ms-debounced update to the browser address bar via the history.replaceState() API. The URL contains the canonicalised parameter set (sorted alphabetically by key) plus a SHA-256-derived 16-character hash for citation integrity. The hash is truncated to 16 hexadecimal characters (64 bits), which provides sufficient collision resistance for citation purposes (search space of order 1.8e19); the full 64-character SHA-256 is available via the JavaScript Web Crypto API (crypto.subtle.digest), which requires HTTPS or localhost. Browsers without history.replaceState or without the Web Crypto API fall back to plain query-string permalinks.

12.2 Reproducibility ID

In parallel with the canonical permalink, the engine generates a deterministic, human-readable Reproducibility ID built from a slug derivation of the inputs:

id_parts = [
  date, time, "tz" + timezone,
  "lat" + lat, "lon" + lon,
  engine + "v" + engine_version,
  frame_mode, nutation_model, precession_model
]
reproducibility_id = "lunar-" + id_parts.join("-")
example: lunar-2026-05-09T22-30Z-saopaulo-de440v86.9.29.0.74-cio-iau2000a-iau2006

The ID is stable across platforms and embeds the engine version, so a citation that quotes the ID can be reproduced bit-for-bit at any future date. If two requests return the same ID, the underlying numerical computation is byte-for-byte identical (modulo non-determinism in concurrent IERS snapshots, ruled out by the same-version constraint). The ID is human-readable and can be inspected directly without decoding.

12.3 JSON schema versioned

The JSON payload carries a schema_version field (currently 3). Schema changes follow semver: additive changes (new fields) bump the minor; backward-incompatible changes (renamed or removed fields) bump the major and emit a deprecation banner for at least one release cycle. The schema shape is summarised below:

{
  "$schema": "https://www.ocalendario.com.br/schema/ocale-output.v3.json",
  "version": "86.9.29.0.74",
  "inputs": { ... canonical-sorted GET params ... },
  "provenance": {
      "engine": "default-mode",
      "engine_version": "...",
      "models": { "nutation": "iau2000a", "precession": "iau2006", ... },
      "iers_snapshot": { "issue_date": "...", "delta_at": 37, "dut1_s": ..., "age_days": ... },
      "reproducibility_id": "lunar-..."
  },
  "time_scales": { "utc": ..., "ut1": ..., "tai": ..., "tt": ..., "tdb": ... },
  "moon": { "ra_deg": ..., "dec_deg": ..., "distance_km": ..., ... },
  "topocentric": { "az_deg": ..., "alt_deg": ..., "rise": "...", "set": "...", ... },
  "phase": { "name": "...", "illumination_pct": ..., "age_days": ..., ... },
  "uncertainty": { "angular_mas_1sigma": ..., "components": { ... } },
  "events": { "next_phases": [...], "next_eclipses": [...], "occultations": [...] },
  "advanced_gr": { "shapiro_ns": ..., "de_sitter_mas_per_yr": ..., ... },
  "self_url": "https://...permalink..."
}

12.4 Horizons regression baseline (status)

The current Horizons cross-check is a smoke fixture set (a small number of representative epochs, with wide tolerances) used to detect gross regressions between releases. It is not yet the broader JPL Horizons baseline that the methodology aims for: the planned baseline covers 200 to 1000 epochs distributed over 1900-2100 across multiple latitudes (including polar cases), stratified by lunar phase, and matched against the engine in the high-precision DE440 plus full IAU 2000A nutation configuration. The smoke fixtures, the broader baseline plan, the public tolerance thresholds, and the current measured residual are documented together in Section 12.8.

12.5 Internal tests

Three test layers run on every release. The status of each layer (pass count, sample size, what it proves) is summarised in the validation table of Section 12.8.

1. A unit-test layer covering the analytical theories and frame transformations against the Meeus 1998 worked examples, plus parameter-influence audits and toggle-matrix audits that exercise every documented input and every engine selector;
2. An algorithm-consistency layer covering engine consistency, bounds invariants, monotonic propagation, identity rotation when EOP is zero, and Pauli-like algebraic identities for the rotation chain;
3. A smoke layer covering a small JPL Horizons cross-check (limited samples; broader baseline pending), the Five-Millennium Canon eclipse contact times across the fixture set defined in Section 12.4, and the LLR retroreflector path against published Apollo 11 round-trip values (Bender et al. 1973).

All tests are pure-PHP and run without external dependencies. The test suite executes in approximately 90 seconds on a typical CI runner, dominated by the Monte Carlo regression (a reduced trial count in the test path versus the production path) and the Horizons fixture comparison.

12.6 IERS data refresh

Three automated refresh procedures maintain the operational data:

• A daily refresh of the IERS Bulletin A data product from the USNO public mirror;
• A weekly refresh of the IERS C04 long-term combined product;
• A weekly refresh of the NAIF binary Earth orientation product.

The data is timestamped; the provenance block exposes the data age and a warning is emitted when it exceeds the configured threshold (default 30 days). A maintenance-key-protected administrative endpoint allows manual on-demand refresh of any data product; this is the recommended path for institutional users who want the latest IERS values without waiting for the next automated tick.

12.7 Scope limits and known issues

This subsection is the explicit declaration of what OcalEN does not do and of the documented approximations behind what it does do. It is the contract with the user: any behaviour outside what is enumerated here is either a bug (please report) or a future feature (not yet implemented). The list is updated whenever a new feature lands or a new limit is documented.

12.7.1 Out of scope (the calculator does not do these)

• Real-time orbit determination. The HTTP request and response lifecycle adds about 50 ms of latency, the engine is stateless, and there is no real-time orbit determination layer.
• Artificial satellites. The engine has no TLE or SP3 ingest, no orbit propagator for spacecraft, and emits no NORAD or COSPAR identifiers.
• Real-time IERS tracking. DUT1, polar motion (x_p, y_p), and Delta-AT come from a daily-cron-refreshed snapshot of IERS Bulletin A, not from a live request. The maximum age between refreshes is 24 hours by default; the engine refuses to evaluate when the snapshot is older than 365 days.
• Sub-millimetre science. Lunar Laser Ranging at the millimetre range level, optical-clock comparisons at the 10^-18 level, and relativistic geodesy require femtosecond-grade timing and station displacement modelling beyond the precision targeted here.
• Safety-of-life or certified navigation. The output carries no SBAS/GBAS integrity flags and is not certified for aviation, marine, or any safety-of-life application.
• Partial date/time inputs. A full date plus time (or a full JD) is required; ranges and astropy-style approximate inputs are not accepted.

12.7.2 Model and ephemeris approximations (documented residuals)

• No DE441. Only DE440 is supported (1550-2650). DE441's extended range (-13200 to +17191) is not implemented; outside the DE440 window the engine falls back to the analytical default path and the additional residual is folded into the uncertainty budget of Section 11.
• DE440 numerical-coefficient quantisation. The bundled DE440 path uses publicly delivered Chebyshev coefficients; the corresponding apparent-direction residual at the Moon is sub-arcsecond at J2000 and grows toward the kernel boundaries (Park et al. 2021 Table 5).
• Truncated ELP-2000/82B in default mode. The default Moon path uses a compact subset of the ELP-2000/82B series. The truncation residual is approximately 3 to 5 arcsec in apparent direction over the 1900-2100 window; outside that window the residual grows because of Delta-T uncertainty in the mean motion arguments. The high-precision DE440 path is the recommended option when sub-arcsecond accuracy is required.
• Truncated VSOP87D in default mode. The default Sun path retains a compact term subset; the residual at J2000 is sub-arcsecond in heliocentric longitude and grows to several arcseconds beyond a 4000-year span from J2000.
• IAU 2000A nutation in the default path. The default path uses a compact subset of the IAU 2000A series, retaining the dominant signal power. The high-precision mode evaluates the full IAU 2000A series; both paths are documented in the engine's provenance field.
• No asteroid perturbations exposed. The four largest asteroids (Ceres, Pallas, Vesta, Hygiea) contribute about 10 microarcsec to the Moon's apparent position. When DE440 is loaded, their effect is absorbed into the integrated lunar ephemeris but is not exposed as a separate output line.
• Outer planets in lite mode. Jupiter, Saturn, Uranus, Neptune are not directly modelled in the lite path; their gravitational perturbations on the Moon are absorbed into the ELP-2000/82B fitted constants. When DE440 is loaded they are available for the relativistic perturbation analysis of Section 8.
• Galactic aberration off by default. The IAU 2006 Resolution B1 recommendation to include galactic aberration is implemented but disabled by default; cumulative effect over the calculator lifetime is below 0.05 mas. The user can enable it from the modern reference systems panel of Section 3.6.
• Free Core Nutation default amplitude. The Lambert 2007 empirical FCN model is implemented with a default amplitude of 0.05 mas; the residual against the IERS solution can reach 0.1 mas in years where the FCN amplitude departs from this default.

12.7.3 Atmosphere and refraction

• Mendes-Pavlis 2004 standard-atmosphere coupling. The Mendes-Pavlis path uses pressure, temperature and humidity at the station as scalar inputs; it does not couple to a radiosonde profile or to a synoptic NWP forecast. Day-to-day boundary-layer variations (subsidence inversions, sea-breeze fronts) are not represented; the residual at 10 deg elevation can exceed the nominal 0.05 mm zenith delay precision when the local boundary layer is anomalous.
• Bennett 1982 below 5 deg. All four refraction models become noticeably inaccurate below about 5 deg of altitude because the column-integrated refractivity becomes a strong function of the local profile (temperature inversions, mirages, anomalous propagation). The calculator emits a warning when h < 5 deg and refraction is applied.
• Chromatic refraction approximation. The Cauchy form (1 + 0.000293 / lambda^2) used for chromatic refraction is accurate at the 0.05 arcsec level over 400 to 800 nm; outside this band (UV or IR) the residual grows and the constant from the Edlen 1966 form should be used instead. Not yet implemented.
• Aerosol optical depth not coupled. The form accepts AOD_550 as a parameter but the value does not currently feed any model; it is reserved for a future extinction-correction layer in the magnitude calculation.

12.7.4 EOP and time-scale snapshots

• Linear extrapolation past snapshot. When the requested epoch falls outside the IERS Bulletin A snapshot, the engine extrapolates DUT1 and polar motion linearly from the most recent two records and stamps an explicit extrapolation = true flag in the JSON output. Linear extrapolation is reasonable for a few weeks past the snapshot end; beyond several months it degrades because the leap-second timing introduces a discontinuity. Hard limit: 365 days outside the snapshot.
• UTC convention before 1972. For epochs before 1972-01-01, UTC is not formally defined; the engine uses the convention DUT1 = 0 and Delta-AT = 10, which matches USNO contemporary practice and the Espenak-Meeus 2006 historical Delta-T extrapolation.
• Espenak-Meeus extrapolation residual. The Espenak-Meeus 2006 polynomial extrapolation for Delta-T outside 1972 to present grows quickly with epoch. At the Late Antique boundary (year 500 CE) the residual reaches approximately 30 arcsec on apparent direction; at the Babylonian boundary (year -500) it reaches several arcminutes. Stephenson, Morrison and Hohenkerk 2016 long-term tabulation would refine this but is not bundled.

12.7.5 Libration, occultations, eclipses

• Eckhardt 1981 free-libration amplitudes when the binary PCK is absent. Without the NAIF lunar binary PCK kernel, the free-libration amplitudes (0.04 deg in longitude, 0.02 deg in latitude over 81-day, 24.2-yr and 75-yr periods) are recovered from the Eckhardt 1981 closed form, with a residual of order 0.1 arcsec on the libration axes. With the binary PCK loaded the residual drops to 0.01 arcsec.
• Occultation catalogue scope. Hipparcos, Tycho-2 and Gaia DR3 are the only stellar catalogues supported (Section 9.3); occultations against extragalactic sources, deep-field surveys, or proper-motion-binary companions are not predicted. Targets outside these catalogues require an external star coordinate as input.
• Eclipse path-summary granularity. The path summary for total/annular eclipses is a coarse 10-point sampling of the umbra ground track; sub-second contact times for academic eclipse reports require Five-Millennium Canon level integration, which is implemented in the regression test but not in the user-facing summary.
• LLR retroreflector array efficiency. The reported laser-ping path uses the geometric distance plus Shapiro plus integrated lunar libration; the array efficiency degradation (Murphy et al. 2010 for Lunokhod 1) is reported when known but is not propagated into the timing budget.

12.7.6 Uncertainty quantification

• Welch-Satterthwaite approximation. The combined standard uncertainty u_c uses the Welch-Satterthwaite formula for the effective degrees of freedom; for strongly non-Gaussian inputs with small nu_eff the approximation degrades. The Monte Carlo path of Section 11.6 is the recommended fallback in those cases, and the engine emits a warning when the GUM-linear and Monte Carlo results disagree by more than 1 percent.
• Monte Carlo sample size at default. The default Monte Carlo trial count is M = 10⁶, which provides sub-percent precision on the 95 percent coverage interval. Educational use with smaller trial counts (10⁴ in the test path) is supported through the form's Monte Carlo selector but the resulting 95 percent interval has a wider sampling error.
• Correlation matrix coverage. The engine ships with a documented correlation matrix for the EOP inputs (DUT1 and polar motion uncorrelated; x_p and y_p with r approximately 0.4 from the IERS LS combination) and the meteorological inputs (P, T, RH with r approximately 0.6 to 0.8). Correlation between epoch and EOP is not currently modelled.

12.7.7 Data freshness and operational

• Data age threshold. A warning is emitted when the IERS Bulletin A data age exceeds the configured threshold (default 30 days, hard limit 365 days). For sub-daily institutional use the maintenance-key-protected administrative endpoint is the recommended path.
• Cache-key edge cases. The per-call cache is keyed on the request epoch in TT plus the active model selection; pathological floating-point edge cases at the rounding boundary may produce a cache miss instead of a hit. Rare and harmless.
• Concurrent IERS refreshes. If the automated daily refresh fires while a request is in flight, the post-refresh request will observe a different data identifier than a request submitted moments earlier. The Reproducibility ID embeds the data identifier so the discrepancy is detectable.

12.8 Validation status

The validation suite is a work in progress. The current state, as of 2026-06-14, is summarised in the table below. Tests verify the contracts they document; broader baselines required to support sub-arcsecond precision claims are listed under "pending baselines". A live, machine-generated mirror of this section that reads the on-disk fixtures and the test directory at request time is published at /calculadora-lunar-cientifica/validation-report; that page is the source of truth for the current measured state.

Pending baselines (not yet implemented). The following extensions are planned and are required before the methodology can defend any sub-arcsecond precision claim. They are listed here so the gap is visible to the reader.

• JPL Horizons regression with 200 to 1000 samples covering 1900-2100, multiple observer latitudes (including polar cases), and stratified by lunar phase;
• USNO Astronomical Almanac cross-check for the canonical mid-latitude ephemeris quantities;
• IMCCE numerical comparison (INPOP21a as the published independent ephemeris) over the DE440 native window;
• Round-trip unit tests across the unit conversions exposed by the form (km / AU / mi, deg / dms / rad, UTC / JD / MJD).

Current measured residual. Spot measurements against the limited Horizons smoke fixtures, in the canonical DE440 + IAU 2006 + IAU 2000A configuration at typical mid-latitude epochs, indicate an apparent-direction residual of order a few arcseconds and a geocentric-distance residual of order a few kilometres. These numbers are consistent with the documented per-source uncertainty budget of Section 11; they are not the millisecond-of-arc figure historically associated with operational reference services such as JPL Horizons or IMCCE. Any future tightening of the published precision target will be added to Section 11.10 only when a matching fixture suite is in place and is checked into the repository.

Honest framing. The combination implemented here (DE440 numerical ephemeris, the IAU 2006 P03 precession, IAU 2000A nutation, IERS 2010 frame chain, GUM 100:2008 uncertainty propagation, CCSDS 502.0-B-3 export) is the same set of standards adopted by the operational reference services. What this methodology does not claim, until the broader baseline above is in place, is observational equivalence with those services at the millisecond-of-arc level. The Reproducibility ID and the JSON schema make any forthcoming validation result citable bit-for-bit against the engine version that produced it.

12.9 Reproducibility via Live Permalink

Implications of the live permalink mechanism (Section 12.1):

• Refreshing the page preserves the exact computational state (the URL contains the full parameter set);
• Sharing a URL (copy from the address bar or via the dedicated button) reproduces the entire scenario, including all Group 1 / Group 2 / Group 3 toggles and theoretical-lab parameters;
• Citation in academic papers can include the URL with hash; the hash is computable from the parameter list and serves as an integrity check against later edits;
• The full 64-character SHA-256 is available via the JavaScript Web Crypto API and can be embedded in a citation alongside the truncated 16-character form when the higher collision resistance is required.

13. Complete changelog

Each entry below records numerical changes to the calculator's scientific computation or new scientific features added. Dates in UTC.

13.1 Version 86.9.16 (2026-05-09)

• Monte Carlo Box-Muller cache reset on seed change so successive runs with the same seed produce bit-equivalent samples (preserves Reproducibility ID guarantee).
• Numeric inputs now reject non-finite values (INF / overflow strings like 1e500) instead of cascading NaN into the engine.
• Apsides parabolic refinement clamped to [-1, +1] so the refined timestamp stays inside the bracket; super/micro Moon distances no longer drift on degenerate triples.
• Eclipse central-duration returns null when the sub-solar point falls off the disc, instead of fabricating a duration from a zero-latitude shadow speed.
• Lunar occultation contact refinement replaces the 1-hour fallback with a bisective expansion (up to 4 doublings) so the emersion is located honestly.
• Bennett 1982 refraction now delegates to Saemundsson 1986 below 5 deg altitude (instead of returning zero), preserving rise/set refraction near the horizon.
• Angular-separation between the Moon and a star uses haversine (2 asin sqrt(a)) instead of spherical law of cosines, preserving sub-arcsec precision required for occultation contact.

13.2 Version 86.9.15 (2026-05-09)

• Principal-phase root finder widened from a 1-day half window to 3 days plus bilateral expansion (up to 4 steps), so First Quarter no longer disappears from the next-events list and Full Moon dates are not skipped past valid cycles. Affects all consumers of the canonical phase engine (calculators, super-Moon search, eclipse pairing, lunar-fase).

13.3 Version 86.9.14 (2026-05-09)

• Popular calculator output extended with Bortle-aware limiting magnitude (Schaefer 1990 plus 5 log10(D/7) aperture term), 24-hour visibility window (above user-set minimum altitude), distance reported in km, Earth radii, lunar diameters and AU, orbital velocity in km/s, apparent diameter percent versus mean (31.1 arcmin), and synodic-month full-Moon name.
• Astrophotography calculator output extended with exposure-motion limits (rule of 500, lunar tracking residual 0.549 arcsec/s, Nyquist sampling check), framing fit (Moon count horizontal/vertical, recommended focal length, ND filter suggestion by apparent magnitude), stacking (SNR sqrt(N), dB equivalent, dynamic range 6.02 bits in dB) and Bortle-aware sky brightness. Citations: Pickering 2002 (Looney 11), Dawes 1867, Bortle 2001, Schaefer 1990, Howell 2006.

13.4 Version 86.8.4 (2026-05-09)

• Eclipse central-duration derivative now uses the same engine as the maximum sample.
• Bennett 1982 refraction implemented natively (was silently delegating to Saemundsson 1986).
• Annual lunar aberration corrected to kappa = 20.49552 arcsec with full Meeus 23.1 / 23.2 (was 20.4898 arcsec, degenerate form).
• Monte Carlo triangular sampling fixed to T = a*(u1+u2-1), Var = a^2/6 (was a^2/24).
• GUM t-Student lookup uses linear interpolation (was nearest-up snap).
• GUM t-Student coverage table densified at nu = 4, 25, 40, 75, 100; max interpolation error below 0.005.
• Refraction zenith threshold tightened from 89.9 deg to 90.0 deg.
• All four refraction models default to 10 deg C standard atmosphere.
• Delta-T Espenak-Meeus polynomial: 10-year linear blend zone 2045-2055 across the year=2050 boundary discontinuity.
• Time-scale output adds delta_t_actual_seconds from canonical chain (Espenak-Meeus estimate kept as separate field).
• Apparent corrections output flags omission of the Meeus 23.1 eccentricity term in Lite mode (~0.34 arcsec, included in Full mode).

13.5 Version 86.8.1 (2026-05-09)

• Nutation IAU 1980 form value remapped to actually reach the Wahr 1981 path (was falling back to Meeus 62-term legacy; ~5.4 arcsec dPsi impact at sample epoch).

13.6 Version 86.7.0 (2026-05-08)

• NAIF DAF binary orientation reader (Type 2 Chebyshev).
• Lense-Thirring secular nodal precession module (~0.001 mas/yr; Iorio 2002, Mashhoon 2001).
• Tycho-2 catalogue reader (2.5 million stars).
• Gaia DR3 Gmag<13 subset reader (18 declination bands).
• Type A and Type B input quantities (normal, rectangular, triangular, U-shape).
• JCGM 100:2008 propagation: independent and correlated combination, Welch-Satterthwaite, t-Student, U = k * u_c.
• JCGM 101:2008 Annex C adaptive Monte Carlo (4 distributions, Box-Muller).
• CCSDS 502.0-B-3 OEM v3.0 export (META, 6D state, 21-element covariance).
• Allan deviation analysis (NIST SP 1065).

13.7 Version 86.6.4 (2026-05-08)

• Mendes-Pavlis 2004 refraction mapping function (J. Geophys. Res. 109, B07404; hydrostatic + wet decomposition).
• Lunar occultation predictor (Meeus 1998 chapter 34, two-stage refinement to about 1 minute precision).

13.8 Version 86.6.2 (2026-05-08)

• Earth Rotation Angle module (IAU 2006 + CIO-based GST + Equation of Origins polynomial).
• TT-TDB module (reduced Fairhead-Bretagnon 1990, ~10 microsecond precision over 1900-2100).
• Annual aberration module (kappa = 20.49552 arcsec, Meeus 1998 ch. 23).
• Diurnal aberration module (k_d = 0.3200 * cos(lat) arcsec, Meeus 23.3).
• Alternate IAU nutation/precession models (IAU 2000B, IAU 1980 / Wahr 1981, IAU 1976) unified via dispatcher.
• Time-scales output exposes ERA, IAU 2006 GST, TDB-TT, JD(TDB).

13.9 Version 86.6.0 (2026-05-08)

• Reproducibility identifier format: deterministic readable slug-key replaces hash-only mechanism.

14. Bibliography

The bibliography is organised in two layers. The first layer is the canonical core bibliography (12 entries) shared with the JSON-LD graph emitted in the page head. The second layer enumerates the additional citations referenced in the text above (relativistic stack, observability, uncertainty quantification, station displacement, IERS bulletins, lunar coordinate time). DOIs are provided when available; missing DOIs reflect pre-DOI publications (Bretagnon and Francou 1988, Meeus 1998 monograph).

14.1 Core bibliography

1. Chapront-Touze, M., & Chapront, J. (1988). ELP 2000-85: a semi-analytical lunar ephemeris adequate for historical times. Astronomy and Astrophysics, 190(1-2), 342-352.
2. Bretagnon, P., & Francou, G. (1988). Planetary theories in rectangular and spherical variables: VSOP87 solution. Astronomy and Astrophysics, 202, 309-315.
3. Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. (2021). The JPL Planetary and Lunar Ephemerides DE440 and DE441. The Astronomical Journal, 161(3), 105. DOI: 10.3847/1538-3881/abd414
4. Konopliv, A. S., Park, R. S., Yuan, D. N., Asmar, S. W., Watkins, M. M., Williams, J. G., Fahnestock, E., Kruizinga, G., Paik, M., Strekalov, D., Harvey, N., Smith, D. E., & Zuber, M. T. (2013). The JPL lunar gravity field to spherical harmonic degree 660 from the GRAIL Primary Mission. Journal of Geophysical Research: Planets, 118(7), 1415-1434. DOI: 10.1002/jgre.20097
5. Capitaine, N., Wallace, P. T., & Chapront, J. (2003). Expressions for IAU 2000 precession quantities. Astronomy and Astrophysics, 412(2), 567-586. DOI: 10.1051/0004-6361:20031539
6. Petit, G., & Luzum, B. (Eds.) (2010). IERS Conventions (2010). IERS Technical Note 36, Verlag des Bundesamts fuer Kartographie und Geodaesie.
7. Wahr, J. M. (1981). The forced nutations of an elliptical, rotating, elastic and oceanless Earth. Geophysical Journal of the Royal Astronomical Society, 64(3), 705-727. DOI: 10.1111/j.1365-246X.1981.tb02691.x
8. International Astronomical Union (1991). Resolution A4: Recommendations from the Working Group on Reference Systems. Proceedings of the XXIst IAU General Assembly, Buenos Aires.
9. International Astronomical Union (2000). Resolution B1.9: Re-definition of Terrestrial Time (TT). Proceedings of the XXIVth IAU General Assembly, Manchester.
10. International Astronomical Union (2006). Resolution B3: Re-definition of Barycentric Dynamical Time (TDB). Proceedings of the XXVIth IAU General Assembly, Prague.
11. Klioner, S. A. (2008). Relativistic scaling of astronomical quantities and the system of astronomical units. Astronomy & Astrophysics, 478(3), 951-958. DOI: 10.1051/0004-6361:20077786
12. Espenak, F., & Meeus, J. (2006). Five Millennium Canon of Solar Eclipses: -1999 to +3000. NASA Technical Publication TP-2006-214141.
13. Meeus, J. (1998). Astronomical Algorithms (2nd ed.). Willmann-Bell, Richmond, Virginia.
14. IAU SOFA Board (2021). IAU SOFA Software Collection: standards of fundamental astronomy. International Astronomical Union, http://www.iausofa.org.
15. International Earth Rotation and Reference Systems Service (2024). IERS Bulletin A: rapid service / prediction of UT1-UTC and polar motion. U.S. Naval Observatory, weekly issues, https://www.iers.org.
16. Bennett, G. G. (1982). The calculation of astronomical refraction in marine navigation. The Journal of Navigation, 35(2), 255-259. DOI: 10.1017/S0373463300022037
17. Saemundsson, T. (1986). Astronomical refraction. Sky and Telescope, 72, 70.
18. Prince, P. J., & Dormand, J. R. (1981). High order embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 7(1), 67-75. DOI: 10.1016/0771-050X(81)90010-3
19. Hairer, E., Norsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed.). Springer-Verlag, Berlin Heidelberg.
20. Tapley, B. D., Schutz, B. E., & Born, G. H. (2004). Statistical Orbit Determination. Academic Press, Burlington, MA.
21. Crassidis, J. L., & Junkins, J. L. (2011). Optimal Estimation of Dynamic Systems (2nd ed.). CRC Press, Boca Raton, FL.
22. Gordon, N. J., Salmond, D. J., & Smith, A. F. M. (1993). Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F (Radar and Signal Processing), 140(2), 107-113.
23. Arulampalam, M. S., Maskell, S., Gordon, N., & Clapp, T. (2002). A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50(2), 174-188. DOI: 10.1109/78.978374
24. Kitagawa, G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. Journal of Computational and Graphical Statistics, 5(1), 1-25.
25. Acton Jr., C. H. (1996). Ancillary data services of NASA's Navigation and Ancillary Information Facility. Planetary and Space Science, 44(1), 65-70. DOI: 10.1016/0032-0633(95)00107-7
26. NAIF (Navigation and Ancillary Information Facility) (2024). SPICE Toolkit Documentation: DAF Required Reading, SPK Required Reading. NASA Jet Propulsion Laboratory, https://naif.jpl.nasa.gov.
27. Joint Committee for Guides in Metrology (2008). Evaluation of measurement data — Supplement 1 to the "Guide to the expression of uncertainty in measurement" — Propagation of distributions using a Monte Carlo method (JCGM 101:2008). Bureau International des Poids et Mesures (BIPM).
28. Murphy, T. W. (2013). Lunar laser ranging: the millimeter challenge. Reports on Progress in Physics, 76(7), 076901. DOI: 10.1088/0034-4885/76/7/076901
29. Viswanathan, V., Fienga, A., Minazzoli, O., Bernus, L., Laskar, J., & Gastineau, M. (2018). The new lunar ephemeris INPOP17a and its application to fundamental physics. Monthly Notices of the Royal Astronomical Society, 476(2), 1877-1888. DOI: 10.1093/mnras/sty096
30. International Laser Ranging Service (2024). Consolidated Ranging Data (CRD) Format Specification v2.01. ILRS, https://ilrs.gsfc.nasa.gov.
31. Pavlis, N. K., Holmes, S. A., Kenyon, S. C., & Factor, J. K. (2012). The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research: Solid Earth, 117(B4), B04406. DOI: 10.1029/2011JB008916
32. Pijpers, F. P. (1998). Helioseismic determination of the solar gravitational quadrupole moment. Monthly Notices of the Royal Astronomical Society, 297(3), L76-L80.
33. de Moura, L., & Ullrich, S. (2021). The Lean 4 theorem prover and programming language. Automated Deduction – CADE 28, LNCS 12699, 625-635. DOI: 10.1007/978-3-030-79876-5_37
34. The mathlib Community (2020). The Lean mathematical library. Proceedings of the 9th ACM SIGPLAN International Conference on Certified Programs and Proofs (CPP 2020), 367-381. DOI: 10.1145/3372885.3373824
35. Sicilia, M. A., Garcia-Barriocanal, E., & Sanchez-Alonso, S. (2017). Community curation in open dataset repositories: insights from Zenodo. Procedia Computer Science, 106, 54-60. DOI: 10.1016/j.procs.2017.03.009

14.2 Extended bibliography

1. Wallace P. T. and Capitaine N., 2006, Earth Rotation Angle and the IAU 2006 precession-nutation framework. Astronomy and Astrophysics, 459, 981. DOI: 10.1051/0004-6361:20065897.
2. McCarthy D. D. and Luzum B. J., 2003, An abridged model of the precession-nutation of the celestial pole (IAU 2000B). Celestial Mechanics and Dynamical Astronomy, 85, 37.
3. Wahr J. M., 1980, The forced nutations of an elliptical, rotating, elastic and oceanless Earth. Geophysical Journal of the Royal Astronomical Society, 64, 705.
4. Fairhead L. and Bretagnon P., 1990, An analytical formula for the time transformation TB-TT. Astronomy and Astrophysics, 229, 240.
5. Hilton J. L. and Hohenkerk C. Y., 2002, Rotational elements of the planets and satellites: a comparison of FB91 to TE405. IAU Joint Discussion 16, Session 5.
6. Mendes V. B. and Pavlis E. C., 2004, High-accuracy zenith delay prediction at optical wavelengths. Geophysical Research Letters, 31, L14602. DOI: 10.1029/2004GL020308.
7. Hog E. et al., 2000, The Tycho-2 catalogue of the 2.5 million brightest stars. Astronomy and Astrophysics, 355, L27.
8. Gaia Collaboration, 2023, Gaia Data Release 3: Summary of the content and survey properties. Astronomy and Astrophysics, 674, A1. DOI: 10.1051/0004-6361/202243940.
9. ESA, 1997, The Hipparcos and Tycho catalogues. ESA SP-1200.
10. Murphy T. W., Nordtvedt K. and Turyshev S. G., 2007, Gravitomagnetic influence on gyroscopes and on the lunar orbit. Physical Review Letters, 98, 071102. DOI: 10.1103/PhysRevLett.98.071102.
11. Iorio L., 2002, On the possibility of measuring the Lense-Thirring effect with a Sun-orbiting spacecraft. Classical and Quantum Gravity, 19, 5473. DOI: 10.1088/0264-9381/19/21/313.
12. Klioner S. A., 2003, A practical relativistic model for microarcsecond astrometry in space. Astronomical Journal, 125, 1580. DOI: 10.1086/367593.
13. Shapiro I. I., 1964, Fourth test of general relativity. Physical Review Letters, 13, 789. DOI: 10.1103/PhysRevLett.13.789.
14. Bertotti B., Iess L. and Tortora P., 2003, A test of general relativity using radio links with the Cassini spacecraft. Nature, 425, 374. DOI: 10.1038/nature01997.
15. de Sitter W., 1916, On Einstein's theory of gravitation, and its astronomical consequences. MNRAS, 77, 155.
16. Williams J. G., Turyshev S. G. and Boggs D. H., 2004, Progress in lunar laser ranging tests of relativistic gravity. Physical Review Letters, 93, 261101. DOI: 10.1103/PhysRevLett.93.261101.
17. Ashby N. and Patla B. R., 2024, A relativistic framework for clocks and astronomical objects on the Moon. Astronomical Journal, 168, 112.
18. Pound R. V. and Rebka G. A., 1960, Apparent weight of photons. Physical Review Letters, 4, 337.
19. Vessot R. F. C. et al., 1980, Test of relativistic gravitation with a space-borne hydrogen maser. Physical Review Letters, 45, 2081.
20. Everitt C. W. F. et al., 2011, Gravity Probe B: final results of a space experiment to test general relativity. Physical Review Letters, 106, 221101.
21. Ciufolini I. and Pavlis E. C., 2004, A confirmation of the general relativistic prediction of the Lense-Thirring effect. Nature, 431, 958.
22. Nordtvedt K., 1968, Equivalence principle for massive bodies. II. Theory. Physical Review, 169, 1014.
23. Hofmann F. and Mueller J., 2018, Relativistic tests with lunar laser ranging. Classical and Quantum Gravity, 35, 035015.
24. Williams J. G., Turyshev S. G. and Boggs D. H., 2012, Lunar laser ranging tests of the equivalence principle. Classical and Quantum Gravity, 29, 184004.
25. Kostelecky V. A. and Russell N., 2011, Data tables for Lorentz and CPT violation. Reviews of Modern Physics, 83, 11.
26. Abbott B. P. et al. (LIGO and Virgo), 2017, GW170104: observation of a 50-solar-mass binary black hole coalescence at redshift 0.2. Physical Review Letters, 119, 161101.
27. Eckhardt D. H., 1981, Theory of the libration of the Moon. The Moon and the Planets, 25, 3.
28. Williams J. G., Konopliv A. S., Boggs D. H. et al., 2014, Lunar interior properties from the GRAIL mission. Journal of Geophysical Research Planets, 119, 1546.
29. Murphy T. W. et al., 2010, Long-term degradation of optical devices on the Moon. Icarus, 211, 1103.
30. JCGM, 2008. Evaluation of measurement data - Guide to the expression of uncertainty in measurement. JCGM 100:2008 (GUM).
31. JCGM, 2008. Evaluation of measurement data - Supplement 1 to the GUM - Propagation of distributions using a Monte Carlo method. JCGM 101:2008.
32. JCGM, 2011. Evaluation of measurement data - Supplement 2 to the GUM - Extension to any number of output quantities. JCGM 102:2011.
33. CCSDS, 2023. Orbit Data Messages, Recommended Standard CCSDS 502.0-B-3. Consultative Committee for Space Data Systems.
34. Riley W. J., 2008. Handbook of Frequency Stability Analysis. NIST Special Publication 1065.
35. Petit G. and Luzum B. (eds), 2010. IERS Conventions (2010). IERS Technical Note 36, Verlag des Bundesamts fur Kartographie und Geodasie.
36. International Earth Rotation and Reference Systems Service, 2024. IERS Bulletin A: rapid service / prediction of UT1-UTC and polar motion. U.S. Naval Observatory weekly issues.
37. Bizouard C., Lambert S., Gattano C., Becker O. and Richard J. Y., 2019, The IERS EOP 14C04 solution for Earth orientation parameters consistent with ITRF 2014. Journal of Geodesy, 93, 621.
38. Stephenson F. R., Morrison L. V. and Hohenkerk C. Y., 2016, Measurement of the Earth's rotation: 720 BC to AD 2015. Proceedings of the Royal Society A, 472, 20160404.
39. Hapke B., 1984, Bidirectional reflectance spectroscopy: 3. Correction for macroscopic roughness. Icarus, 59, 41.
40. Pickering K. A., 2002, The southern limits of the ancient star catalog. DIO, 12, 1.
41. Lambert S. B., 2007, Empirical model of the Free Core Nutation. Astronomy and Astrophysics, 471, 1097.
42. Irwin A. W. and Fukushima T., 1999, A numerical time ephemeris of the Earth. Astronomy and Astrophysics, 348, 642.
43. Acton C. H., 1996, Ancillary data services of NASA's Navigation and Ancillary Information Facility. Planetary and Space Science, 44, 65 (SPICE introduction).
44. Saastamoinen J., 1972, Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. Bulletin Geodesique, 105, 279.

15. Glossary

Definitions of the technical terms used throughout the document. The same glossary is exposed in machine-readable form via the schema:DefinedTermSet JSON-LD graph emitted in the page head.

ICRS

International Celestial Reference System; the kinematically non-rotating reference system aligned with extragalactic radio sources (IAU 1997).

ICRF

International Celestial Reference Frame; realisation of ICRS by quasar positions (IERS).

GCRS

Geocentric Celestial Reference System; geocentric counterpart of ICRS (IERS Conventions 2010 Section 5.3).

ITRS

International Terrestrial Reference System; Earth-fixed, co-rotating frame (IERS Conventions 2010 Chapter 4).

CIO

Celestial Intermediate Origin; reference origin used in the new "CIO-based" Earth-rotation transformation (Capitaine 2003).

CIP

Celestial Intermediate Pole; instantaneous pole of Earth rotation (IAU 2000).

ERA

Earth Rotation Angle; the angle between the CIO and the Terrestrial Intermediate Origin, linear function of UT1.

UT1

Universal Time, Earth-rotation-based; differs from UTC by DUT1 (less than 0.9 s by definition).

TAI

International Atomic Time; the proper-time scale of the SI second realised at the geoid.

TT

Terrestrial Time; TT = TAI + 32.184 s; the independent variable for terrestrial ephemerides.

TDB

Barycentric Dynamical Time; differs from TT by periodic terms below ~2 ms.

Delta-T

Delta-T = TT - UT1; tabulated by IERS, modelled by Espenak and Meeus 2006 outside the IERS coverage window.

ELP-2000/82B

Semi-analytical lunar ephemeris by Chapront-Touze and Chapront (1988); a compact truncation retaining the dominant signal terms is used in the default mode.

VSOP87D

Planetary theory by Bretagnon and Francou (1988) in heliocentric ecliptic spherical coordinates of date.

DE440

JPL Development Ephemeris 440 (Park et al. 2021); covers 1550-2650; used as a higher-precision option when the corresponding numerical-ephemeris data is available locally.

NPB

Frame-Bias x Precession x Nutation matrix; transforms vectors from GCRS to the true equator and equinox of date.

IAU 2006 P03

Precession model of Capitaine, Wallace and Chapront (2003), adopted by the IAU in 2006.

IAU 2000A

High-precision nutation model of Mathews, Herring and Buffett (2002), expressed as a long luni-solar plus planetary series.

SOFA

Standards of Fundamental Astronomy: the IAU-endorsed reference implementation of these models.

GUM

JCGM 100:2008, the Guide to the Expression of Uncertainty in Measurement; canonical reference for combined uncertainty u_c via RSS with sensitivity coefficients.

CCSDS OEM

CCSDS 502.0-B-3 Orbit Ephemeris Message; standard format for spacecraft and natural-body orbit data exchange.

16. Advanced Layer 5 modules (interactive)

Layer 5 is the interactive research tier sitting on top of the production ephemeris pipeline documented in sections 1 to 12. Each module is a self-contained tile that exposes a working numerical method, a deterministic worked example, and an honest statement of what is and is not implemented. The modules are reachable from the scientific calculator UI (cientifica.php) and cross-link back into this methodology page through the anchors listed below.

16.1 Cowell propagator

The Cowell tile integrates the Moon's geocentric Cartesian state with a Dormand-Prince RK8(7) 13-stage embedded pair (Prince and Dormand 1981; Hairer, Norsett and Wanner, Solving Ordinary Differential Equations I, Table 5.2), using adaptive step control driven by the embedded error estimate. The force model is Earth point-mass acceleration plus the Earth J2 zonal term taken from EGM2008, plus the Sun as a third body in the standard direct-plus-indirect Cowell formulation. Sun coordinates are produced by the low-precision Meeus analytic series (accuracy near 0.5 degree), which is sufficient for the 24-hour propagation drift demonstration the tile is built around. A pure two-body 24-hour run shows an energy drift near 4.9e-15 (machine precision); enabling the solar third body produces a tidal position drift near 58 km over the same interval, consistent with the analytic estimate. The UI exposes three tolerance presets (1e-9, 1e-12, 1e-14).

Not implemented: planetary perturbations from Venus, Mars, Jupiter and Saturn; post-Newtonian corrections inside the integrator (the PN stack is exposed as a separate Layer 5 tile and as section 8 in this document); ocean and solid-Earth tide accelerations; lunar libration coupling.

16.2 Kalman filter

Two filters share the same matrix backend. The linear KF runs the textbook prediction x' = F x, P' = F P F^T + Q and update K = P H^T (H P H^T + R)^-1 with the Joseph-free standard form (I - K H) P. The extended KF uses a nonlinear state transition f(x, dt) with a Jacobian F propagating the covariance. Two demos are wired in: a one-dimensional constant-velocity tracker and a two-dimensional circular orbit observed only in range. References: Tapley, Schutz and Born, Statistical Orbit Determination; Crassidis and Junkins, Optimal Estimation of Dynamic Systems. Sample run with sigma_z = 1.0 produces a constant-velocity RMS near 0.497 (about 50% noise reduction over the raw measurements); the EKF orbit demo with sigma_R = 0.05 gives a position RMS near 0.115 over 100 steps. Matrix inversion uses Gauss-Jordan with partial pivoting; singular systems throw RuntimeException.

Not implemented: Unscented Kalman Filter (UKF), square-root covariance form, adaptive process-noise tuning of Q.

16.3 Particle filter

A bootstrap Sampling Importance Resampling (SIR) particle filter with systematic resampling (Kitagawa 1996) using the canonical stratified draw u_i = (i - 1 + u_0) / N with u_0 ~ U[0, 1/N). Effective sample size ESS = 1 / sum(w_i^2) drives an automatic resample when ESS < N * threshold (default 0.5 N). The wired benchmark is the Gordon-Salmond-Smith 1993 nonlinear test x_{k+1} = 0.5 x + 25 x / (1 + x^2) + 8 cos(1.2 k) + w, z = x^2 / 20 + v. References: Arulampalam et al. 2002, IEEE Transactions on Signal Processing 50:174; Doucet and Johansen 2009. Sample RMS at N = 2000 is about 4.11, inside the literature range of 4 to 6.

Not implemented: auxiliary particle filter, Rao-Blackwellised particle filter, regularised resampling.

16.4 SPICE kernel writer

The kernel writer emits a NAIF DAF binary SPK of type 9 (Lagrange polynomial, unequal time steps). The DAF architecture follows the NAIF spec: a 1024-byte file record, an optional comment area, summary records, name records, and element records, all encoded as little-endian IEEE 754 doubles. The 28-byte FTP-string validation sequence from the NAIF DAF Required Reading is written verbatim: FTPSTR:\r:\n:\r\n:\r\x00:\x81:\x10\xce:ENDFTP. Round-trip validation against NAIF CSPICE via spiceypy succeeds: after furnsh, spkobj reports the target body, spkcov reports the correct coverage window, and spkez recovers the Moon state at the midpoint epoch. A 24-hour Moon kernel writes to 4096 bytes with the recovered geocentric distance equal to 384,400 km at midpoint. Kernels are served from /lua/calculadora-lunar/_download-bsp.php with SHA validation and auto-cleanup after one hour. Reference: NAIF "DAF Required Reading" and the "spk.req" type 9 specification.

Not implemented: Type 13 (Hermite), Type 18 (variable-degree Lagrange/Hermite), PCK orientation kernels, populated comment area.

16.5 Numerical witnesses for Lean/Coq proofs

Two Lean 4 / Mathlib draft files live under core/astro/proofs/lean/: NpbDeterminant.lean states the theorem that the determinant of the NPB chain equals 1, with a proof skeleton using Matrix.det_mul applied to the product of three rotation determinants; and RotationOrthogonality.lean states the orthogonality of the elementary rotation matrices R_x, R_y, R_z. Both proof bodies are sorry placeholders; rigorous human-readable proofs are written out in PAPER_PROOFS.md as Theorem 1 (NPB determinant equals 1), Theorem 2 (Kepler's Third Law from Newtonian gravity), Theorem 3 (IAU rotation orthogonality), and Theorem 4 (ELP-2000 Poisson series convergence via the Weierstrass M-test). The PHP layer produces numerical witnesses: the NPB determinant computed from the IAU 2006 / 2000A matrices is expected to equal 1.0 to within 1e-16, and the rotation orthogonality residual ||R R^T - I||_F stays below 1e-14.

Not implemented: full Lean typecheck inside the toolchain (requires elan plus a Mathlib install of approximately 2 GB).

16.6 Monte Carlo uncertainty propagation (JS client-side)

A pure JavaScript Web Worker performs Monte Carlo uncertainty propagation entirely client-side, with no WebAssembly dependency. Normal samples come from the Box-Muller transform; uniform samples are drawn directly for inputs declared as uniform. Sample sizes of 10,000, 100,000, 500,000 and 1,000,000 are selectable from the UI. The output is a 50-bin histogram and the percentiles p2.5, p50 and p97.5. The propagated formula is provided as an editable JavaScript function body, function(x) { ... }, where x is the input record. Progress is streamed back via postMessage every 10 percent. Reference: JCGM 101:2008, Evaluation of Measurement Data — Supplement 1 to the GUM, Propagation of Distributions Using a Monte Carlo Method.

Not implemented: optional WebAssembly kernel for speedup; correlation matrix between inputs (current version assumes independence).

16.7 LLR residuals and CRD parser

The Lunar Laser Ranging tile predicts the round-trip light time 2 r / c for a given station-retroreflector pair, with a nominal value near 2.564 seconds. The station catalogue ships APOLLO (Apache Point), OCA (Cote d'Azur, Grasse), McDonald and Matera, each with ILRS coordinates. The retroreflector catalogue ships A11, A14, A15 (Apollo), L17 (Lunokhod 2, Le Monnier 1973) and L21 (Lunokhod 1, Mare Imbrium 1970), each with selenodetic coordinates. Typical literature RMS residuals are 1.5 cm for APOLLO (Murphy 2013) and 2.5 cm for OCA (Viswanathan 2018). The CRD v2.01 parser accepts H1-H4 header records, record type 10 (full-rate), record type 11 (normal points), and record type 99 (end of session); H8, H9 and configuration records are skipped silently. A synthesizer produces realistic CRD streams with optional Gaussian noise (Box-Muller, sigma expressed in centimetres). Documented error budget contributions: atmosphere near 10 cm, retroreflector array signature near 1 cm, Shapiro relativistic delay near 10 ns. Reference: ILRS Consolidated Laser Ranging Data (CRD) format specification at ilrs.gsfc.nasa.gov.

Not implemented: ingestion of real ILRS data files (requires NASA Earthdata credentials); orbit fitting from residuals.

16.8 Citation snapshot / Zenodo DOI

The citation tile assembles a Zenodo Deposition API metadata payload as JSON and constructs a dry-run curl command for the user to inspect. Actual upload requires a user-provided ZENODO_ACCESS_TOKEN exported into the shell; the methodology page never holds the token. The APA-format bibliography rendered by Layer 4 (sections 14.1 and 14.2 of this document) is reused as the citation block. Reference: Zenodo Deposition REST API documentation.

Not implemented: automatic publishing of a new DOI (requires the user's Zenodo token and is intentionally left as a human-in-the-loop step).

16.9 Engine quick selector

At the top of the scientific calculator form, a three-way radio control switches between Auto, Lite and DE440 engines. Auto selects DE440 when the SPK kernel is available on disk and silently falls back to Lite otherwise. Lite uses ELP-2000/82B truncated for the Moon, VSOP87D for the Sun and inner planets, and the SOFA polyfill for frame transformations, with no kernel files required. DE440 reads the NAIF SPK kernel through the custom PHP DAF reader described in section 4.1. The active state is rendered using the CSS :has(input[type="radio"]:checked) selector, with no JavaScript needed.

17. Honest limitations of the calculator

This section consolidates the things the calculator explicitly does not do, complementing the per-module "Not implemented" notes in section 16 and the scope limits in section 12.7. The intent is to keep the methodology calibrated to what production actually ships, so that downstream citations stay honest.

• Automatic DOI publication. The Zenodo tile (16.8) builds the deposition payload and the curl command, but never uploads. Publishing is a deliberate human-in-the-loop step requiring the user's own Zenodo access token.
• WebAssembly Monte Carlo kernel. The Monte Carlo tile (16.6) is pure JavaScript inside a Web Worker. A WASM kernel is documented as an optional upgrade path but is not shipped.
• Lean / Coq machine-checked proofs. The Lean files (16.5) carry sorry proof bodies; the rigorous proofs live in PAPER_PROOFS.md and are reproduced numerically in PHP, but no machine typecheck is performed in the build.
• Real ILRS data ingestion. The LLR tile (16.7) parses synthetic CRD streams it generates itself. Real NASA Earthdata-gated ILRS streams require credentials and are not pulled.
• Planetary perturbations in the Cowell propagator. The Layer 5 Cowell tile (16.1) models Earth point-mass plus J2 plus the Sun as a third body. Venus, Mars, Jupiter and Saturn, ocean and solid-Earth tides, and post-Newtonian corrections inside the integrator are out of scope.
• UKF, square-root KF and adaptive Q. The Kalman tile (16.2) ships only the linear KF and the EKF with the textbook standard-form covariance update.
• 3D visualisation. All result blocks are tabular or 2D static plots. No WebGL or three.js scene is rendered, and no animated globe or selenographic globe is produced.

Found a methodology issue? Want to cite this work?

Email rcgwebsites@gmail.com for:

• Errors in algorithm description, equations, or references
• Missing citations or methodology gaps
• Validation data against JPL Horizons / SPICE / Skyfield (we welcome cross-validation reports)
• Citation requests (BibTeX, APA, plaintext blocks)
• Suggestions for additional algorithmic models (DE441, asteroid perturbations, IERS daily updates)

We reply within ~5 business days. For peer-review collaboration or independent code audit, please indicate institutional affiliation.