I'm writing a Python script to calculate the effect of atmospheric drag on a spacecraft in LEO by evaluating the change in the square root of the semi-major axis over time. The governing equation I'm using is:
np.sqrt(a_final) - np.sqrt(a_initial) = (- np.sqrt(mu) * C_d * A * p * delta_time) / (2 * m)
Where: a_initial = initial semi major axis, a_final = final semi major axis, mu is the gravitational parameter (m³/s²), C_d is the drag coefficient (dimensionless), A is the cross-sectional area (m²), m is the spacecraft mass (kg), p is the atmospheric density at altitude (kg/m³), delta_time is the time interval in seconds.
Ballistic_coefficient = (C_d * A) / m # (m²/kg)
B* = (Ballistic_coefficient * p_0) / 2 #(1/m)
where: p_0 is the reference atmospheric density.
Here's the confusion: Celestrak defines p₀ as a reference atmospheric density, definition not clear. However, the Spacetrack Report (PDF, page 31) uses p₀ in kg/m², leading to a mismatch in units and confusion when extracting Ballistic_coefficient from B*.
What is the correct value and unit of p₀ that I should use in the equation to extract a ballistic coefficient in SI units (m²/kg) from a given TLE B* value, so that I can use it correctly in the semi-major axis decay formula?
