Fix DDH analytic STIGMA delay derivative scaling by T_⊙ ​

[vhaasteren]

Summary

• Remove the spurious / Tsun.value factors from analytic ∂delay/∂STIGMA in DDHmodel.d_delayS_d_par (bug since Maple rewrite 733cf704, PINT ≥ 1.0).
• Rewrite the STIGMA branch in terms of the shared Shapiro logNum / geo / dSINI/dSTIGMA so it cannot drift from the delay expression again.
• Add test_ddh_stigma_derivative_matches_finite_difference: analytic vs d_delay_d_param_num RMS ratio for STIGMA and H3 (control).

The delay itself, ∂/∂H3, and ELL1H STIGMA derivatives were already correct. This only repairs the analytic DDH STIGMA design-matrix column (too large by 1/T⊙, which had inflated Fisher information and produced absurdly small formal σ(STIG) in GLS fits with free STIGMA.

Ping here for @dlakaplan, who originally derived the derivative so best to confirm

NANOGrav was not using this derivative in their work. Other PTAs would, but they are using the tempo2 equivalent.

Test plan

• pytest tests/test_ddh.py::test_ddh_stigma_derivative_matches_finite_difference -v
• black --check on touched files
• CI green

Minimal reproduction

import io
import numpy as np
from astropy import units as u
import pint.binaryconvert
import pint.derived_quantities
import pint.simulation
from pint import Tsun
from pint.models import get_model
Mp, Mc, i, PB = 1.4 * u.Msun, 1.1 * u.Msun, 85 * u.deg, 0.5 * u.day
A1 = pint.derived_quantities.a1sini(Mp, Mc, PB, i)
par = f"""
PSRJ X
RAJ 18:57:36
DECJ +09:43:17
F0 186.5
PEPOCH 49453
DM 13
BINARY DD
PB {PB.value}
A1 {A1.value}
T0 49452.94
OM 276.55
ECC 0.1
SINI {np.sin(i).value}
M2 {Mc.value}
TZRMJD 54177.5
TZRFRQ 424
TZRSITE ao
CLK TT(TAI)
UNITS TDB
EPHEM DE405
TIMEEPH FB90
"""
m = get_model(io.StringIO(par))
t = pint.simulation.make_fake_toas_uniform(
    50000, 51000, 200, m, add_noise=False, error=1 * u.us
)
m2 = pint.binaryconvert.convert_binary(m, "DDH")
an = m2.d_delay_d_param(t, "STIGMA")
num = m2.d_delay_d_param_num(t, "STIGMA")
ratio = np.sqrt(np.mean(an**2)) / np.sqrt(np.mean(num**2))
print("analytic/numerical", ratio)
print("1/Tsun.value", 1 / Tsun.value)
# Expect (buggy): ratio ≈ 1/Tsun.value ≈ 2.03e5
# After fix: ratio ≈ 1

[dlakaplan]

I could believe that there was an error there, although I don't have my old derivations handy. But especially if this passes the analytic vs numerical test I think it's likely in good shape.