rho_high Function

   function rho_high(me, height, temperature, t_500, sun_dec, geo_lat) result(density)
      class(jacchia_roberts_type), intent(inout) :: me
      real(dp), intent(in) :: height !! Spacecraft altitude (km)
      real(dp), intent(in) :: temperature !! Exospheric temperature at spacecraft altitude (K)
      real(dp), intent(in) :: t_500 !! Exospheric temperature at altitude of 500 km (K)
      real(dp), intent(in) :: sun_dec !! Declination of the sun in TOD coordinates (radians)
      real(dp), intent(in) :: geo_lat !! Geodetic latitude of spacecraft (radians)
      real(dp) :: density

      real(dp) :: f, log_di, gamma, exp1, r
      real(dp) :: di, polar125
      integer(ip) :: i, j

      real(dp), parameter :: SIN_PI_OVER_4_CUBED = sin(0.25_dp * PI)**3

      density = 0.0_dp
      polar125 = me%cb_polar_radius + 125.0_dp

      do i = 1, 6
         ! Compute constituent density sum for this atmospheric component
         if (i <= 5) then  ! Skip hydrogen (i=6) for initial calculation
            log_di = CON_DEN(i,7)
            do j = 6, 1, -1
               log_di = log_di * me%t_infinity + CON_DEN(i,j)
            end do
            di = 10.0_dp**log_di / AVOGADRO
         end if

         ! Compute second exponent in density expression
         gamma = 35.0_dp * MOL_MASS(i) * G_ZERO * me%cb_polar_squared * &
                 (me%t_infinity - me%tx) / &
                 (GAS_CON * me%sum * me%t_infinity * &
                  (me%tx - TZERO) * polar125)

         ! Compute first exponent in density expression
         exp1 = 1.0_dp + gamma

         ! A factor which is non-unity only for helium
         f = 1.0_dp

         ! Compute corrections for helium
         if (i == 3) then  ! Helium
            exp1 = exp1 - 0.38_dp
            ! Handle sun_dec = 0 case to avoid NaN from 0/0
            if (abs(sun_dec) < 1.0e-15_dp) then
               ! When sun_dec = 0, the formula simplifies to f = 10^0 = 1
               f = 1.0_dp
            else
               f = 4.9914_dp * abs(sun_dec) * &
                  (sin(0.25_dp * PI - 0.5_dp * geo_lat * sun_dec / abs(sun_dec))**3 - &
                  SIN_PI_OVER_4_CUBED) / PI
               f = 10.0_dp**f
            end if
         end if

         ! Add corrections to density, skip hydrogen unless above 500 km
         if (height > 500.0_dp .and. i == 6) then
            r = MOL_MASS(6) * 10.0_dp**(73.13_dp - (39.4_dp - 5.5_dp * log10(t_500)) * &
                log10(t_500)) * (t_500 / temperature)**exp1 * &
                ((me%t_infinity - temperature) / (me%t_infinity - t_500))**gamma / &
                AVOGADRO
            density = density + r
         else if (i <= 5) then
            r = f * MOL_MASS(i) * di * (me%tx / temperature)**exp1 * &
                ((me%t_infinity - temperature) / (me%t_infinity - me%tx))**gamma
            density = density + r
         end if
      end do

   end function rho_high