Conversation
…ynamicMassFit computeFragEndParams() already returns the final lat/lon/height; it now also returns the apparent ground-fixed radiant azimuth and elevation there, which dark flight needs as input. The local horizon rotates along the path, so the values at the reference point are not valid at the end. The derotation follows SampleTrajectoryPositions and Orbit.calcOrbit, without a J2000 -> epoch-of-date precession step: the trajectory ECI coordinates are already in the epoch of date. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…tion The trajectory ECI coordinates are in the epoch of date (stations from geo2Cartesian, lines of sight from altAz2RADec) and raDec2AltAz expects epoch-of-date coordinates. Precessing the derotated radiant from J2000 to the date applied precession a second time, biasing the ground-fixed azimuth/elevation by the general precession accumulated since J2000 (~50"/yr). Without that step the output matches the apparent ground-fixed radiant reported by Orbit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ssFit Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
runFragSim() received gamma_a but ignored it, propagating the fragment with a hardcoded Gamma*A = 0.7*1.21 = 0.847 while the dynamic mass was computed with gamma_a (default 0.55). The simulation therefore started with a deceleration 0.847/gamma_a times the measured one. It now uses gamma = 1 and shape_factor = gamma_a, as dynamicMass() does. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…odel The covariance was built from least_squares' res.jac, which is the Jacobian reweighted by the soft_l1 loss, not the Jacobian of the linear model, and then scaled by the robust (MAD) residual variance. This inflated the deceleration uncertainty (by ~25% in a test case), and with it the +/-2 sigma dynamic mass interval, since the mass scales as deceleration^-3. It now uses the linear model Jacobian. Also drops a duplicated MAD = 0 fallback. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
fitVelocity() already starts with its own robust fit, so the preceding _robust_linear_fit() call only provided a starting point for an identical fit. The result is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The outlier plot label said 5 sigma regardless of --sigma_clip (default 3). The sigma_clip docstring described a slope-uncertainty threshold, but the clipping uses the robust (MAD) residual standard deviation. The --maxvel help said there is no filter by default, but velocities above 73 km/s are always removed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The +/-2 sigma dynamic mass interval only propagates the deceleration uncertainty. Velocity, atmospheric density, bulk density and Gamma*A are not included, although the mass scales as v^6 and rho_atm^3. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
When the evaluation velocity is already below 3 km/s the fragment simulations are skipped, but the plot label and the final summary still used the final mass and coordinates, raising a NameError. They now default to NaN. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The radiant fitted by the trajectory solver is the initial tangent of the path (the solver models gravity as a drop from that line), and drag does not rotate the direction of motion, but gravity steepens it at a rate g*cos(elev)/v. The final elevation is now steepened by that turn, integrated with the average speed over the observed part and the simulated speeds after it. In a test case this adds 0.27 deg, matching an integration over the observed velocity profile to 0.005 deg. The azimuth is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This reverts commit 7808ca1. The fixed Gamma*A in runFragSim() is on purpose. The single-body equations only depend on Gamma*A/m^(1/3), which the dynamic mass fixes from the measured deceleration, so simulating with the same Gamma*A as the dynamic mass gives an identical velocity and path for any gamma_a (only the mass scales as Gamma*A^3). A fixed simulation Gamma*A is what lets gamma_a be calibrated against the observed final velocity. Adds a comment stating this. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The solver models gravity as a drop from the fitted line, so traj.radiant_eci_mini is the tangent of the path at its beginning, and Orbit.calcOrbit() derotates it with v_init. In an end-to-end test (a fireball simulated in ECI with drag against the co-rotating air and gravity, observed from three stations and solved with Trajectory), derotating with v_init gave a mean final-direction error of 0.13 deg over seven geometries, against 0.18 deg with v_avg_norot. The difference is within the solver's own radiant error for slow, long fireballs, and negligible for fast ones. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The solver only models the gravity drop when gravity_correction is set; the fitted radiant is then the path tangent at its beginning. Without it, the straight-line fit absorbs part of the curvature and is closer to the tangent at the middle of the observed path, so integrating the turn from the beginning would count part of it twice. The turn now starts at t = 0 or at the middle of the observations accordingly. End-to-end test with the full solver (including the timing re-run where the gravity drop is applied), seven synthetic fireballs with drag against the co-rotating air and gravity: mean final-direction error 0.081 deg with gravity_correction (v_init derotation, turn from t = 0; 0.071 deg at the start, i.e. the solver's own radiant error), and 0.102 deg without it (turn from the middle, against 0.130 deg from t = 0). v_avg_norot gave 0.110 and 0.120 deg, supporting the previous commit's use of v_init. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…y model The evaluation height came from a straight line fitted to the measured heights versus time. At the centre of the window that line is biased upwards by the curvature of the decelerating path, and it carries the measurement noise. computeFragEndParams() then mapped that height back to a time through the solver's model heights, so the simulation started at a different instant than the velocity it was given. The height is now taken at the evaluation time from the same model interpolation computeFragEndParams() uses (which includes the gravity drop), so the height, velocity and time of the starting point coincide exactly. Synthetic fireballs solved with the full solver: the error of the evaluation height against the true height dropped from +522..+780 m to -19..+37 m for a 15 km window, and from a mean of 338 m to 51 m over all tested cases. In short windows at the very end of slow, long fireballs the error is bounded by the solver's own model heights there (up to ~130 m). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The final point was placed on the fitted straight line, without the gravity drop that the solver applies to its model points (and that the evaluation height now includes), so the fragment started on the curved trajectory and ended on the straight line. When gravity_correction is set, the final point now gets the solver's drop at its time, keeping only the component perpendicular to the line, as in the solver's model heights (its re-projection onto the lines of sight absorbs the along-track part into the length). Checks: at the observed points of real trajectories this reproduces the solver's model heights to <= 10 m in 6 of 7 cases (56 m for a near-vertical one), against up to 195 m without the drop and 118 m with the full vertical drop. Against the true final position of synthetic fireballs solved with the full solver, the mean height error dropped from 319 m to 201 m. The residual comes from the solver's drop model itself, which assumes a constant vertical velocity and overestimates the drop at the end of slow, long fireballs; its own model heights show the same error there. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replaces the hardcoded 9.81 m/s^2 with MetSim's G0 scaled to the final height, G0/(1 + h/r_earth)^2, where most of the turn accumulates (the turn rate goes as 1/v). The effect is below 1% of the turn (~0.001 deg in a test case); it keeps the constant consistent with the simulation used in the same function. Also fixes a missing space. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The length came from interpolating height against length, which mixes the savgol-smoothed heights with the raw lengths. At shallow entry angles a small height difference becomes a large length difference, so the starting length of the fragment, and with it the final point, could be off by hundreds of metres, or kilometres for grazing trajectories. It is now obtained by inverting the trajectory geometry: the length whose point on the fitted line, with the solver's gravity drop, has the evaluation height. The drop moves into pointOnTrajectory(), which is now used both here and for the final point. Checks on seven real trajectories: the geometric length reproduces the solver's own stored lengths at every observed point (median ~0, rms <= 35 m), the inversion recovers the evaluation height to 1e-9 m, and it differs from the interpolated length by -78 to -2130 m near the evaluation region. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
DynamicMassFit.computeFragEndParams()returns the final lat/lon/height of the propagatedfragment. It now also returns the apparent ground-fixed radiant azimuth and elevation at that
point, which dark flight needs as input. Both are printed in the CLI summary for the three mass
cases.
SampleTrajectoryPositions.sampleTrajectory()precessed the derotated radiant from J2000 to theepoch of date before converting it to azimuth/elevation. The radiant is already in the epoch of
date, so precession was applied twice. That step and its imports are removed.
Why the coordinates are already in the epoch of date
geo2Cartesian()with its defaultprecess_j2000=False.altAz2RADec(): RA = LST - HA, with apparent sidereal time.raDec2AltAz()uses LST directly and expects epoch-of-date RA/Dec.Orbit.calcOrbit()performs the same derotation and conversion without precession, and labelsthe result "apparent ground-fixed, epoch of date".
SampleTrajectoryPositionsitself prints "The horizontal coordinates are apparent above a fixedground, not topocentric in J2000!".
Why at the final point
The local horizon rotates along the path by ~0.009 deg per km of ground track, so
traj.orbit.azimuth/elevation_apparent_norot, evaluated at the reference point, is not valid atthe end. The difference is tenths of a degree for final points tens of km downrange, and several
degrees for shallow trajectories propagated hundreds of km.
Method
DynamicMassFituses the same derotation asSampleTrajectoryPositionsandOrbit.calcOrbit(),evaluated at the final point and time. It uses the Earth rotation velocity and east direction at
that point, the radiant from
traj.radiant_eci_mini, and the final JD and lat/lon for theconversion, without precession.
The derotation uses the average speed
v_avg_norot, asSampleTrajectoryPositionsdoes. Drag actson the air-relative velocity, so the air-relative direction is conserved and the path is straight in
the Earth-fixed frame. Recovering that direction from a straight-line ECI fit requires an average
speed over the fitted segment. This was tested on synthetic fireballs with a known air-relative
direction, fitted as the solver does. Average speeds recover the direction to <= 0.02 deg for a slow,
strongly decelerating fireball, and to <= 0.005 deg for medium and fast ones.
v_initor theterminal ~3 km/s speed would be off by ~1 deg in the slow case.
Validation
The reference was computed independently of wmpl with pyerfa (IAU SOFA), from kinematics only: the
inertial velocity rotated to Earth-fixed with the Greenwich apparent sidereal time, minus w x r,
projected onto the local geodetic ENU basis, with no RA/Dec conversions and no derotation formula.
It agrees with astropy to 0.4" and with a finite difference of Earth-fixed positions to 3e-6 deg.
On several real trajectory solutions (slow and fast fireballs, both hemispheres):
0.14-0.34 deg (16-37 km) away, so the trajectory frame is epoch of date.
traj.orbit.azimuth/elevation_apparent_norotmatches the reference to < 5e-6 deg.DynamicMassFit: the returned final azimuth/elevation matches the reference to <= 1e-4 deg.SampleTrajectoryPositions: before this PR it is off by 0.13-0.54 deg in azimuth or elevation atevery sample. That output matches reading the radiant as J2000 (to ~0.01 deg, i.e. nutation),
not as epoch of date. After this PR it matches the reference to <= 1e-4 deg.
The removed error has the signature of precession. For a fixed ground-fixed direction evaluated
between 1980 and 2050, it vanishes at J2000, changes sign across it, and grows at 48.65"/yr and
35.45"/yr for two test directions. General precession predicts 48.64"/yr and 35.51"/yr (50.29"/yr
times the sine of the angle to the ecliptic pole). Without the step, the error stays <= 4e-5 deg over
the whole range.
TrajectoryKMLonly uses the sampled lat/lon/height, so it is unaffected. Both CLIs run end to end.🤖 Generated with [Claude Code]