rho_cor Function

   pure function rho_cor(height, utc_mjd, geo_lat, geo) result(correction)
      real(dp), intent(in) :: height !! Spacecraft altitude (km)
      real(dp), intent(in) :: utc_mjd !! UTC Modified Julian Date
      real(dp), intent(in) :: geo_lat !! Geodetic latitude of spacecraft (radians)
      type(geoparms_type), intent(in) :: geo !! Geomagnetic parameters
      real(dp) :: correction !! Correction factor (multiplicative)

      real(dp) :: geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha
      real(dp) :: sin_lat, eta_lat

      ! Compute geomagnetic activity correction
      if (height < 200.0_dp) then
         geo_cor = 0.012_dp * geo%tkp + 0.000012_dp * exp(geo%tkp)
      else
         geo_cor = 0.0_dp
      end if

      ! Compute semiannual variation correction
      f = (5.876e-7_dp * height**2.331_dp + 0.06328_dp) * exp(-0.002868_dp * height)

      ! Compute years since Jan 1, 1958 midnight (MJD 36204.0)
      ! Note: The C++ code uses 6204.5 because GMAT uses GSFC MJD (reference epoch
      ! Jan 5, 1941 noon = MJD 29999.5). In GSFC MJD, Jan 1, 1958 midnight = 6204.5.
      ! This Fortran implementation uses standard MJD where Jan 1, 1958 midnight = 36204.0.
      ! Verified: 36204.0 - 29999.5 = 6204.5
      day_58 = (utc_mjd - 36204.0_dp) / 365.2422_dp

      tausa = day_58 + 0.09544_dp * &
              ((0.5_dp * (1.0_dp + sin(2.0_dp * PI * day_58 + 6.035_dp)))**1.65_dp - 0.5_dp)
      alpha = sin(4.0_dp * PI * tausa + 4.259_dp)
      g = 0.02835_dp + (0.3817_dp + 0.17829_dp * sin(2.0_dp * PI * tausa + 4.137_dp)) * alpha
      semian_cor = f * g

      ! Compute seasonal latitudinal variation
      sin_lat = sin(geo_lat)
      eta_lat = sin(2.0_dp * PI * day_58 + 1.72_dp) * sin_lat * abs(sin_lat)
      slat_cor = 0.014_dp * (height - 90.0_dp) * eta_lat * &
                 exp(-0.0013_dp * (height - 90.0_dp)**2)

      correction = 10.0_dp**(geo_cor + semian_cor + slat_cor)

   end function rho_cor