program compressible_rayleigh_flow_heat
    implicit none
    integer :: iostat_val, iter
    double precision :: M1_mach, T01_K, q_heat_kJkg, cp_kJkgK, P01_bar, gamma_ratio
    double precision :: R_gas, T0_T0star_1, T0_star_K, q_max_kJkg, T02_K, T0_T0star_2
    double precision :: M2_mach, M_low, M_high, M_mid, T0_mid
    double precision :: T_Tstar_1, T_Tstar_2, P_Pstar_1, P_Pstar_2, P0_P0star_1, P0_P0star_2
    double precision :: P1_bar, P2_bar, P02_bar, T1_K, T2_K, delta_s_JkgK
    logical :: is_choked

    ! Read inputs
    read(*,*,iostat=iostat_val) M1_mach       ! Inlet Mach Number (e.g. 0.30)
    read(*,*,iostat=iostat_val) T01_K         ! Inlet Total Temp [K] (e.g. 300.0)
    read(*,*,iostat=iostat_val) q_heat_kJkg   ! Heat Added [kJ/kg] (e.g. 400.0)
    read(*,*,iostat=iostat_val) cp_kJkgK      ! Specific Heat [kJ/kg.K] (e.g. 1.005)
    read(*,*,iostat=iostat_val) P01_bar       ! Inlet Total Pressure [bar] (e.g. 5.0)
    read(*,*,iostat=iostat_val) gamma_ratio   ! Specific Heat Ratio (e.g. 1.40)

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

    if (M1_mach <= 0.01d0 .or. T01_K <= 0.0d0 .or. cp_kJkgK <= 0.0d0 .or. P01_bar <= 0.0d0) then
        write(*,*) 'ERROR: Mach, temperature, cp, and pressure must be positive.'
        stop
    end if

    R_gas = cp_kJkgK * (1.0d0 - 1.0d0 / gamma_ratio) * 1000.0d0 ! J/kg.K

    ! Inlet Total Temperature Ratio
    T0_T0star_1 = (2.0d0 * (gamma_ratio + 1.0d0) * (M1_mach**2) * &
                  (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))) / &
                  ((1.0d0 + gamma_ratio * (M1_mach**2))**2)

    T0_star_K = T01_K / T0_T0star_1
    q_max_kJkg = cp_kJkgK * (T0_star_K - T01_K)

    if (q_heat_kJkg >= q_max_kJkg) then
        is_choked = .true.
        M2_mach = 1.0d0
        T02_K = T0_star_K
        T0_T0star_2 = 1.0d0
    else
        is_choked = .false.
        T02_K = T01_K + (q_heat_kJkg / cp_kJkgK)
        T0_T0star_2 = T02_K / T0_star_K

        ! Solve for M2 via Bisection
        if (M1_mach < 1.0d0) then
            M_low = M1_mach
            M_high = 0.9999d0
        else
            M_low = 1.0001d0
            M_high = M1_mach
        end if

        do iter = 1, 100
            M_mid = 0.5d0 * (M_low + M_high)
            T0_mid = (2.0d0 * (gamma_ratio + 1.0d0) * (M_mid**2) * &
                     (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M_mid**2))) / &
                     ((1.0d0 + gamma_ratio * (M_mid**2))**2)
            if (M1_mach < 1.0d0) then
                if (T0_mid < T0_T0star_2) then
                    M_low = M_mid
                else
                    M_high = M_mid
                end if
            else
                if (T0_mid < T0_T0star_2) then
                    M_high = M_mid
                else
                    M_low = M_mid
                end if
            end if
        end do
        M2_mach = M_mid
    end if

    ! Static and Stagnation Ratios
    T_Tstar_1 = (M1_mach**2) * ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2)))**2
    T_Tstar_2 = (M2_mach**2) * ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2)))**2

    P_Pstar_1 = (gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2))
    P_Pstar_2 = (gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2))

    P0_P0star_1 = ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2))) * &
                  (((2.0d0 + (gamma_ratio - 1.0d0) * (M1_mach**2)) / (gamma_ratio + 1.0d0))**(gamma_ratio / (gamma_ratio - 1.0d0)))
    P0_P0star_2 = ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2))) * &
                  (((2.0d0 + (gamma_ratio - 1.0d0) * (M2_mach**2)) / (gamma_ratio + 1.0d0))**(gamma_ratio / (gamma_ratio - 1.0d0)))

    P1_bar = P01_bar * (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))**(-gamma_ratio / (gamma_ratio - 1.0d0))
    P2_bar = P1_bar * (P_Pstar_2 / P_Pstar_1)
    P02_bar = P01_bar * (P0_P0star_2 / P0_P0star_1)
    T1_K = T01_K / (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))
    T2_K = T1_K * (T_Tstar_2 / T_Tstar_1)

    delta_s_JkgK = (cp_kJkgK * 1000.0d0) * log(T2_K / T1_K) - R_gas * log(P2_bar / P1_bar)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — COMPRESSIBLE RAYLEIGH HEAT FLOW ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Inlet Mach / Stagnation T01 = ', M1_mach, ' / ', T01_K, ' K'
    write(*,'(A,F10.2,A,F10.2,A)') 'Heat Added / Max Choking q  = ', q_heat_kJkg, ' kJ/kg / ', q_max_kJkg, ' kJ/kg'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4)')   'EXIT MACH NUMBER (M2)       = ', M2_mach
    if (is_choked) then
        write(*,'(A)') 'Flow Status                 = THERMALLY CHOKED (M2 = 1.0)'
    else
        write(*,'(A)') 'Flow Status                 = UNCHOKED HEAT ADDITION'
    end if
    write(*,'(A,F10.1,A,F10.1,A)') 'Total Temperature T01 -> T02= ', T01_K, ' K -> ', T02_K, ' K'
    write(*,'(A,F10.3,A,F10.3,A)') 'Total Pressure P01 -> P02   = ', P01_bar, ' bar -> ', P02_bar, ' bar'
    write(*,'(A,F10.1,A,F10.1,A)') 'Static Temp T1 -> T2        = ', T1_K, ' K -> ', T2_K, ' K'
    write(*,'(A,F10.2,A)')  'Entropy Generation Delta s  = ', delta_s_JkgK, ' J/(kg.K)'
    write(*,'(A)') '============================================================'

end program compressible_rayleigh_flow_heat