space::irbem

Whole strips of igrf_strip_points go through detail::igrf_evaluate_strip; the remainder (at most a strip minus one) runs the scalar kernel, which computes the identical bits, so a caller can neither observe nor need to care where the boundary fell. This is what field_batch's host lane calls; it is public so a caller who KNOWS the batch is host-bound can skip the device-crossover reasoning entirely.

Measured by bench/irbem_bench.cpp: 100 ns/point against the scalar loop's 362 at 2¹⁶ points, degree 13, pinned to a P-core — 3.6×; see the file brief for the methodology and the variants that measured worse.

Computed rather than tabulated: Frame's values run GDZ=0, GEO..MAG=1..6, SPH=7, RLL=8, HEE..HEEQ=9..11, so removing the three angular frames is two contiguous shifts. A switch would need a default arm that no correct call can reach, and unreachable code is code that cannot be tested.

The only function permitted to build one.

The only function permitted to build one, so every context that exists has been validated.

The checks run in a fixed order — epoch, tilt, drivers, rotations — so a caller with several problems is told about them in that sequence rather than in an order that depends on the build.

Every driver is checked for finiteness and none is range-checked. That is deliberate: a caller legitimately leaves the slots its chosen external model ignores filled with a placeholder, so refusing a negative solar-wind speed here would reject a perfectly valid internal-field-only evaluation. Whether a driver is inside the range its model was fitted over is that model's question, asked where the driver is actually consumed.

The exact closed form, e.g. Torge, Geodesy (3rd ed.) §4.1, or Hofmann-Wellenhof et al., GPS: Theory and Practice, eq. (10.1): x = (N+h)cosφ cosλ, y = (N+h)cosφ sinλ, z = (N(1-e²)+h) sinφ, with the prime-vertical radius N = a/W. The (1-e²) on z is the ellipsoid: the surface normal at latitude φ does not pass through the centre, so the vertical drops short of the equatorial plane by exactly that factor.

Latitude comes from detail::geodetic_latitude (Bowring 1976, iterated). Altitude then uses the projection form h = p cosφ + z sinφ - aW rather than Bowring's h = p/cosφ - N. The two are algebraically identical — substituting the forward transform gives p cosφ + z sinφ = N(1 - e² sin²φ) + h = aW + h — but this one has no cosine in a denominator, so it is exact on the polar axis and loses no digits approaching it.

The textbook spherical-polar relation with latitude in place of colatitude: x = r cosφ cosλ, y = r cosφ sinλ, z = r sinφ. No ellipsoid is involved — SPH is GEO re-expressed, nothing more — so this is exact and its inverse car_to_sph is a true inverse.

Latitude is taken as 90° - θ with the colatitude θ = atan2(sqrt(x²+y²), z), and longitude as atan2(y, x). Both are atan2 rather than acos/atan for the usual reason — full range, no cancellation as the argument approaches ±1, and no division by a vanishing radius.

RLL is the awkward frame: a geocentric radius carried alongside a geodetic latitude (IRBEM's sysaxes table calls this out — "the latitude is still geodetic latitude and is therefore not interchangeable with SPH"). So the job is not a spherical-to-geodetic conversion; it is to solve for the altitude h that puts the point at geocentric radius r along the ellipsoid normal at a given geodetic latitude. Squaring the forward transform makes that a quadratic in h: r² = (N + h)²cos²φ + (N(1-e²) + h)²sin²φ = h² + 2Bh + C, B = N(1 - e² sin²φ) = aW, C = N²(cos²φ + (1-e²)² sin²φ), h = -B + sqrt(B² - C + r²).

and the discriminant simplifies exactly, since B² - C = -N² e⁴ sin²φ cos²φ: h = -B + sqrt(r² - N² e⁴ sin²φ cos²φ).

The positive root is the one taken: the other places the point on the far side of the axis.

λ_geo and r₀ come from the Earth–Moon barycentre's mean elements (F&H Table 4) fed through the equation-of-centre approximation of F&H eqn. 36; Ω from eqn. 14; and θ from eqn. 17, evaluated as a two-argument arctangent of (cos i · sin(λ−Ω), cos(λ−Ω)) so that the quadrant follows from the geometry instead of from a sign convention. F&H's prose ("the quadrant of θ is opposite that of λ−Ω") does not describe this construction; the two-argument form is what reproduces F&H's own worked example (their Table 8 implies θ = 259.899186°, against the 259.89919° they print), so the prose is read as a slip and the arithmetic is trusted.

Both frames are centred on the Sun, so this is a pure rotation for a position and for a field alike — T(HAE_D→HEE_D) = E(0, 0, λ_geo), F&H §3.2.2.

T(HAE_D→HEEQ) = E(Ω, i, θ) (F&H §3.2.2): swing +X onto the ascending node of the solar equator, tip by the 7.25° inclination, then rotate within the equator until +X lies under the Earth. Both frames are heliocentric, so again there is no translation.

The series (Aoki et al., Astron. Astrophys. 105, 359, 1982) is conventionally written for 0ʰ UT1 and then advanced by the elapsed day; the equivalent whole-instant form used here is GMST [s] = 67310.54841 + (876600ʰ·3600 + 8640184.812866)·T + 0.093104·T² − 6.2e-6·T³

with T in Julian centuries of UT1 from J2000.0. The leading 67310.54841 s is the familiar 24110.54841 s at 0ʰ plus the half day between midnight and the J2000.0 noon epoch.

The large linear coefficient is evaluated exactly rather than as written: 876600·3600 = 3155760000 is precisely 36525 · 86400, so that term is exactly 86400 · d for d days since J2000.0, and since the result is reduced modulo a day only its FRACTIONAL part matters. Dropping the whole days before multiplying keeps the sum from carrying an integer of order 10⁹ that is about to be discarded. Measured against an extended-precision evaluation of the same series over ±100 years, that is worth a factor of ~300: 1.5e-11° against 4.6e-9°, pinned by IrbemGmst.TheFractionalDayFoldingIsWorthAFactorOfHundreds.

Neither figure is anywhere near mattering physically — both are microarcseconds, and the coordinate-transform line of docs/ERROR_BUDGET.md is 1e-10 RELATIVE. It is done because it is free, not because the naive form would break anything.

(2): θ = 100.461 + 36000.770·T0 + 15.04107·H degrees, with T0 and H as in detail::hapgood_epoch.

This is a linear truncation of the IAU series and drops its quadratic term, so it drifts slowly away from gmst_iau1982_degrees — about 2 arcseconds near J2000, growing to a few tens of arcseconds a century out. It exists so a result can be reproduced against codes that use it.

(5), itself the Almanac for Computers low-precision series: M = 357.528 + 35999.050·T0 + 0.04107·H Λ = 280.460 + 36000.772·T0 + 0.04107·H λ☉ = Λ + (1.915 − 0.0048·T0)·sin M + 0.020·sin 2M ε = 23.439 − 0.013·T0

all in degrees, with T0 and H as in detail::hapgood_epoch. The sin M term is the equation of centre truncated at second order in the eccentricity; the series is accurate to about 0.01° in λ☉ over 1950–2050, which propagates to 0.01° in the GSE/GSM axes — far inside the budget for a frame whose defining direction is the Sun.

The Sun sits on the ecliptic by definition, so its ecliptic latitude is zero and the standard ecliptic-to-equatorial rotation by the obliquity reduces to (cos λ☉, cos ε · sin λ☉, sin ε · sin λ☉).

The transposed-direction case of rotation_matrix; see it for the full description.

The stored-direction case of rotation_matrix; see it for the full description.

Three constrained overloads share this name and this documentation, because they are three cases of one operation and a caller never chooses between them:

• To equal to From is the identity. It exists so that generic code parameterized on a target frame compiles for every frame including its own, instead of needing a special case at every call site.

• A pair Rotations stores in that direction is a read of nine already-computed numbers.

• The remaining pairs are stored the other way round, and are returned as the TRANSPOSE. A frame rotation is orthogonal — no scale, no shear — so its inverse is exactly its transpose. That is why the reverse direction is not stored: the transpose is a re-indexed read that cannot disagree with the forward matrix, whereas a second stored copy could, and a general inversion would cost a determinant and nine divisions to reproduce numbers already held exactly.

Any other pair fails RotationAvailable and does not compile.

The whole cost of a transform, once rotations exists: one 3×3 matrix–vector product, no branches, no allocation, and the frame checked entirely at compile time. Call it as transform<Frame::GSM>(p, rotations) — the source frame and whether this is a position or a field are deduced from the argument, so a GEO-to-GSM transform simply cannot be handed a GSM input.

Supported pairs are the eight Rotations stores and their eight inverses; anything else is a compile error, and must be composed through GEO deliberately.

The batch form is not a convenience wrapper. It is the only shape in which the problem is parallel: one L* point is Nder independent root-finds of a few traces each, and a device dispatch does not pay for Nder = 25 lines (gpu/dispatch.hpp measures the crossover at ~512). Hand over ntime points at once and every stage becomes ntime × Nder wide, which is where the measured 48.9× on the trace kernel actually lands. A loop calling make_lstar per point cannot be accelerated — the same fact lstar.hpp states about trace_invariant.

Each point is treated as a LOCALLY MIRRORING particle at pitch_angles_deg, matching IRBEM's make_lstar (90°) and make_lstar_shell_splitting (arbitrary). The shell that particle drifts on is the set of field lines carrying its (B_m, I), and Φ is the flux through the polar cap those lines' footpoints enclose.

The reference lane, and the one to reach for when there is a single point to compute — but NOT the fast one. Nder = 25 root-finds is a batch of 25, an order of magnitude below the ~512-line crossover gpu/dispatch.hpp measures, so a single L* runs almost entirely on the host however much hardware is present. That is a property of the problem, not of this implementation: the parallelism in L* lives across POINTS, and make_lstar_batch is where it is taken.

Four groups: {0, 0+}, {1-, 1, 1+, 2-}, {2, 2+, 3-} and Kp >= 3. In the Kp x 10 scaling the attainable values are 0, 3, 7, 10, 13, 17, 20, 23, 27, 30, ..., so the group edges are 4 (between 0+ = 3 and 1- = 7), 20 (2- = 17 from 2 = 20) and 30 (3- = 27 from 3 = 30). The edges are inclusive on the upper side — 20 is bin 3 and 30 is bin 4 — which the oracle switch points in tools/oracle/mead_diff.cpp confirm. A negative or NaN Kp (which check_validity reports separately) yields bin 1 and anything at or above 30 yields bin 4.

The float instantiation is how the host lane is made to run the same arithmetic as the device kernel: rounding the coefficients on the way in, rather than evaluating in double and rounding the answer, is what makes a host-vs-device disagreement attributable to the DEVICE.

This is the whole model. The order is: GSM to SM by the tilt (a rotation about y, the same one t89_components makes and for the same reason — the paper's frame has the dipole axis as z); the SM position to its aberrated twin (x_m, y_m) by a rotation about z through 4 degrees; the three polynomials of the file brief; and the SM components back to GSM. The aberration rotates the POSITION only — see the brief for the measurement that says so.

The tilt arrives three times over: as its sine and cosine for the frame rotation, and as the angle in degrees for the polynomials. All three are properties of the epoch, not of the point, so they are paid once per timestamp by the caller rather than recomputed here — there is no trigonometry in this function at all. The three MUST describe the same angle; the entry points below guarantee it by deriving all three from one tilt_rad.

The value is always returned, including when the status is Status::OutOfValidityRange — status.hpp's standing rule. What is refused outright is input that has no meaning: a non-finite coordinate, tilt or Kp, a point inside the Earth, or a tilt beyond a right angle. The last is not a singularity as it is in t89_field (there is no tan(psi) here, and psi = 90 deg itself has a value) but a definition: the geodipole tilt is the angle between the dipole axis and the GSM z axis and lies in [-90, 90] degrees, the same bound make_field_context enforces. Within it, and with a finite radius, the three quadratics cannot overflow — every squared coordinate is bounded by the finite r^2 and every coefficient by 30 — so no check on the OUTPUT is needed, and none is made.

What IS checked against the published envelope, both through status.hpp so the rules live in one place: Kp against 0 <= Kp <= 9 (in IRBEM's Kp x 10 scaling) and the position against the paper's r <= 17 R_E. The Kp check happens before the binning, so a caller who passes Kp = 12 gets the Kp >= 3 set AND is told the model was never fitted there.

The host twin of irbem_mead_f32: the same expressions, in the same order, in the same precision, from coefficients rounded to float FIRST — what makes a disagreement between the two lanes attributable to the device. It is also the lane a machine with no GPU actually runs.

The layout is the kernel's ABI and is stated in exactly two places — here and the comment above irbem_mead_f32 in irbem.slang: [0] sin psi, [1] cos psi, [2] psi in degrees, [3..9] a_1..a_7, [10..12] b_1..b_3, [13..19] c_1..c_7. A test evaluates both lanes on the same points, which is what actually keeps the two statements in step.

The shape mirrors t89_field_batch exactly — device above the registry's measured crossover, the fp64 host loop otherwise, one folded Status for the batch, and a returned value that says which lane actually served the call. What differs is the arithmetic intensity: this model is ~50 flops for 24 bytes in and 12 out, about 1.4 flops/byte, which is within a factor of three of the streaming dipole kernel that LOSES on this seam and an order of magnitude below T89's ~11. The irbem_mead_f32 row of gpu/dispatch.hpp carries the measurement and the verdict.

The batch reports the same caveats the scalar entry point does, folded over the whole batch. If any point is beyond the published r <= 17 R_E the batch says Status::OutOfValidityRange and is still computed in full; if any point is inside the Earth or not finite it says Status::DomainError and every output is zeroed.

P = m_p n V^2, the ram pressure of a proton wind. This is the ONLY combination of the two solar-wind drivers the model reads — a fact about the published idea (Chapman-Ferraro pressure balance knows nothing but P) that the oracle turns out to share bit for bit.

Pressure balance at the subsolar magnetopause, P = B_mp^2 / 2 mu_0 with B_mp ∝ R_mp^-3, puts the boundary at R_mp ∝ P^(-1/6) (Chapman & Ferraro 1931; the textbook statement is Kivelson & Russell, Introduction to Space Physics, 1995, §6.2). A self-similar compression of the whole boundary-current system by the factor s then scales its field as s^3 B(s r) (Mead, JGR 69, 1181, 1964, whose boundary-field coefficients all carry R_mp^-3). This function is the s; the evaluator applies the transform.

O'Brien & McPherron (2000), with b = 7.26 nT/sqrt(nPa) and c = 11 nT: the index with the magnetopause currents' contribution removed, leaving the part the ring current is responsible for. This is the "modified Dst" the model's description names as its ring-current driver.

The quiet field already contains a quiet ring current (T89's C_5 term), so what the Dst driver adds is the DIFFERENCE from quiet — which is what makes the model reduce to the quiet field exactly at the reference conditions rather than double-counting a ring. The offset c cancels in the difference; only Dst and the pressure correction survive.

The ring is the finite-moment disc potential of Tsyganenko (1989) eq. (9), with radius a and half-thickness D. On its axis at the origin its field is B_z(0) = 2 C / (a + D)^3 (set rho = z = 0 in eq. 16-17's B_z = C (2 u^2 - rho^2) / S^5 with u = a + D, S = u) — z being the DIPOLE axis, since the ring lies in the dipole equator; in GSM the central field points along (sin psi, 0, cos psi). The Dessler-Parker-Sckopke relation (Dessler & Parker, JGR 64, 2239, 1959; Sckopke, JGR 71, 3125, 1966) says the symmetric ring's field at the centre IS the pressure-corrected Dst, so the amplitude that makes B_z(0) = ΔDst* is C = ΔDst* (a + D)^3 / 2.

Computing in double and rounding the RESULT is what keeps the float block the nearest float to the true parameters rather than the product of float intermediates: s is a sixth root and C a cube, and both amplify a rounding in their argument.

A flat, symmetric A^(3) disc in the dipole equator: Tsyganenko (1989) eqs. (16)-(17) with the sheet unwarped (z_s = 0) and its thickness constant, so every dz_s/dx, dz_s/dy and dD/dx term of the published form is identically zero and what is left is q = 3 u / (xi S^5), B_x = q x z, B_y = q y z, B_z = (2 u^2 - rho^2) / S^5,

with xi = sqrt(z^2 + D^2), u = a + xi, S = sqrt(rho^2 + u^2), all in SM, and the result rotated back to GSM. Split out of opd_components so that the unit shape — the model's ∂B/∂Dst — can be tested and normalised on its own.

The one line of the file brief: s^3 B_q(s r) + C R(r). The quiet field is T89's evaluator on the SCALED position — a point at r in a magnetosphere compressed by s sees what the quiet magnetosphere has at s r — with the s^3 that keeps the transform a similarity of the currents; the ring is opd_ring_unit at the UNSCALED position, because the Dst-driven increment is a property of the inner magnetosphere and not of the boundary. Both are solenoidal, so their sum is.

The value is always returned, including under Status::OutOfValidityRange — status.hpp's standing rule. The oracle, by contrast, returns baddata strictly outside the same envelope (measured: n = 50.0001, V = 500.001, Dst = 20.001 and -100.001, r = 60.05 all refuse; the bounds themselves are accepted), and a caller porting from it should expect a number with a caveat where it got a sentinel.

Refused outright, as Status::DomainError with a zero field, is arithmetic with no answer: a non-finite input; a tilt of |psi| >= pi/2, at which the quiet field's tan(psi) does not exist; a radius inside the Earth; a density or speed that is not positive, for which there is no pressure and no magnetopause to compress (s would be zero or complex); and a non-finite answer, which the quiet field's exp(x / dx) produces far enough sunward that only a unit-confusion bug gets there, and which must not be handed to a tracer as a NaN.

The host twin of irbem_opd_f32: the same expressions, in the same order, in the same precision, from a parameter block rounded to float FIRST. What makes a host-vs-device disagreement attributable to the device.

The first thirty floats are EXACTLY t89_param_block for the quiet bin, because the kernel hands them to the shared t89_eval untouched; the last two are the scale and the ring amplitude. Stated here and above irbem_opd_f32 in irbem.slang, and kept in step by a test that runs both lanes on the same points.

The shape of t89_field_batch exactly: one status for N points, the worst of them; an out-of-validity batch computed in full and a domain-error batch zeroed in full; the returned value true exactly when the device serviced the call. The kernel is T89's plus a ring, at the same ~11 flops/byte, and runs the SAME shared t89_eval — so the device wins by the same margin and the crossover is the T89 row's (see gpu/dispatch.hpp).

TT = {1, t, t^2, t^3} with t the tilt in DEGREES; an even term takes pair[0] + pair[1] t^2, an odd term pair[0] t + pair[1] t^3. Done in double whatever T is, then rounded: the rounding on the way in is what makes the float lane the device's twin.

This is the report's subroutine, block for block: the radial rules first (zero inside 2 R_E, a linear-in-r^2 taper to 2.5 R_E, the zero template beyond 15 R_E), then the envelope E = exp(-0.06 r^2), then the series in the report's enumeration. That enumeration is a fixed pattern that reads as three nested loops with two truncation rules — the trip count depends on nothing but the loop counters, so the work is the same at every point and a device lane does not diverge:

• x power p = 0..4 outermost; y^2 power q = 0..2 next, cut off when p + 2q > 5;

• z power s innermost, from 0, cut off when s > 4 OR the running index p + 2q + s passes 5 — with B_y getting one term FEWER than B_x/B_z in each block, because the overall factor of y costs it one order.

Written with the report's own bookkeeping counters (ijk, k) rather than the closed-form rules above, so that the 32 + 22 + 32 terms come out in exactly the order the tables are laid out in — the counts are asserted, not assumed.

The rotation in and out is the one every Tsyganenko-family evaluator makes: about y by the tilt, GSM to SM, because the report's series is defined in SM (its z IS the dipole axis).

Three checks, through status.hpp so the rules live in one place: finiteness of the inputs, the tilt against |psi| <= pi/2 (the model has no tan(psi) and is total at a right angle, so the bound is the one FieldContext already guarantees and nothing stricter), and the position against the published r_GEO <= 15 R_E. There is NO driver check because the model reads no driver; check_validity is still consulted so that the envelope table, not this file, is the authority on that.

The value is always returned, status.hpp's standing rule — and beyond 15 R_E that value is a ZERO external field, because that is what the published model says there: the report's template "sets the field to zero" where "the power series diverges". A caller who extrapolates this model gets the published extrapolation, which is no external field at all, and is told so. Below 2 R_E the field is likewise the published zero, and the status is Status::Ok because that region is inside the envelope: the report chose zero there on purpose.

The host twin of irbem_opq_f32: the same expressions, in the same order, in the same precision, from coefficients folded in double and rounded to float FIRST. That is what makes a disagreement between the two lanes attributable to the device — a contraction, a driver's exp — rather than to the arithmetic having been written differently on the two sides.

The layout is the kernel's ABI and is stated in exactly two places — here and the comment above irbem_opq_f32 in irbem.slang. A test evaluates both lanes on the same points, which is what actually keeps the two statements in step.

This is the routine to call for more than a handful of points. opq_field is the reference lane: it is what the batch is verified against, and what runs when there is no device or the batch is too small to pay for one.

The arithmetic is T89's regime: 86 fused polynomial terms and one exp for 24 bytes in and 12 out, ~15 flops/byte, no data-dependent branch, and a 174-float parameter block read identically by every lane in a workgroup — see the irbem_opq_f32 row of gpu/dispatch.hpp for the measurement and the crossover derived from it.

The batch reports the same caveats the scalar entry point does, folded over the batch. If any point is beyond 15 R_E the batch says Status::OutOfValidityRange and is still computed in full (those points come back as the published zero); if any point is inside the Earth or not finite it says Status::DomainError and every output is zeroed.

The shape mirrors the internal field's trace_invariant_batch, and the return value says which lane served the call for the same reason: a silent fallback is what makes a performance claim worthless. This overload is the host lane only. The device tracer that composes an external model (irbem_trace_total_f32) is written for T89's parameter block; composing OP-77 on the device would need its own tracer kernel, and until one exists the honest answer is false here every time, never a quiet substitution of the internal field. Single traces and drift shells still run — in fp64, on the host, at the reference lane's cost.

(2)'s inner sum.

Done once per epoch rather than per point, because the drivers are a property of the epoch — the same economy HotState makes with the tilt's trigonometry. The result is what the evaluator, the host batch lane and the device kernel all consume, so the three lanes cannot disagree about what the drivers meant.

Computed in double and rounded to T at the end — the fp32 lane wants the amplitudes the kernel receives, and the kernel receives them as floats made from the best available doubles.

Every row is the printed cylindrical entry with cos phi = x / rho, sin phi = y / rho substituted and rho cancelled, which turns each into a polynomial: B_x = b_rho cos phi - b_phi sin phi, B_y = b_rho sin phi + b_phi cos phi. For a row of the form (f cos phi, g sin phi, h cos phi) that is B_x = (f x^2 - g y^2) / rho^2, B_y = (f + g) x y / rho^2, B_z = h x / rho, and in every row the numerators carry the rho^2 or rho needed to cancel. The tilt rows carry their own sin psi, as printed.

The four irrational leading factors are the Schmidt-normalized Legendre coefficients: sqrt(3) (row 7, n = 2), sqrt(10) (row 8, n = 4), sqrt(3/2) (row 16, n = 3) and sqrt(15) / 8 (row 17, n = 5 — see the file brief for why the printed 1 / (8 sqrt(15)) is a misprint of that).

The whole model in one straight line: rotate GSM to SM about y by the tilt (SM's z IS the dipole axis, which is the frame the harmonics are defined in), divide by 10 R_E, evaluate the 17 harmonics, weight them by the amplitudes, rotate back. No loop with a data-dependent trip count, no branch, no transcendental: the tilt arrives as its sine and cosine.

Separate from status.hpp's check_position because the box is stated in SM, which that function does not know about, and because IRBEM's table — which that function reads — publishes no limit for this model at all. The bounds are closed: a point exactly on the boundary is inside.

Two layers, in status.hpp's order. First check_validity for this model's row, which catches a non-finite driver among the four this model reads (and, the IRBEM table being unbounded for kext = 8, nothing else). Then the paper's own caveat, which the table does not carry: Dst < -200 nT is declared invalid. The interval is closed, so -200 itself is inside.

The value is always returned, including when the status is Status::OutOfValidityRange — status.hpp's standing rule. What is refused outright is arithmetic that has no answer: a non-finite input, a point inside the Earth, a tilt that is not an angle to an axis, or an extrapolation so far that the fourth-order polynomial overflows (a unit-confusion bug's signature — kilometres where Earth radii were meant — and the one way this evaluator can produce a non-finite number from finite input). Unlike t89_field there is no refusal at |psi| = pi/2: the model carries sin psi and the rotation, never tan psi, so a right-angle tilt is merely unphysical, not undefined.

What is checked against the paper: the position against the fitted box (om97_check_fitted_region) and the drivers against its Dst floor (om97_check_drivers). Both report, neither suppresses.

The host twin of irbem_om97_f32: the same expressions, in the same order, in the same precision, from the same float amplitudes. That is what makes a disagreement between the two lanes attributable to the device rather than to the arithmetic having been written differently on the two sides. It is also the lane a machine with no GPU actually runs for float batches.

The layout is the kernel's ABI and is stated in exactly two places — here and the comment above irbem_om97_f32 in irbem.slang. A test evaluates both lanes on the same points, which is what actually keeps the two statements in step.

This is the routine to call for more than a handful of points. om97_field is the reference lane: it is what the batch is verified against, and what runs when there is no device or the batch is too small to pay for one.

The amplitudes are computed ONCE for the batch, which is the whole economy of the epoch model: per point the kernel does ~250 flops of polynomial over 24 bytes in and 12 out, a little under T89's ~11 flops/byte and well above the streaming dipole's 0.5. Its measured crossover is in the irbem_om97_f32 row of gpu/dispatch.hpp.

The batch reports the same caveats the scalar entry point does, folded over the whole batch. A DOMAIN-ERROR batch is zeroed in full and never reaches the device; an out-of-validity batch is computed in full.

False only for bin 7, which repeats bin 6's published Kp >= 5- column. Exposed rather than left to a comment because a caller running a storm-time study deserves to know that the two highest bins are not distinguished by this implementation.

The float instantiation is not a convenience: it is how the host lane is made to run the same arithmetic as the device kernel, which receives these as float. Rounding the coefficients on the way in — rather than evaluating in double and rounding the answer — is what makes a host-vs-device disagreement attributable to the DEVICE.

This is the whole model, and it is written so that each block is one equation of the paper. The order is the paper's order, and the two coordinate systems are kept apart deliberately (see the file brief): the tail and ring current are evaluated in SM and rotated to GSM; the closure and Chapman-Ferraro terms are evaluated in GSM directly.

On the derivatives. Every d/dx and d/dy below is the analytic derivative of the expression above it, and none is optional: the field is the curl of a vector potential whose x and y dependence runs through W, through the sheet surface z_s and through the thickness D, so dropping any of them does not merely approximate the answer, it breaks div B = 0 — which is why a finite-difference divergence test catches exactly this class of mistake and is the main test this file carries.

With P = A_phi / rho the potential-over-radius, the curl of A = P * (-y, x, 0) is B_x = -x dP/dz, B_y = -y dP/dz, B_z = 2P + x dP/dx + y dP/dy — which is where the paper's eqs. (14)-(17) and (19) come from, and the form in which they are evaluated here.

The value is always returned, including when the status is Status::OutOfValidityRange — that is status.hpp's standing rule and the whole reason a Result exists: extrapolating an empirical fit is a decision only the caller can make. What is refused outright is arithmetic that has no answer: a non-finite input, or a tilt of exactly ±90 degrees, at which eq. (11)'s tan(psi) does not exist.

Two things are checked against the published envelope, both through status.hpp so the rules live in one place: Kp against 0 <= Kp <= 9 (in IRBEM's Kp x 10 scaling), and the position against T89's r_GEO <= 70 R_E. The Kp check happens before the binning, so a caller who passes Kp = 12 gets the most disturbed published set AND is told the model was never fitted there.

There is one more refusal and it is about the OUTPUT, not the input. Eq. (20) carries exp(x / dx) with dx ~ 20 R_E, so a caller who extrapolates far enough sunward — past about 1e4 R_E, which is nothing a trace can reach but is exactly what a unit-confusion bug produces when kilometres arrive where Earth radii were meant — overflows it. The field is then infinite, and inf * 0 in the B_y assembly makes a NaN out of a component that is EXACTLY zero everywhere the model means anything (IrbemT89.OnTheNoonMidnightMeridianTheFieldStaysInThatPlane). A NaN that escapes here surfaces a hundred RK4 steps later with nothing left pointing at its cause, so a non-finite answer is refused the same way a non-finite input is.

The host twin of irbem_t89_f32: the same expressions, in the same order, in the same precision, evaluated from coefficients rounded to float FIRST. That is what makes a disagreement between the two lanes attributable to the device — a contraction, a driver's exp — rather than to the arithmetic having been written differently on the two sides. It is also the lane a machine with no GPU actually runs.

The layout is the kernel's ABI and is stated in exactly two places — here and the comment above irbem_t89_f32 in irbem.slang. A test evaluates both lanes on the same points, which is what actually keeps the two statements in step.

This is the routine to call for more than a handful of points. t89_field is the reference lane: it is what the batch is verified against, and what runs when there is no device or the batch is too small to pay for one.

T89 sits in the regime where the device wins. It is ~400 flops for 24 bytes in and 12 out — ~11 flops/byte, an order of magnitude above the streaming dipole kernel that LOSES 0.69x on this same seam and within a factor of two of IGRF, which wins 8.96x. Nothing is stored per point beyond the answer, there is no branch on data, and the 30-float parameter block is read identically by every lane in a workgroup.

The batch reports the same caveats the scalar entry point does, folded over the whole batch. One status for N points can only be the worst of them, so that is what it is: if any point is beyond T89's published r_GEO <= 70 R_E, the batch says Status::OutOfValidityRange, and if any point is inside the Earth or not finite it says Status::DomainError. An out-of-validity batch is still computed in full — status.hpp's standing rule — and a DOMAIN-ERROR batch is zeroed in full, which is exactly what t89_field does with its one point. Checking positions costs the batch two comparisons per point and two envelope lookups per CALL: the radii are folded as squares, so there is no per-point sqrt and the device lane's measured throughput is unchanged.

Proportional to the radius on purpose: truncation error scales as h/r and cancellation as r/h, so their crossing point moves with r and a single absolute step is optimal at exactly one altitude. IRBEM takes an absolute dX for the whole batch, which is why this is a default rather than the only option — bderivs and bderivs_batch accept an explicit step, and a differential comparison must use one so both implementations are differenced identically.

Four field evaluations: the base point, then x + h·ê_x, x + h·ê_y, x + h·ê_z. Forward, not central, because that is what the reference does — established black-box, and exactly (the measured relative difference is 0.0, not "small"). A central difference would cost seven evaluations and buy two orders of magnitude; it is offered by nobody here because agreeing with IRBEM is the contract, and the accuracy that matters is quantified in the file brief instead of silently improved.

The gradient of the magnitude and the Jacobian of the vector are differenced from the SAME four evaluations. That is not merely an economy: it is what makes grad_b_mag and diff_b mutually consistent to the same order, which is what grad_curv_curl's grad_par and curvature silently assume when they compute the same parallel derivative two different ways.

No field model, no evaluations, ~50 flops. This is the routine the file brief measures at 19 ns/point on the host and rules off the device permanently.

The one subtlety is curvature. The identity (B̂·∇)B̂ = [ − (·B̂)B̂] with  = (B̂·∇)B/|B| holds because B̂·B̂ = 1 forces the derivative of B̂ to be perpendicular to it. Writing it that way — projecting out Â's own parallel part rather than subtracting grad_par·B̂/|B| — is what the reference does, and the difference is not cosmetic: the two agree only when grad_b_mag and diff_b are mutually consistent, which finite differences make them only to first order. The form used here needs grad_b_mag not at all, so curvature and r_curv are exactly perpendicular to B̂ by construction rather than approximately.

The criterion is the sign of d|B|/ds along +B̂: the field falls to a minimum at the magnetic equator and rises toward both feet, and B points from the southern foot over the equator to the northern one, so a rising field in the direction B points means the equator is behind you and you are north of it. That is one signed scalar and it is exactly GradCurvCurl::grad_par, which is why this routine and the derivative routines are in the same file.

Computed with a central difference along B̂ rather than the forward differences of bderivs — three evaluations instead of four, and symmetric, which matters here specifically: near the equator the true derivative passes through zero, and a one-sided step can overshoot the minimum and report the wrong side of a boundary the answer is a discrete function of. A derivative that is merely inaccurate is fine; a sign that is wrong is a different hemisphere.

The reference is a bare loop over points, and it is the most trivially parallel routine in the whole library: one point in, one vector out, no state carried between iterations. It is also the least arithmetically intense thing here worth offloading — ~1 900 flops for 24 bytes in and 24 out, about 20 flops/byte — which is exactly why the crossover is consulted rather than assumed. Measured on an RTX 3070 Ti against this file's own fp64 host lane, -O3 -march=native -ffp-contract=off, best of five, transfers included:

— 14.2 ns/point on the device at 2¹⁴ against the host's 309. Below field_batch_crossover the ~115 µs dispatch floor alone exceeds the whole computation and the host wins; a per-point loop calling Igrf::evaluate can never reach the device at all, which is the reason this entry point takes the whole batch. The curve turns over slightly at 2²⁰ because the routine is staging-bound by then, not kernel-bound: the same dispatch measured in isolation costs 5.9 ns/point at 2²¹, so the remaining ~12 ns is this file converting fp64 positions down to fp32 and fp32 fields back up, which is inherent to the typed API and not to the device.

The device lane returns fp32. Measured maximum relative deviation against the fp64 host lane: 7.3 × 10⁻⁷ over 2 000 points and 1.1 × 10⁻⁶ over 2²⁰ — so at the larger sample it just exceeds the 1 × 10⁻⁶ Bgeo budget of docs/ERROR_BUDGET.md §4. That is stated rather than rounded away: it is the cost of summing 105 harmonic terms in fp32, it grows like the tail of a distribution as the sample grows, and the fix if a caller needs the budget honoured at every point is the host lane, not a different kernel. Nothing here accumulates across points — one point, one thread, no reduction — so the budget's reduction concern still bites only the integrals in lstar.hpp.

One dispatch, not four. Each point needs four field evaluations — the base and three one-sided neighbours — and the obvious implementation issues four batched dispatches of N points each. That pays the ~115 µs dispatch floor four times for work that has no dependency between the four groups whatsoever. Building the 4N points up front and dispatching once pays it once, and at the same time quadruples the occupancy of the launch — which is what pulls the crossover down to bderivs_batch_crossover, a quarter of field_batch_crossover. Measured: 0.64× at 64 points, 1.27× at 128, 2.49× at 256, 8.44× at 1 024, 26.1× at 2¹⁴ (49.9 ns/point against the host's 1.30 µs). Four dispatches instead of one would have moved that crossover to 512 points and left everything below it on the host.

The layout is [all N base points][all N +x][all N +y][all N +z] rather than four consecutive points per input point. Same dispatch either way, but this way each of the four groups is a contiguous run whose lanes read neighbouring coefficients in the same order, and the host-side differencing walks four unit-stride streams instead of one stride-4 one.

On the device the differenced values are fp32, and the file brief's tables are the consequence: the achievable accuracy is ~1.5 × 10⁻³ relative against ~1 × 10⁻⁷ on the host, and the step that achieves it is eight thousand times larger. That four-order gap is why step_re defaults to the lane's step rather than to a constant, and why a caller comparing the two lanes must pass an explicit step to both — at a matched dX = 10⁻² they agree to three digits, and at a matched dX = 10⁻⁵ the device is wrong by 25%. It is a real limitation, quantified here rather than hidden behind an average.

One detail that is worth more than it looks: the differencing divides by the step the device actually took, recovered from the fp32 position it was handed, not by the step that was intended. x + h is not representable in fp32, and at h = 4 × 10⁻⁴·r the rounding is a relative error on h of order ε₃₂·r/h ≈ 1.5 × 10⁻⁴ — comparable to everything else in this routine's budget, and free to remove.

A host loop, permanently, and the number that settles it. 136 bytes would have to cross the bus per point (16 fp32 in, 18 out) for 19 ns of host arithmetic — measured, 17.7 ns/point at 2¹⁰ and 19.9 at 2¹⁸, i.e. 50–57 Mpts/s, -O3 -march=native -ffp-contract=off. This seam's payload bandwidth is 7.15 GB/s, measured directly by dispatching a degree-1 IGRF kernel over 2²¹ points — 30 flops of kernel, so what is left is the round trip — which puts 136 B/point at ≥ 19.0 ns/point of transfer alone. The copy costs exactly what the computation costs, before the kernel runs and before the ~115 µs dispatch floor is paid. Arithmetic intensity is ~0.4 flops/byte, below even the dipole kernel's 0.5 — and gpu/dispatch.hpp already records that the dipole kernel LOSES 0.69× at every batch size it was measured at. There is no size at which this one wins: the ratio is fixed, not amortizable. So no kernel for it exists in gpu/irbem.slang and none should be written.

The routine is still worth having as a batch entry point: it keeps the loop in one place, it vectorizes, and it lets a caller hand over an ephemeris rather than write the loop themselves.

Three field evaluations per point — the base, and one step either way along B̂ — and the same one-dispatch discipline as bderivs_batch, with one difference that is forced by the physics: B̂ is not known until the base evaluation has come back, so this cannot be a single 3N dispatch. It is two: one of N points to get the directions, then one of 2N to get the neighbours. Two dispatches for three evaluations is the best a data dependency of this shape allows, and it is still a factor of N better than the per-point loop the name suggests.

fp32 costs almost nothing here, because only the SIGN of the difference survives. Measured over 2¹⁸ random points spanning r = 1.5 … 8.5, the device lane and the fp64 host lane disagree on 11 points in 262 144 — 4 × 10⁻⁵ — and every one of them is a point sitting within a step of the magnetic equator, where d|B|/ds passes through zero and the two hemispheres are genuinely adjacent. No precision decides that question; a point on the equator is in neither hemisphere. Measured speedup: 0.51× at 128 points, 0.97× at 256, 3.76× at 1 024, 25.4× at 2¹⁴.

Three stages, because I is defined between the two mirror points and a one-directional walk from an arbitrary start point covers only part of that path:

1. walk in the direction of DECREASING field to the magnetic equator, recording B_min;

2. from the equator, integrate outward to the mirror point;

3. from the equator, integrate outward the OTHER way, to the conjugate mirror point.

I is the sum. Doing only the first half and doubling it would be right for a centred dipole and wrong for every real field, which is asymmetric about the equator — and wrong by an amount that grows exactly where the models matter most.

The magnitude of the degree-1 field: √(g₁⁰² + g₁¹² + h₁¹²). Epoch-dependent — the Earth's moment has fallen roughly 6% over the last century — and L ∝ M^(1/3), so taking it from the model rather than from a constant is what keeps L_m right across epochs.

McIlwain defined L through a tabulated function; Hilton (1971) fitted a closed form accurate to about one part in 10⁴, which is what every implementation since has used. With X = I³·B_m/M and M the dipole moment, L³·B_m/M = 1 + 1.35047·X^(1/3) + 0.465376·X^(2/3) + 0.0475455·X

Higham §4.2 gives |E| ≤ (n−1)u / (1 − (n−1)u) for recursive summation of n terms. What is returned is the numerator (n−1)u, so it understates that bound by the denominator — which differs from one by (n−1)u itself, a part in ~10⁴ for the fp32 cases here and a part in ~10¹³ for fp64, both far under the factor-of-two the budget's own spelling already carries (see the note below). It is a worst case — every rounding error the same sign and maximal — and the realistic figure is the random walk of random_walk_estimate.

This is the √N·ε row of ERROR_BUDGET §3 (which, as in naive_sum_bound, states it in ε rather than u, so its ~4 × 10⁻⁶ at n = 10³ is twice what this returns). It is an estimate, rather than a bound — a summation whose errors correlate (a monotone integrand, which a field-line quadrature very much is) can exceed it, up to naive_sum_bound. Both are quoted in the budget because the honest margin lies between them: one to three orders above the discretization floor, not the two to three the planning estimate assumed.

Answering it is what decides where a reduction may run. At the XJ budget of 1 × 10⁻⁴ relative (ERROR_BUDGET §4) this returns ~1.7 × 10³ for float — the same order as the term count of a single field-line quadrature, i.e. the reduction would consume its whole budget — against ~9 × 10¹¹ for double, where the question stops mattering. That contrast is the accumulator invariant, in one number.

Spelled as a function rather than left to == Status::Ok at every call site so that "success" has one definition — if a future status ever has to be treated as benign, it is added here.

A caller that must satisfy several envelopes at once (drivers, position, coefficient files) wants one status out, and wants the first reason rather than the last, because the checks are ordered cheapest-and-most-fundamental first.

It is the enumerator's own value, so the kernel needs no table: the device-side code assigns the same small integers by hand and this function is the host-side statement of that contract.

Validated rather than cast, because the bytes come from a buffer a driver wrote: an uninitialised slot, a kernel that never ran, or a device-side bug produces an integer with no meaning, and static_casting it into a Status would be undefined behaviour dressed up as a result.

This is a lossy, one-way bridge: six statuses map onto one number and the reason is gone. There is deliberately no inverse — a caller who wants to know why keeps the Result.

Note that Status::OutOfValidityRange collapses to the sentinel here even though the value was computed, because IRBEM's convention has no way to say "here is the number, with a caveat". That loss is precisely the cost of the boundary, and precisely why nothing inside this library crosses it.

Exact equality, matching IRBEM's own .eq. baddata test. A NaN is not the sentinel — it compares unequal to everything, this function included — which is correct: a NaN that escaped a computation is a different failure from one the library deliberately reported.

This never suppresses a value. It is called beside the evaluation, not instead of it: the caller gets the number the functional form produces AND the knowledge that the number is an extrapolation. Deciding what to do about that — clamp, refuse, or proceed and annotate — is a scientific judgement this library will not make on the caller's behalf.

Only the drivers the model actually reads are examined, so a maginput vector full of fill values in the slots a model ignores is still Ok for that model. That matters in practice: a caller running the same vector through T89 and T96 has, at most, the union of what both need.

The two failures are ordered, most fundamental first: every driver the model reads is checked for finiteness before any is checked against its range, so a NaN anywhere in the used set reports Status::DomainError rather than whichever came first in the list.

The spatial envelope is stated two different ways by the two model families, and both are checked: the pre-2001 models publish a maximum geocentric radius (they stop having a tail), while T01 and later publish a minimum GSM x (their tail current sheet is what runs out). A model publishing neither accepts every finite point.

Only TS07D needs them: it is not a closed-form fit but a per-interval expansion whose radial basis function coefficients are downloaded per six-minute interval (IRBEM ships setup_ts07d_files.sh to fetch a Coeffs/ and a TAIL_PAR/ tree for exactly this reason). Without them the model has no parameters, which is a different failure from having bad ones — hence its own status rather than Status::DomainError.

T89 is not continuous in Kp. Tsyganenko (Planet. Space Sci. 37, 5, 1989) sorted the data into seven Kp intervals and fitted a separate coefficient set to each: {0, 0+}, {1-, 1, 1+}, {2-, 2, 2+}, {3-, 3, 3+}, {4-, 4, 4+}, {5-, 5, 5+}, and {>= 6-}. The returned index is the 1-based bin number, which is the iopt argument the published T89 interface takes.

The thresholds therefore fall at Kp x 10 = 5, 15, 25, 35, 45, 55 — between the values Kp can actually take, since the third-of-a-unit steps land on 0, 3, 7, 10, 13, 17, ... in this scaling. There is consequently no "exact boundary" case for a real Kp to sit on, and the bin edges are chosen midway rather than at a bin's own extreme so that the mapping is robust to a caller who has already rounded. The convention where it cannot matter is nonetheless decided and uniform: a threshold value belongs to the bin below it, so the bins are the half-open intervals (-inf, 5], (5, 15], ... (55, +inf) in this scaling. A negative Kp — which check_validity reports separately — still yields bin 1, and any Kp above 5.5 yields bin 7.

Walks downhill in |B| until the field stops falling, then fits a parabola in arc length through the three samples that bracket the minimum and takes a partial RK4 step to its vertex. The coarse sample is not the answer and returning it would be a silent O(ds) error in the position: measured on an L≈4 line at the default step, the coarse minimum sits 5.5 × 10⁻⁶ relative above the oracle's Bmin and 4.4 × 10⁻³ R_E from its position.

A start point that is ALREADY at the minimum is not a special case in the physics and is not one here: both neighbours are probed, and if both are higher the bracket is centred on the start.

Measured against the oracle over 84 start points × 4 epochs at matched IGRF: Bmin to 1.9 × 10⁻⁶ relative (budget 10⁻⁵) and the position to 4.6 × 10⁻⁵ R_E. The Bmin residual is the ORACLE's, not ours — our value at that resolution is within 10⁻⁸ of our own converged value, and the oracle sits the same 1.9 × 10⁻⁶ from it.

The mirror field is not searched for, it is KNOWN: the first adiabatic invariant makes it B_m = B_local / sin²α exactly. What the trace finds is where along the line that value is reached, by walking in the direction of INCREASING |B| — which is what puts the mirror point in the particle's own magnetic hemisphere, as the reference does.

α = 90° is returned without tracing at all: the particle mirrors where it is, so the answer is the input point verbatim and B_m = B_local. That is not an optimisation, it is the only answer that is exactly right — a trace would return a point one refinement away from where the particle demonstrably is.

A line that reaches PathTraceOptions::r0 before B_m is a particle in the LOSS CONE: it hits the atmosphere instead of mirroring. That is physics, so it is Status::OpenFieldLine and not an error — the reference answers baddata here, which cannot distinguish it from a bad input.

The termination is on GEODETIC altitude, which is why this converts inside the step loop rather than comparing radii and converting once at the end. The WGS-84 ellipsoid's polar semi-axis is 21.4 km shorter than its equatorial one and field-line feet sit at |latitude| > 45°, so a geocentric stand-in would land up to ~15 km from the requested surface — some 30× the ~0.5 km the oracle's own iteration converges to, and in the one place where an altitude error maps directly onto an atmospheric density.

hemisphere follows IRBEM's hemi_flag. Hemisphere::Same walks in the direction of increasing |B|, which is the input point's own side of the magnetic equator; Hemisphere::Opposite walks the other way. Hemisphere::North and Hemisphere::South name a side outright: the same-hemisphere foot is traced first and, if its geodetic latitude has the wrong sign, the other direction is traced instead. Classifying by the FOOT's latitude rather than the start's is what makes the answer right for a start point below the magnetic equator but above the geographic one — the two differ by up to 11° and the dip equator wanders further still.

The half-line the input point sits on, sampled uniformly for a plot. The direction is the input point's own magnetic hemisphere — the direction of increasing |B| — which is the reference's behaviour and the only one that makes "toward the Earth" well defined for a point off the magnetic equator. Sample 0 is the input point verbatim; the last sample is the first one inside PathTraceOptions::r0, so the path visibly crosses the surface rather than stopping short of it, which is what the reference does and what a plot wants.

The step is PathTraceOptions::step_size when positive and L/steps_per_l otherwise — the reference takes ds as a required argument here and nowhere else, and this is why.

Samples run along +B̂: index 0 is the foot the field points away from (the southern one, for the real geomagnetic field), the last index is the foot it points into, and the input point sits verbatim at TracedLine::start_index. The reference orders its posit the other way; the choice is arbitrary, so it is stated rather than inherited, and a caller who cares should read start_index instead of assuming an end.

Both ends are capped ON the reference surface by regula falsi inside the last step, so the path's extent is a property of the field line and not of the step size.

Lm, Bmin and XJ come from the samples already in the caller's buffer rather than from a second trace — one pass over the line, not two. XJ is I for a particle mirroring at the input point, integrated by the same right-endpoint rule lstar.hpp uses outward from the equator in both directions, so this and trace_invariant are the same quadrature on the same grid rather than two approximations that happen to be close. (IrbemTraceApi.XjAgreesWithTheInvariantTracer holds them to 1 × 10⁻⁴ relative; what separates them is only that RK4 is not exactly reversible, so the two grids drift apart at O(ds⁵) per step.)

The batch form of trace_field_line_toward_earth, and the only shape a device can accelerate: one trace is a serial RK4 chain, so no hardware makes it faster, and the parallelism is entirely ACROSS lines.

This kernel's crossover is 256 lines — MEASURED, and lower than the invariant tracer's 512 rather than higher, which is not what the output shape predicts. Both run the same RK4 chain over the same field. irbem_trace_i_f32 returns four floats per line — ~9 400 flops per byte moved. This one returns four floats per STEP, so at ~250 steps a line it moves ~250× as many bytes for the same arithmetic: ~125 flops/byte. The transfer term the invariant tracer never pays is real, but it does not move the crossover, because a crossover is a RATIO and this routine's host lane is ~2.6× dearer per line too. Where the transfer shows up instead is the ceiling — see the two bullets under the table. That is why it is a separate kernel with a separate registry row and a separate measured threshold, and why merging the two would be a mistake in both directions.

Measured on an RTX 3070 Ti against this header's own fp64 host lane (-O3 -march=native -ffp-contract=off), fixed step ds = 0.02 R_E, 512 samples of headroom per line, starts spread over L = 2…8 (IrbemTraceApiGpu.PathKernelCrossover):

against a host lane flat at 154–165 µs/line. Every cell is a range over FOUR full runs, not a best-of: below 1 024 lines the spread is a few per cent and the crossover reproduces at 256 every time, while above 4 096 the same size has measured 9.6× and 23.5× on different runs. Two things in that row are worth stating plainly because they contradict the obvious expectation:

• The crossover came out LOWER than the invariant tracer's, not higher. Not because this kernel is cheaper — because the HOST lane is dearer. A fixed ds = 0.02 R_E line is ~250 steps against the invariant tracer's ~120 at ds = L/50, so the host costs 154 µs/line against 60, and the same submit floor is paid off in half the batch. A crossover is a ratio.

• Where the bandwidth shows is the CEILING. The invariant tracer's speedup is still climbing at 65 536 lines (48.9×); this one flattens in the low-to-mid twenties and goes run-to-run noisy above 4 096, where a batch stages 100 MB–1 GB of results through host memory. The ranges above are quoted as ranges on purpose: at 16 384 lines the spread across four runs is a factor of 2.4, so any single number from that end of the curve is a performance claim the next run will not support.