rho_125 Function

   function rho_125(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) :: f4, q1, q2, q3, q4, q5, q6
      real(dp) :: rho_prime, x_star, u(2), w(2), v, factor_k, t_100
      real(dp) :: rho_sum, rhoi, roots_2, log_f3
      integer(ip) :: i

      ! Compute base density polynomial
      rho_prime = ZETA_CON(7)
      do i = 6, 1, -1
         rho_prime = rho_prime * me%t_infinity + ZETA_CON(i)
      end do

      ! Compute base temperature
      t_100 = me%tx + OMEGA * (me%tx - TZERO)

      ! 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 powers
      q2 = 1.0_dp / u(1)
      q3 = -1.0_dp / u(2)
      q5 = 1.0_dp / v
      q4 = (1.0_dp + w(1) * q2 + w(2) * q3 + me%root1 * me%root2 * &
           (me%cb_polar_squared - roots_2) * q5) / x_star
      q1 = -2.0_dp * q4 - q3 - q2
      q6 = -q5 - 2.0_dp * (me%x_root + me%cb_polar_radius) * q4 - &
           (me%root2 + me%cb_polar_radius) * q3 - &
           (me%root1 + me%cb_polar_radius) * q2

      ! Compute log of f3 function
      log_f3 = q1 * log((height + me%cb_polar_radius) / (100.0_dp + me%cb_polar_radius)) + &
               q2 * log((height - me%root1) / (100.0_dp - me%root1)) + &
               q3 * log((height - me%root2) / (100.0_dp - me%root2)) + &
               q4 * log((height**2 - 2.0_dp * me%x_root * height + roots_2) / &
                       (1.0e4_dp - 200.0_dp * me%x_root + roots_2))

      ! Compute f4 function
      f4 = (height - 100.0_dp) * q5 / ((height + me%cb_polar_radius) * &
           (100.0_dp + me%cb_polar_radius)) + &
           q6 * atan(me%y_root * (height - 100.0_dp) / &
           (me%y_root**2 + (height - me%x_root) * (100.0_dp - me%x_root))) / &
           me%y_root

      ! Compute factor_k (diffusion scale height coefficient)
      factor_k = -1500625.0_dp * G_ZERO * me%cb_polar_squared / &
                 (GAS_CON * CON_C(5) * (me%tx - TZERO))

      ! Compute mass-dependent sum
      rho_sum = 0.0_dp
      do i = 1, 5
         rhoi = MOL_MASS(i) * NUM_DENS(i) * exp(MOL_MASS(i) * factor_k * (f4 + log_f3))

         if (i == 3) then  ! Helium
            rhoi = rhoi * (t_100 / temperature)**(-0.38_dp)
         end if

         rho_sum = rho_sum + rhoi
      end do

      density = rho_sum * rho_prime * t_100 / temperature

   end function rho_125