program aerodynamic_heating_reentry
    implicit none
    integer :: iostat_val
    double precision :: altitude_km, V_inf_ms, Rn_m, emissivity
    double precision :: alt_m, T_inf_K, P_inf_Pa, rho_inf_kgm3, a_sound_ms, Mach_inf
    double precision :: q0_Wm2, q0_kWm2, q0_MWm2, Trad_K, Trad_C, h0_MJkg, P02_kPa
    double precision, parameter :: SIGMA = 5.670374d-8  ! W/(m2.K4)
    double precision, parameter :: K_SG  = 1.7415d-4   ! Sutton-Graves Earth constant
    double precision, parameter :: GAMMA = 1.4d0, R_AIR = 287.05d0

    ! Read inputs
    read(*,*,iostat=iostat_val) altitude_km   ! Flight Altitude [km] (e.g. 50.0)
    read(*,*,iostat=iostat_val) V_inf_ms      ! Reentry Velocity [m/s] (e.g. 6500.0)
    read(*,*,iostat=iostat_val) Rn_m          ! Nose Radius [m] (e.g. 0.50)
    read(*,*,iostat=iostat_val) emissivity    ! Surface Emissivity (e.g. 0.85)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for aerodynamic heating calculation.'
        stop
    end if

    if (altitude_km < 0.0d0 .or. altitude_km > 100.0d0 .or. V_inf_ms <= 0.0d0 .or. Rn_m <= 0.0d0) then
        write(*,*) 'ERROR: Altitude (0-100 km), velocity, and nose radius must be positive.'
        stop
    end if

    alt_m = altitude_km * 1000.0d0

    ! US Standard Atmosphere Model (Troposphere up to Mesosphere)
    if (altitude_km <= 11.0d0) then
        T_inf_K = 288.15d0 - 6.5d-3 * alt_m
        P_inf_Pa = 101325.0d0 * (T_inf_K / 288.15d0)**5.25588d0
    else if (altitude_km <= 20.0d0) then
        T_inf_K = 216.65d0
        P_inf_Pa = 22632.0d0 * exp(-9.80665d0 * 0.0289644d0 * (alt_m - 11000.0d0) / (8.31446d0 * 216.65d0))
    else if (altitude_km <= 32.0d0) then
        T_inf_K = 216.65d0 + 1.0d-3 * (alt_m - 20000.0d0)
        P_inf_Pa = 5474.87d0 * (216.65d0 / T_inf_K)**34.1632d0
    else if (altitude_km <= 47.0d0) then
        T_inf_K = 228.65d0 + 2.8d-3 * (alt_m - 32000.0d0)
        P_inf_Pa = 868.015d0 * (228.65d0 / T_inf_K)**12.2011d0
    else if (altitude_km <= 51.0d0) then
        T_inf_K = 270.65d0
        P_inf_Pa = 110.906d0 * exp(-9.80665d0 * 0.0289644d0 * (alt_m - 47000.0d0) / (8.31446d0 * 270.65d0))
    else if (altitude_km <= 71.0d0) then
        T_inf_K = 270.65d0 - 2.8d-3 * (alt_m - 51000.0d0)
        P_inf_Pa = 66.9389d0 * (T_inf_K / 270.65d0)**12.2011d0
    else
        T_inf_K = 214.65d0 - 2.0d-3 * (alt_m - 71000.0d0)
        if (T_inf_K < 165.0d0) T_inf_K = 165.0d0
        P_inf_Pa = 3.9564d0 * (T_inf_K / 214.65d0)**17.0816d0
    end if

    rho_inf_kgm3 = P_inf_Pa / (R_AIR * T_inf_K)
    a_sound_ms = sqrt(GAMMA * R_AIR * T_inf_K)
    Mach_inf = V_inf_ms / a_sound_ms

    ! Sutton-Graves Stagnation Point Convective Heat Flux [W/m2]
    q0_Wm2 = K_SG * sqrt(rho_inf_kgm3 / Rn_m) * (V_inf_ms**3)
    q0_kWm2 = q0_Wm2 / 1000.0d0
    q0_MWm2 = q0_Wm2 / 1.0d6

    ! Radiative Equilibrium Wall Temperature [K]
    Trad_K = (q0_Wm2 / (max(0.1d0, emissivity) * SIGMA))**0.25d0
    Trad_C = Trad_K - 273.15d0

    ! Stagnation Enthalpy [MJ/kg]
    h0_MJkg = (1004.5d0 * T_inf_K + 0.5d0 * (V_inf_ms**2)) / 1.0d6

    ! Rayleigh Pitot / Hypersonic Shock Stagnation Pressure [kPa]
    P02_kPa = (0.92d0 * rho_inf_kgm3 * (V_inf_ms**2)) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — HYPERSONIC AERODYNAMIC HEATING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Altitude / Velocity       = ', altitude_km, ' km / ', V_inf_ms, ' m/s'
    write(*,'(A,F10.2,A,F10.2,A)') 'Free-Stream Mach / Sound  = ', Mach_inf, ' / ', a_sound_ms, ' m/s'
    write(*,'(A,F10.3,A,F10.4,A)') 'Air Density / Pressure    = ', rho_inf_kgm3, ' kg/m3 / ', P_inf_Pa, ' Pa'
    write(*,'(A,F10.3,A,F10.2)')   'Nose Radius / Emissivity  = ', Rn_m, ' m / ', emissivity
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'STAGNATION HEAT FLUX (q0) = ', q0_kWm2, ' kW/m2 (', q0_MWm2, ' MW/m2)'
    write(*,'(A,F10.1,A,F10.1,A)') 'Radiative Wall Temp Trad  = ', Trad_K, ' K (', Trad_C, ' deg C)'
    write(*,'(A,F10.2,A)')  'Stagnation Enthalpy (h0)  = ', h0_MJkg, ' MJ/kg'
    write(*,'(A,F10.2,A)')  'Post-Shock Stagnation P02 = ', P02_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program aerodynamic_heating_reentry