rho_100 Function

   function rho_100(me, height, temperature) result(density)
      class(jacchia_roberts_type), intent(inout) :: me
      real(dp), intent(in) :: height !! Spacecraft altitude (km)
      real(dp), intent(in) :: temperature !! Exospheric temperature (K)
      real(dp) :: density !! Atmospheric density (g/cm^3)

      real(dp) :: m_poly, s_poly, f2, p1, p2, p3, p4, p5, p6, b(6)
      real(dp) :: x_star, u(2), w(2), v, factor_k, roots_2, log_f1
      integer(ip) :: i

      logical,parameter :: use_vallado = .false. !! Flag to control whether to use Vallado's
                                                 !! `factor_k` equation instead of original GDTS one.

      ! Compute M(z) polynomial
      m_poly = M_CON(7)
      do i = 6, 1, -1
         m_poly = m_poly * height + M_CON(i)
      end do

      ! Compute temperature dependent coefficients
      do i = 1, 6
         b(i) = S_CON(i) + S_BETA(i) * me%tx / (me%tx - TZERO)
      end do

      ! Compute functions of auxiliary temperature values
      roots_2 = me%x_root**2 + me%y_root**2
      x_star = -2.0_dp * me%root1 * me%root2 * me%cb_polar_radius * &
               (me%cb_polar_squared + 2.0_dp * me%cb_polar_radius * &
                me%x_root + roots_2)
      v = (me%cb_polar_radius + me%root1) * &
          (me%cb_polar_radius + me%root2) * &
          (me%cb_polar_squared + 2.0_dp * me%cb_polar_radius * &
           me%x_root + roots_2)
      u(1) = (me%root1 - me%root2) * &
             (me%root1 + me%cb_polar_radius)**2 * &
             (me%root1**2 - 2.0_dp * me%root1 * me%x_root + roots_2)
      u(2) = (me%root1 - me%root2) * &
             (me%root2 + me%cb_polar_radius)**2 * &
             (me%root2**2 - 2.0_dp * me%root2 * me%x_root + roots_2)
      w(1) = me%root1 * me%root2 * me%cb_polar_radius * &
             (me%cb_polar_radius + me%root1) * &
             (me%cb_polar_radius + roots_2 / me%root1)
      w(2) = me%root1 * me%root2 * me%cb_polar_radius * &
             (me%cb_polar_radius + me%root2) * &
             (me%cb_polar_radius + roots_2 / me%root2)

      ! Compute S(z) polynomial for z = root1
      s_poly = b(6)
      do i = 5, 1, -1
         s_poly = s_poly * me%root1 + b(i)
      end do
      p2 = s_poly / u(1)

      ! Compute S(z) polynomial for z = root2
      s_poly = b(6)
      do i = 5, 1, -1
         s_poly = s_poly * me%root2 + b(i)
      end do
      p3 = -s_poly / u(2)

      ! Compute S(z) polynomial for z = negative earth radius
      s_poly = b(6)
      do i = 5, 1, -1
         s_poly = -s_poly * me%cb_polar_radius + b(i)
      end do
      p5 = s_poly / v

      ! Compute power of fourth quantity in f1 function
      p4 = (b(1) - me%root1 * me%root2 * me%cb_polar_squared * &
           (b(5) + b(6) * (2.0_dp * me%x_root + me%root1 + &
            me%root2 - me%cb_polar_radius)) + w(1) * p2 + w(2) * p3 - &
            me%root1 * me%root2 * b(6) * me%cb_polar_radius * roots_2 + &
            me%root1 * me%root2 * (me%cb_polar_squared - roots_2) * p5) / x_star

      ! Compute power of first quantity in f1 function
      p1 = b(6) - 2.0_dp * p4 - p3 - p2

      ! Compute p6 factor in f2 function
      p6 = b(5) + b(6) * (2.0_dp * me%x_root + me%root1 + &
           me%root2 - me%cb_polar_radius) - p5 - &
           2.0_dp * (me%x_root + me%cb_polar_radius) * p4 - &
           (me%root2 + me%cb_polar_radius) * p3 - &
           (me%root1 + me%cb_polar_radius) * p2

      ! Compute natural log of f1 function
      log_f1 = p1 * log((height + me%cb_polar_radius) / (90.0_dp + me%cb_polar_radius)) + &
               p2 * log((height - me%root1) / (90.0_dp - me%root1)) + &
               p3 * log((height - me%root2) / (90.0_dp - me%root2)) + &
               p4 * log((height**2 - 2.0_dp * me%x_root * height + roots_2) / &
                       (8100.0_dp - 180.0_dp * me%x_root + roots_2))

      ! Compute f2 function
      f2 = (height - 90.0_dp) * (M_CON(7) + p5 / ((height + me%cb_polar_radius) * &
           (90.0_dp + me%cb_polar_radius))) + &
           p6 * atan(me%y_root * (height - 90.0_dp) / &
           (me%y_root**2 + (height - me%x_root) * (90.0_dp - me%x_root))) / &
           me%y_root

      ! Compute f1 power
      if (use_vallado) then
         ! ==== GMAT FORMULA (INCORRECT?) ====
         ! Based on Vallado (3rd Ed, p 951),
         ! simplifies this equation by removing the factor
         ! flc4 = 1500625*cb_polar_squared/CON_C(5).
         !
         ! WHY IT'S WRONG:
         ! Roberts (1971) had this factor INSIDE the F2 formula:
         !   F2 = (h-90) * (1500625*Ra²/Ca[5] * M_CON[7] + ...)
         ! GTDS refactored by factoring it OUT of F2 and multiplying it into k:
         !   factor_k = fkl * flc4 = (-g0/(R*(Tx-T0))) * (1500625*Ra²/CON_C(5))
         !   F2 = (h-90) * (M_CON[7] + ...)  [factor removed]
         ! These are mathematically equivalent.
         !
         ! Removing the factor from BOTH places leaves it completely absent.
         ! This causes density errors of 143-244% at 95-97 km vs numerical integration.
         factor_k = -G_ZERO / (GAS_CON * (me%tx - TZERO))
      else
         ! ==== ORIGINAL GTDS FORMULA (CORRECT?) ====
         ! This is the formula from GTDS (it is commented out in GMAT) and matches Roberts
         ! when accounting for the refactoring (factor moved from F2 to factor_k).
         ! Empirically validated: matches INPE 5-point Newton-Cotes numerical
         ! integration to within 0.03-0.4% throughout the 90-100 km range.
         factor_k = -1500625.0_dp * G_ZERO * me%cb_polar_squared / &
                    (GAS_CON * CON_C(5) * (me%tx - TZERO))
      end if

      ! Compute final density
      density = RHO_ZERO * TZERO * m_poly * exp(factor_k * (log_f1 + f2)) / &
                (MZERO * temperature)

   end function rho_100