exotherm Function

   function exotherm(me, position, sun_vector, geo, height, sun_dec, geo_lat) &
                     result(exotemp)
      class(jacchia_roberts_type), intent(inout) :: me
      real(dp), intent(in) :: position(3) !! Spacecraft position (km)
      real(dp), intent(in) :: sun_vector(3) !! Sun unit vector
      type(geoparms_type), intent(in) :: geo !! Geomagnetic parameters
      real(dp), intent(in) :: height !! Spacecraft height (km)
      real(dp), intent(in) :: sun_dec !! Sun declination (radians)
      real(dp), intent(in) :: geo_lat !! Geodetic latitude (radians)
      real(dp) :: exotemp !! Local exospheric temperature (K)

      real(dp) :: expkp, hour_angle, cross_denom
      real(dp) :: theta, eta, tau, th22, t1, sun_denom, cos_denom
      real(dp) :: c_star(5)
      real(dp) :: cosAlpha
      integer(ip) :: i

      real(dp), parameter :: error_tolerance = 1.0e-14_dp
      real(dp), parameter :: real_tol = 1.0e-15_dp
      real(dp), parameter :: d37 = 37.0_dp * RAD_PER_DEG  !! 37 degrees in radians
      real(dp), parameter :: d6  = 6.0_dp  * RAD_PER_DEG  !! 6 degrees in radians
      real(dp), parameter :: d43 = 43.0_dp * RAD_PER_DEG  !! 43 degrees in radians
      real(dp), parameter :: a    = 371.6678_dp !! Constant for temperature calculation
      real(dp), parameter :: b    = 0.0518806_dp !! Constant for temperature calculation
      real(dp), parameter :: c    = -294.3505_dp !! Constant for temperature calculation
      real(dp), parameter :: kbar = -0.00216222_dp !! Constant for temperature calculation

      logical,parameter :: use_old_root_method = .false. !! to use the old (gdts) root method,
                                                         !! which seems to have some problems.

      ! Compute hour angle of the sun
      sun_denom = sqrt(sun_vector(1)**2 + sun_vector(2)**2)
      cos_denom = sqrt(position(1)**2 + position(2)**2)

      ! Check for invalid position (zero magnitude in xy-plane)
      if (cos_denom < real_tol) then
         write(*,*) 'ERROR: Position vector has zero magnitude in xy-plane'
         exotemp = TZERO
         return
      end if

      ! Check for degenerate sun vector (sun at coordinate pole, XY magnitude = 0).
      ! For Earth-Sun geometry this cannot happen (sun declination <= 23.4 deg gives
      ! sun_denom >= 0.92), but guard against it to avoid NaN from 0/0.
      if (sun_denom < real_tol) then
         hour_angle = 0.0_dp
      else

         ! Calculate cosine of angle between vectors
         cosAlpha = (sun_vector(1) * position(1) + sun_vector(2) * position(2)) / &
                  (sun_denom * cos_denom)

         ! Compute cross product for sign determination
         cross_denom = abs(sun_vector(1) * position(2) - sun_vector(2) * position(1))

         ! Calculate hour angle with boundary checking
         ! When vectors are nearly collinear, cosAlpha determines the angle directly
         if (cosAlpha >= 1.0_dp - error_tolerance) then
            ! Vectors aligned: hour angle = 0
            hour_angle = 0.0_dp
         else if (cosAlpha <= -1.0_dp + error_tolerance) then
            ! Vectors opposite: hour angle = +/- PI
            ! If cross product is too small (collinear case), default to +PI
            if (cross_denom < real_tol) then
               hour_angle = PI
            else
               hour_angle = sign(1.0_dp, sun_vector(1) * position(2) - &
                              sun_vector(2) * position(1)) * PI
            end if
         else
            ! General case: use cross product to determine sign
            ! If vectors are collinear (shouldn't happen here due to cosAlpha checks),
            ! default to 0
            if (cross_denom < real_tol) then
               hour_angle = 0.0_dp
            else
               hour_angle = sign(1.0_dp, sun_vector(1) * position(2) - &
                           sun_vector(2) * position(1)) * acos(cosAlpha)
            end if
         end if

      end if  ! sun_denom guard

      ! Compute sun and spacecraft position dependent part of temperature
      theta = 0.5_dp * abs(geo_lat + sun_dec)
      eta = 0.5_dp * abs(geo_lat - sun_dec)
      tau = hour_angle - d37 + d6 * sin(hour_angle + d43)
      if (tau < -PI) then
         tau = tau + 2.0_dp * PI
      else if (tau > PI) then
         tau = tau - 2.0_dp * PI
      end if

      th22 = sin(theta)**2.2_dp
      t1 = geo%xtemp * (1.0_dp + 0.3_dp * (th22 + cos(0.5_dp * tau)**3.0_dp * &
                        (cos(eta)**2.2_dp - th22)))
      expkp = exp(geo%tkp)

      ! Compute t_infinity based on altitude
      if (height < 200.0_dp) then
         me%t_infinity = t1 + 14.0_dp * geo%tkp + 0.02_dp * expkp
      else
         me%t_infinity = t1 + 28.0_dp * geo%tkp + 0.03_dp * expkp
      end if

      me%tx = a + b * me%t_infinity + c * exp(kbar * me%t_infinity)

      ! Compute temperature based on altitude regime
      if (height < 125.0_dp) then
         ! Compute height dependent polynomial
         me%sum = CON_C(5)
         do i = 4, 1, -1
            me%sum = CON_C(i) + me%sum * height
         end do

         ! Compute temperature
         exotemp = me%tx + (me%tx - TZERO) * me%sum / 1.500625e6_dp

      else if (height > 125.0_dp) then
         ! Compute temperature dependent polynomial
         me%sum = CON_L(5)
         do i = 4, 1, -1
            me%sum = CON_L(i) + me%sum * me%t_infinity
         end do

         ! Compute temperature
         exotemp = me%t_infinity - (me%t_infinity - me%tx) * &
                   exp(-(me%tx - TZERO) / (me%t_infinity - me%tx) * &
                       (height - 125.0_dp) / 35.0_dp * me%sum / &
                       (me%cb_polar_radius + height))
      else
         exotemp = me%tx
      end if

      ! Compute auxiliary quantities for heights <= 125 km
      if (height <= 125.0_dp) then
         ! Obtain coefficients of polynomial for auxiliary quantities
         c_star(1) = CON_C(1) + 1500625.0_dp * me%tx / (me%tx - TZERO)
         do i = 2, 5
            c_star(i) = CON_C(i)
         end do

         if (use_old_root_method) then
            ! original routine (gmat/gdts)
            call find_cstar_roots_original(c_star, me%tx, me%root1, me%root2, me%x_root, me%y_root)
         else
            ! this one seems better behaved
            call find_cstar_roots(c_star, me%tx, me%root1, me%root2, me%x_root, me%y_root)
         end if
      end if

   end function exotherm