jacchia_roberts_density Function

   function jacchia_roberts_density(me, height, position, sun_vector, geo_lat, &
                                    utc_mjd) result(density)
      class(jacchia_roberts_type), intent(inout) :: me
      real(dp), intent(in) :: height !! Spacecraft height above reference ellipsoid (km)
      real(dp), intent(in) :: position(3) !! Spacecraft position vector (km, TOD/GCI)
      real(dp), intent(in) :: sun_vector(3) !! Unit vector to Sun (TOD/GCI)
      real(dp), intent(in) :: geo_lat !! Geodetic latitude (degrees)
      real(dp), intent(in) :: utc_mjd !! UTC Modified Julian Date
      real(dp) :: density !! Atmospheric density (kg/m^3)

      real(dp) :: sun_dec !! Sun declination (radians)
      real(dp) :: temperature !! Local exospheric temperature (K)
      real(dp) :: t_500 !! Temperature at 500 km (K)
      real(dp) :: geo_lat_rad !! Geodetic latitude in radians
      real(dp) :: raw_density !! Density before unit conversion (g/cm^3)
      real(dp) :: f107_daily !! Daily F10.7 value (adjusted for measurement time)
      real(dp) :: f107a_avg !! 81-day average F10.7a (adjusted for measurement time)
      real(dp) :: tkp !! Kp index
      type(geoparms_type) :: geo !! Geomagnetic parameters
      type(flux_data_type) :: flux_data !! Space weather flux data
      logical :: sw_status !! Status for space weather data retrieval

      ! sun declination: atan2(Z, sqrt(X^2 + Y^2))
      sun_dec = atan2(sun_vector(3), sqrt(sun_vector(1)**2 + sun_vector(2)**2))

      if (me%use_sw_file) then
         ! Get space weather data for this date
         call me%sw_data%get_flux_data(utc_mjd, flux_data, sw_status)

         if (.not. sw_status) then
            write(*,*) 'WARNING: Using nominal space weather values'
            ! Use nominal values as fallback
            geo%xtemp = me%parms_to_temp(nominalF107, nominalF107a, nominalKp)
            geo%tkp = nominalKp
         else
            ! Prepare flux data with proper timing adjustments.
            ! Kp: 6.7hr lag (matches PrepareKpData). F10.7: previous day (matches PrepareKpData).
            ! F10.7a: detected day (matches PrepareApData intent; see prepare_flux_data docstring).
            call prepare_flux_data(me%sw_data, flux_data, utc_mjd, tkp, f107_daily, f107a_avg)
            ! Calculate exospheric temperature from adjusted F10.7 data
            geo%xtemp = me%parms_to_temp(f107_daily, f107a_avg, tkp)
            ! Set Kp index in geomagnetic parameters:
            geo%tkp = tkp
         end if

      else
         ! Use specified override values
         geo%xtemp = me%parms_to_temp(me%f107_override, me%f107a_override, me%kp_override)
         geo%tkp = me%kp_override
      end if

      ! Convert geodetic latitude to radians
      geo_lat_rad = geo_lat * RAD_PER_DEG

      ! Compute height-dependent density
      if (height <= 90.0_dp) then
         ! At or below 90 km: use constant density
         raw_density = RHO_ZERO
      else if (height < 100.0_dp) then
         ! 90-100 km: use rho_100
         temperature = me%exotherm(position, sun_vector, geo, height, &
                                   sun_dec, geo_lat_rad)
         raw_density = me%rho_100(height, temperature)
      else if (height <= 125.0_dp) then
         ! 100-125 km: use rho_125
         temperature = me%exotherm(position, sun_vector, geo, height, &
                                   sun_dec, geo_lat_rad)
         raw_density = me%rho_125(height, temperature)
      else if (height <= 2500.0_dp) then
         t_500 = me%exotherm(position, sun_vector, geo, 500.0_dp, &
                             sun_dec, geo_lat_rad)
         temperature = me%exotherm(position, sun_vector, geo, height, &
                                   sun_dec, geo_lat_rad)
         raw_density = me%rho_high(height, temperature, t_500, sun_dec, geo_lat_rad)
      else
         raw_density = 0.0_dp
      end if

      ! Apply density correction
      raw_density = raw_density * rho_cor(height, utc_mjd, geo_lat_rad, geo)

      ! Convert from g/cm^3 to kg/m^3
      density = raw_density * 1.0e3_dp

   end function jacchia_roberts_density