find_cstar_roots_original Subroutine

public pure subroutine find_cstar_roots_original(c_star, tx, root1, root2, x_root, y_root)

Original version of root-finding routine using the roots subroutine. This is based on the routine from GMAT.

Problems with this method:

1. Sequential deflation accumulates errors: After finding each real root, the polynomial is deflated (reduced in degree). Each deflation introduces numerical errors that compound, corrupting the polynomial coefficients.

2. Massive polynomial errors: While the real roots (root1, root2) found first are accurate (~1e-8 error), the complex conjugate pair computed from the twice-deflated polynomial has errors of 50-700 (absolute polynomial evaluation), despite the roots routine converging.

3. roots_2 errors of 70-76%: At extreme temperatures (Tx ~ 300 K or Tx ~ 550-600 K), the quantity roots_2 = x_root² + y_root² differs by 70-76% from the correct value due to deflation corruption.

4. Severe density errors in 100-125 km range: The incorrect complex roots propagate through exponential terms in rho_100 and rho_125, causing density overestimation by factors of 2-22× at low solar activity (F10.7 ~ 50). At 125 km, the old method predicts 22× higher density than reality.

5. Discontinuity at 125 km boundary: The density jumps by 37× in 5 km (125->130 km) because rho_high (>125 km) doesn't use complex roots and produces correct values, while rho_125 uses corrupted roots.

Why the real roots still work:

The real roots are found BEFORE deflation corrupts the polynomial, so they remain accurate. Only the complex roots, found last from the twice-deflated polynomial, are severely corrupted.

Solution:

Use find_cstar_roots instead, which uses simultaneous Newton-Raphson with temperature-dependent initial guesses (Kuga/INPE 1985) and computes complex roots analytically via Vieta's formulas, avoiding deflation entirely. Achieves ~1e-8 accuracy for all roots across all physical temperatures.

Arguments

   pure subroutine find_cstar_roots_original(c_star, tx, root1, root2, x_root, y_root)

      use jacchia_roberts_utilities, only: roots, deflate_polynomial

      real(dp), intent(inout)  :: c_star(5) !! Quartic polynomial coefficients `c_star(1)..c_star(5)` [note: modified in-place by deflation]
      real(dp), intent(in)  :: tx        !! Temperature at 125 km (K) [not used here]
      real(dp), intent(out) :: root1     !! Larger real root (> 125 km)
      real(dp), intent(out) :: root2     !! Smaller real root (< 100 km)
      real(dp), intent(out) :: x_root    !! Real part of complex conjugate root pair
      real(dp), intent(out) :: y_root    !! Imaginary part of complex conjugate root pair

      real(dp) :: aux(4,2)
      integer(ip) :: na

      na = 5

      ! Get 1st real root
      aux(1,1) = 125.0_dp
      aux(1,2) = 0.0_dp
      call roots(c_star, na, aux, 1)
      root1 = aux(1,1)

      call deflate_polynomial(c_star, na, root1, c_star)

      ! Get 2nd real root
      aux(1,1) = 200.0_dp
      aux(1,2) = 0.0_dp
      call roots(c_star, na-1, aux, 1)
      root2 = aux(1,1)

      call deflate_polynomial(c_star, na-1, root2, c_star)

      ! Get remaining complex roots
      aux(1,1) = 10.0_dp
      aux(1,2) = 125.0_dp
      call roots(c_star, na-2, aux, 1)
      x_root = aux(1,1)
      y_root = abs(aux(1,2))

   end subroutine find_cstar_roots_original