Speed up propagation step time by 3x

[scottshambaugh]

I've been trying out different numerical propagators, and this is the fastest so far, with results aligning with both Orekit and Tudat. Nice project!

This PR rolls in the following performance improvements:

• Compile gravity twice, the second time with fma (fused multiply add) compiler support, and use fma at runtime if the machine supports it. Cuts the gravity evaluation by ~5x
• Pull out some MSIS calculations that were being computed multiple times per step
• The 5-iteration Bowring geodetic conversion converged in a single refinement step in my testing, so drop the loop
• Reuse already fetched space weather
• Faster version string import (skip importlib)
• Skip fully parsing gravity term data higher than the requested order

One more performance issue to flag: the initial timestep PID convergence is slow, and is the dominant time on short-duration propagation - about half the evaluations for a one-hour horizon. A warm start or forced initial step would drastically reduce the needed calls there. I believe this would require a coordinated numeris change.

benchmark	before	after
earthgravity/accel 40x40	13.6 us	2.6 us
nrlmsise/density_400km	2.48 us	1.80 us
30 days prop (the script below)	4.8 s	1.5 s

import time
import numpy as np
import satkit as sk
settings = sk.propsettings()
settings.gravity_model, settings.gravity_degree, settings.gravity_order = sk.gravmodel.egm96, 40, 40
settings.use_sun_gravity = settings.use_moon_gravity = settings.use_relativistic_correction = True
settings.tide_model, settings.integrator = sk.tidemodel.solid_step1, sk.integrator.rkv98
settings.abs_error = settings.rel_error = 1e-9
settings.use_spaceweather = True
mass_kg, area_m2, cd, cr = 500.0, 5.0, 2.2, 1.2
props = sk.satproperties(cdaoverm=cd * area_m2 / mass_kg, craoverm=cr * area_m2 / mass_kg)
epoch = sk.time(2025, 1, 1, 12, 0, 0)
altitude_m, inclination_rad = 550e3, np.radians(51.6)
r = sk.consts.earth_radius + altitude_m
v = np.sqrt(sk.consts.mu_earth / r)
state = np.array([r, 0.0, 0.0, 0.0, v * np.cos(inclination_rad), v * np.sin(inclination_rad)])
sk.propagate(state, epoch, duration_secs=60.0, propsettings=settings, satproperties=props)
t = time.perf_counter()
for day in range(30):
    res = sk.propagate(state, epoch + sk.duration.from_seconds(day * 86400.0), duration_secs=86400.0, propsettings=settings, satproperties=props)
    state = np.asarray(res.state)
print(f"30 days in {time.perf_counter() - t:.2f} s, final state {state.tolist()}")

[scottshambaugh]

Space weather timing bug broken out to #176 instead of being lumped in with this.

[ssmichael1]

This is a very impressive speedup. I'll take a closer look over the weekend, and merge if performance/accuracy is expected. Thank you.

[ssmichael1]

For your comment about the initial timestep PID convergence: I'll have to think about this one. If I remember correctly, I had a tough time figuring this out and ended up copying the algorithm for initial step size from the ODE module in the Julia programming language. I could:

1. Allow user-settable initial step; this is an easy win.
2. Look at adjusting the PID coefficients (but could have big downstream effects; nervous about tuning)
3. Modify initial step size calculation to be less conservative.

I'll look at this over the weekend as well, and thanks again for your contributions.

[ssmichael1]

All looks good. Merged. Thank you!