! ============================================================================
! ThermoFluidCalc — Rayleigh Flow Solver
! Reference: Shapiro, Dynamics and Thermodynamics of Compressible Fluid Flow
! ============================================================================
program rayleigh_flow
    implicit none
    
    ! Input variables
    double precision :: M1 ! Inlet Mach number
    double precision :: gamma ! Specific heat ratio (default 1.4)
    double precision :: T01 ! Inlet Stagnation Temp [K] (0 = skip)
    double precision :: P01 ! Inlet Stagnation Pressure [kPa] (0 = skip)
    double precision :: q ! Heat added per unit mass [kJ/kg] (0 = skip)
    
    ! Constants
    double precision, parameter :: R_air = 287.05d0 ! Gas constant [J/kg-K]
    
    ! Intermediate and output variables
    double precision :: Cp, Cp_kJ
    double precision :: T0star, q_max, q_Joules
    double precision :: M2 ! Exit Mach number
    double precision :: T2_T1, T02_T01, P2_P1, P02_P01, rho2_rho1
    double precision :: T1, P1, rho1, V1, T2, P2, rho2, V2
    double precision :: P02, T02 ! exit stagnation properties
    
    double precision :: T_Tstar1, T0_T0star1, P_Pstar1, P0_P0star1, V_Vstar1, rho_rhostar1
    double precision :: T_Tstar2, T0_T0star2, P_Pstar2, P0_P0star2, V_Vstar2, rho_rhostar2
    double precision :: delta_s_over_R
    double precision :: phi_target
    double precision :: M1_orig
    double precision :: P01_eff
    
    integer :: status ! 1=Subsonic unchoked, 2=Supersonic unchoked, 3=Thermally choked
    
    ! Read input from stdin
    read(*,*) M1
    read(*,*) gamma
    read(*,*) T01
    read(*,*) P01
    read(*,*) q
    
    ! Set defaults and bounds
    if (gamma <= 1.0d0) gamma = 1.4d0
    if (M1 <= 0.0d0) M1 = 0.5d0
    if (q < 0.0d0) q = 0.0d0
    
    ! If heat is added, we must have reference stagnation temperature and pressure
    ! If the user did not provide them, we default them
    if (T01 <= 0.0d0) T01 = 300.0d0
    if (P01 <= 0.0d0) P01 = 101.325d0
    
    M1_orig = M1
    Cp = gamma * R_air / (gamma - 1.0d0) ! J/kg-K
    Cp_kJ = Cp / 1000.0d0 ! kJ/kg-K
    
    T0star = T01 / rayleigh_t0_ratio(M1, gamma)
    q_max = Cp_kJ * (T0star - T01)
    
    ! Initialize shock effects for supersonic choking
    P01_eff = P01
    
    ! ── DETERMINE FLOW STATUS AND SOLVE EXIT MACH ────────────
    if (q <= q_max) then
        ! FLOW IS UNCHOKED
        if (M1 < 1.0d0) then
            status = 1 ! Subsonic unchoked
        else
            status = 2 ! Supersonic unchoked
        end if
        
        T02 = T01 + q / Cp_kJ
        phi_target = T02 / T0star
        
        if (M1 < 1.0d0) then
            M2 = solve_mach_from_t0_ratio(phi_target, gamma, .true., M1)
        else
            M2 = solve_mach_from_t0_ratio(phi_target, gamma, .false., M1)
        end if
    else
        ! FLOW IS THERMALLY CHOKED
        status = 3
        M2 = 1.0d0
        T02 = T01 + q / Cp_kJ ! Exit stagnation temperature
        
        ! The new inlet flow must become subsonic and adapt to a lower Mach number.
        ! This holds true for both subsonic and supersonic initial states.
        phi_target = T01 / T02
        M1 = solve_mach_from_t0_ratio(phi_target, gamma, .true., 0.5d0)
        
        ! Stagnation temperature at inlet remains constant
        ! Stagnation pressure at inlet remains constant if subsonic, or suffers shock loss if supersonic
        if (M1_orig > 1.0d0) then
            P01_eff = P01 * shock_p0_ratio(M1_orig, gamma)
        else
            P01_eff = P01
        end if
        
        ! Re-calculate choking conditions for the new Rayleigh line
        T0star = T01 / rayleigh_t0_ratio(M1, gamma)
        q_max = Cp_kJ * (T0star - T01)
    end if
    
    ! ── COMPUTE RAYLEIGH RATIOS AT INLET AND EXIT ────────────
    T_Tstar1 = ((gamma + 1.0d0) * M1 / (1.0d0 + gamma * M1**2))**2
    T0_T0star1 = rayleigh_t0_ratio(M1, gamma)
    P_Pstar1 = (gamma + 1.0d0) / (1.0d0 + gamma * M1**2)
    P0_P0star1 = P_Pstar1 * ((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0))**(gamma / (gamma - 1.0d0))
    V_Vstar1 = (gamma + 1.0d0) * M1**2 / (1.0d0 + gamma * M1**2)
    rho_rhostar1 = 1.0d0 / V_Vstar1
    
    T_Tstar2 = ((gamma + 1.0d0) * M2 / (1.0d0 + gamma * M2**2))**2
    T0_T0star2 = rayleigh_t0_ratio(M2, gamma)
    P_Pstar2 = (gamma + 1.0d0) / (1.0d0 + gamma * M2**2)
    P0_P0star2 = P_Pstar2 * ((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0))**(gamma / (gamma - 1.0d0))
    V_Vstar2 = (gamma + 1.0d0) * M2**2 / (1.0d0 + gamma * M2**2)
    rho_rhostar2 = 1.0d0 / V_Vstar2
    
    ! Property changes across duct
    T2_T1 = T_Tstar2 / T_Tstar1
    T02_T01 = T0_T0star2 / T0_T0star1
    P2_P1 = P_Pstar2 / P_Pstar1
    P02_P01 = (P01_eff / P01) * (P0_P0star2 / P0_P0star1)
    rho2_rho1 = rho_rhostar2 / rho_rhostar1
    delta_s_over_R = (gamma / (gamma - 1.0d0)) * log(T2_T1) - log(P2_P1)
    
    ! ── COMPUTE ACTUAL PHYSICAL PROPERTIES ───────────────────
    T1 = T01 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)
    V1 = M1 * sqrt(gamma * R_air * T1)
    P1 = P01_eff / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)**(gamma / (gamma - 1.0d0))
    rho1 = P1 / (R_air * T1 / 1000.0d0)
    
    T2 = T1 * T2_T1
    V2 = M2 * sqrt(gamma * R_air * T2)
    P2 = P1 * P2_P1
    rho2 = P2 / (R_air * T2 / 1000.0d0)
    P02 = P01_eff * (P0_P0star2 / P0_P0star1)
    
    ! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
    write(*, '(A, I2)') "Status Code = ", status
    select case (status)
        case (1)
            write(*, '(A)') "Status = Subsonic Flow"
        case (2)
            write(*, '(A)') "Status = Supersonic Flow"
        case (3)
            if (M1_orig < 1.0d0) then
                write(*, '(A)') "Status = Thermally Choked Subsonic Flow"
            else
                write(*, '(A)') "Status = Thermally Choked Supersonic Flow (Shocked)"
            end if
    end select
    
    write(*, '(A, F14.6)') "Inlet Mach (M1) = ", M1
    write(*, '(A, F14.6)') "Exit Mach (M2) = ", M2
    write(*, '(A, F14.6)') "Original Inlet Mach = ", M1_orig
    write(*, '(A, F14.6)') "Heat Added (q) = ", q
    write(*, '(A, F14.6)') "Max Heat Added (q_max) = ", q_max
    write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma
    
    write(*, '(A, F14.6)') "Inlet T/Tstar = ", T_Tstar1
    write(*, '(A, F14.6)') "Inlet T0/T0star = ", T0_T0star1
    write(*, '(A, F14.6)') "Inlet P/Pstar = ", P_Pstar1
    write(*, '(A, F14.6)') "Inlet P0/P0star = ", P0_P0star1
    write(*, '(A, F14.6)') "Inlet V/Vstar = ", V_Vstar1
    write(*, '(A, F14.6)') "Inlet rho/rhostar = ", rho_rhostar1
    
    write(*, '(A, F14.6)') "Exit T/Tstar = ", T_Tstar2
    write(*, '(A, F14.6)') "Exit T0/T0star = ", T0_T0star2
    write(*, '(A, F14.6)') "Exit P/Pstar = ", P_Pstar2
    write(*, '(A, F14.6)') "Exit P0/P0star = ", P0_P0star2
    write(*, '(A, F14.6)') "Exit V/Vstar = ", V_Vstar2
    write(*, '(A, F14.6)') "Exit rho/rhostar = ", rho_rhostar2
    
    write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
    write(*, '(A, F14.6)') "Stagnation Temp Ratio (T02/T01) = ", T02_T01
    write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
    write(*, '(A, F14.6)') "Stagnation Pressure Ratio (P02/P01) = ", P02_P01
    write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
    write(*, '(A, F14.6)') "Entropy Change (delta_s/R) = ", delta_s_over_R
    
    write(*, '(A, F14.4)') "Inlet Temperature (T1) = ", T1
    write(*, '(A, F14.4)') "Exit Temperature (T2) = ", T2
    write(*, '(A, F14.2)') "Inlet Velocity (V1) = ", V1
    write(*, '(A, F14.2)') "Exit Velocity (V2) = ", V2
    write(*, '(A, F14.4)') "Inlet Pressure (P1) = ", P1
    write(*, '(A, F14.4)') "Exit Pressure (P2) = ", P2
    write(*, '(A, F14.4)') "Inlet Stagnation Temp (T01) = ", T01
    write(*, '(A, F14.4)') "Exit Stagnation Temp (T02) = T02"
    write(*, '(A, F14.4)') "Inlet Stagnation Pres (P01) = ", P01
    write(*, '(A, F14.4)') "Exit Stagnation Pres (P02) = ", P02
    write(*, '(A, F14.6)') "Inlet Density (rho1) = ", rho1
    write(*, '(A, F14.6)') "Exit Density (rho2) = ", rho2

contains

    ! Stagnation temperature ratio function T0/T0*
    double precision function rayleigh_t0_ratio(M, g)
        double precision, intent(in) :: M, g
        rayleigh_t0_ratio = (2.0d0 * (g + 1.0d0) * M**2 * (1.0d0 + (g - 1.0d0) / 2.0d0 * M**2)) / (1.0d0 + g * M**2)**2
    end function rayleigh_t0_ratio

    ! Stagnation pressure ratio across a normal shock P02/P01
    double precision function shock_p0_ratio(M, g)
        double precision, intent(in) :: M, g
        double precision :: t1, t2
        t1 = (g + 1.0d0) / (2.0d0 * g * M**2 - (g - 1.0d0))
        t2 = ((g + 1.0d0) * M**2) / (2.0d0 + (g - 1.0d0) * M**2)
        shock_p0_ratio = t1**(1.0d0 / (g - 1.0d0)) * t2**(g / (g - 1.0d0))
    end function shock_p0_ratio

    ! Solve Mach number from Stagnation Temperature Ratio using Newton-Raphson
    double precision function solve_mach_from_t0_ratio(target_phi, g, is_sub, init_guess)
        double precision, intent(in) :: target_phi, g, init_guess
        logical, intent(in) :: is_sub
        double precision :: M, M_next, f_val, df_val
        double precision :: num, den, dnum, dden
        integer :: iter
        
        M = init_guess
        if (is_sub .and. M >= 1.0d0) M = 0.5d0
        if (.not. is_sub .and. M <= 1.0d0) M = 2.0d0
        
        ! Safeguard target_phi
        if (target_phi >= 0.99999d0) then
            solve_mach_from_t0_ratio = 1.0d0
            return
        end if
        
        do iter = 1, 100
            ! We solve N(M) - target_phi * D(M) = 0
            ! N(M) = 2(g+1)M^2 + (g^2 - 1)M^4
            ! D(M) = (1 + g M^2)^2
            num = 2.0d0 * (g + 1.0d0) * M**2 + (g**2 - 1.0d0) * M**4
            den = (1.0d0 + g * M**2)**2
            f_val = num - target_phi * den
            
            ! Derivatives:
            ! N'(M) = 4(g+1)M + 4(g^2 - 1)M^3
            ! D'(M) = 4 g M (1 + g M^2)
            dnum = 4.0d0 * (g + 1.0d0) * M + 4.0d0 * (g**2 - 1.0d0) * M**3
            dden = 4.0d0 * g * M * (1.0d0 + g * M**2)
            df_val = dnum - target_phi * dden
            
            if (abs(df_val) < 1.0d-12) exit
            
            M_next = M - f_val / df_val
            
            ! Keep guess in bounds
            if (is_sub) then
                if (M_next <= 0.0d0) M_next = M / 2.0d0
                if (M_next >= 1.0d0) M_next = (M + 1.0d0) / 2.0d0
            else
                if (M_next <= 1.0d0) M_next = (M + 1.0d0) / 2.0d0
            end if
            
            if (abs(M_next - M) < 1.0d-8) then
                M = M_next
                exit
            end if
            
            M = M_next
        end do
        
        solve_mach_from_t0_ratio = M
    end function solve_mach_from_t0_ratio

end program rayleigh_flow