! ============================================================================
! ThermoFluidCalc — Fanno Flow Solver
! Reference: Shapiro, Dynamics and Thermodynamics of Compressible Fluid Flow
! ============================================================================
program fanno_flow
    implicit none
    
    ! Input variables
    double precision :: M1 ! Inlet Mach number
    double precision :: f ! Darcy friction factor
    double precision :: D ! Pipe diameter [mm]
    double precision :: L ! Pipe length [m]
    double precision :: gamma ! Specific heat ratio (default 1.4)
    double precision :: T0 ! Inlet Stagnation Temp [K] (0 = skip)
    double precision :: P0 ! Inlet Stagnation Pressure [kPa] (0 = skip)
    
    ! Constants
    double precision, parameter :: R_air = 287.0d0 ! Gas constant [J/kg-K]
    
    ! Intermediate and output variables
    double precision :: D_m ! Diameter in meters
    double precision :: fLstar_over_D_inlet, fL_over_D
    double precision :: Lstar ! Choking length [m]
    double precision :: Lstar_max_shock
    double precision :: M2 ! Exit Mach number
    double precision :: T2_T1, P2_P1, rho2_rho1, P02_P01
    double precision :: T1, P1, rho1, V1, T2, P2, rho2, V2
    double precision :: P02, T02 ! exit stagnation properties
    double precision :: T_Tstar1, P_Pstar1, rho_rhostar1, V_Vstar1, P0_P0star1
    double precision :: T_Tstar2, P_Pstar2, rho_rhostar2, V_Vstar2, P0_P0star2
    double precision :: delta_s_over_R
    double precision :: M_shock1, M_shock2, x_shock
    integer :: status ! 1=Subsonic, 2=Supersonic, 3=Choked Subsonic, 4=Choked Shock in Duct, 5=Choked Shock Pushed Upstream
    
    ! Read input from stdin
    read(*,*) M1
    read(*,*) f
    read(*,*) D
    read(*,*) L
    read(*,*) gamma
    read(*,*) T0
    read(*,*) P0
    
    ! Set defaults and bounds
    if (gamma <= 1.0d0) gamma = 1.4d0
    if (f <= 0.0d0) f = 0.02d0
    if (D <= 0.0d0) D = 50.0d0
    if (L < 0.0d0) L = 0.0d0
    if (M1 <= 0.0d0) M1 = 0.5d0
    
    D_m = D / 1000.0d0
    fLstar_over_D_inlet = fanno_param(M1, gamma)
    Lstar = fLstar_over_D_inlet * D_m / f
    fL_over_D = f * L / D_m
    
    ! Initialize shock variables
    M_shock1 = 0.0d0
    M_shock2 = 0.0d0
    x_shock = 0.0d0
    
    ! ── DETERMINE FLOW STATUS AND SOLVE EXIT MACH ────────────
    if (M1 < 1.0d0) then
        ! SUBSONIC INLET FLOW
        if (L <= Lstar) then
            status = 1 ! Subsonic unchoked
            M2 = solve_mach_from_fanno(fLstar_over_D_inlet - fL_over_D, gamma, M1, 1.0d0)
        else
            status = 3 ! Choked Subsonic (Inlet Mach must decrease)
            ! Solve new inlet Mach such that Lstar = L
            M1 = solve_mach_from_fanno(fL_over_D, gamma, 0.001d0, 1.0d0)
            fLstar_over_D_inlet = fanno_param(M1, gamma)
            Lstar = L
            M2 = 1.0d0
        end if
    else
        ! SUPERSONIC INLET FLOW
        if (L <= Lstar) then
            status = 2 ! Supersonic unchoked
            M2 = solve_mach_from_fanno(fLstar_over_D_inlet - fL_over_D, gamma, 1.0d0, M1)
        else
            ! Check if shock can fit inside the duct
            ! Shock at inlet means flow becomes subsonic at inlet.
            ! Post-shock Mach of M1:
            M_shock2 = sqrt((2.0d0 + (gamma - 1.0d0) * M1**2) / (2.0d0 * gamma * M1**2 - (gamma - 1.0d0)))
            Lstar_max_shock = Lstar + fanno_param(M_shock2, gamma) * D_m / f
            
            if (L <= Lstar_max_shock) then
                status = 4 ! Shock inside the duct
                M2 = 1.0d0
                ! Solve for shock Mach number M_shock1 in range [1.0, M1]
                call solve_shock_position(M1, f, D_m, L, gamma, M_shock1, M_shock2, x_shock)
            else
                status = 5 ! Shock pushed back upstream
                ! Flow is completely subsonic in the duct, choked at exit.
                ! Inlet Mach is reduced to subsonic value
                M1 = solve_mach_from_fanno(fL_over_D, gamma, 0.001d0, 1.0d0)
                fLstar_over_D_inlet = fanno_param(M1, gamma)
                Lstar = L
                M2 = 1.0d0
            end if
        end if
    end if
    
    ! ── COMPUTE FANNO RATIOS AT INLET AND EXIT ───────────────
    T_Tstar1 = (gamma + 1.0d0) / (2.0d0 + (gamma - 1.0d0) * M1**2)
    P_Pstar1 = (1.0d0 / M1) * sqrt(T_Tstar1)
    rho_rhostar1 = (1.0d0 / M1) * sqrt((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0))
    V_Vstar1 = 1.0d0 / rho_rhostar1
    P0_P0star1 = (1.0d0 / M1) * ((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0))**((gamma + 1.0d0) / (2.0d0 * (gamma - 1.0d0)))
    
    T_Tstar2 = (gamma + 1.0d0) / (2.0d0 + (gamma - 1.0d0) * M2**2)
    P_Pstar2 = (1.0d0 / M2) * sqrt(T_Tstar2)
    rho_rhostar2 = (1.0d0 / M2) * sqrt((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0))
    V_Vstar2 = 1.0d0 / rho_rhostar2
    P0_P0star2 = (1.0d0 / M2) * ((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0))**((gamma + 1.0d0) / (2.0d0 * (gamma - 1.0d0)))
    
    ! Property changes across duct
    T2_T1 = T_Tstar2 / T_Tstar1
    P2_P1 = P_Pstar2 / P_Pstar1
    rho2_rho1 = rho_rhostar2 / rho_rhostar1
    P02_P01 = P0_P0star2 / P0_P0star1
    delta_s_over_R = log(P0_P0star1 / P0_P0star2)
    
    ! ── COMPUTE OPTIONAL ACTUAL PROPERTIES ───────────────────
    T1 = 0.0d0; P1 = 0.0d0; rho1 = 0.0d0; V1 = 0.0d0
    T2 = 0.0d0; P2 = 0.0d0; rho2 = 0.0d0; V2 = 0.0d0
    P02 = 0.0d0; T02 = 0.0d0
    
    if (T0 > 0.0d0) then
        T1 = T0 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)
        V1 = M1 * sqrt(gamma * R_air * T1)
        T2 = T1 * T2_T1
        V2 = M2 * sqrt(gamma * R_air * T2)
        T02 = T0 ! Stagnation temperature is constant in adiabatic Fanno flow
    end if
    
    if (P0 > 0.0d0) then
        P1 = P0 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)**(gamma / (gamma - 1.0d0))
        P2 = P1 * P2_P1
        P02 = P0 * P02_P01
        
        if (T0 > 0.0d0) then
            rho1 = P1 / (R_air * T1 / 1000.0d0) ! P1 is in kPa, convert R to kJ/kg-K
            rho2 = P2 / (R_air * T2 / 1000.0d0)
        end if
    end if
    
    ! ── 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)
            write(*, '(A)') "Status = Choked Subsonic Flow (Inlet Mach Reduced)"
        case (4)
            write(*, '(A)') "Status = Choked Supersonic Flow (Normal Shock in Duct)"
        case (5)
            write(*, '(A)') "Status = Choked Supersonic Flow (Shock Pushed Upstream)"
    end select
    
    write(*, '(A, F14.6)') "Inlet Mach (M1) = ", M1
    write(*, '(A, F14.6)') "Exit Mach (M2) = ", M2
    write(*, '(A, F14.6)') "Friction Factor (f) = ", f
    write(*, '(A, F14.6)') "Pipe Diameter (D) = ", D
    write(*, '(A, F14.6)') "Pipe Length (L) = ", L
    write(*, '(A, F14.6)') "Choking Length (Lstar) = ", Lstar
    write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma
    
    write(*, '(A, F14.6)') "Inlet T/Tstar = ", T_Tstar1
    write(*, '(A, F14.6)') "Inlet P/Pstar = ", P_Pstar1
    write(*, '(A, F14.6)') "Inlet rho/rhostar = ", rho_rhostar1
    write(*, '(A, F14.6)') "Inlet V/Vstar = ", V_Vstar1
    write(*, '(A, F14.6)') "Inlet P0/P0star = ", P0_P0star1
    
    write(*, '(A, F14.6)') "Exit T/Tstar = ", T_Tstar2
    write(*, '(A, F14.6)') "Exit P/Pstar = ", P_Pstar2
    write(*, '(A, F14.6)') "Exit rho/rhostar = ", rho_rhostar2
    write(*, '(A, F14.6)') "Exit V/Vstar = ", V_Vstar2
    write(*, '(A, F14.6)') "Exit P0/P0star = ", P0_P0star2
    
    write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
    write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
    write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
    write(*, '(A, F14.6)') "Stagnation Pressure Ratio (P02/P01) = ", P02_P01
    write(*, '(A, F14.6)') "Entropy Change (delta_s/R) = ", delta_s_over_R
    
    if (status == 4) then
        write(*, '(A, F14.6)') "Shock Pre-Mach (Ms1) = ", M_shock1
        write(*, '(A, F14.6)') "Shock Post-Mach (Ms2) = ", M_shock2
        write(*, '(A, F14.6)') "Shock Position (x_shock) = ", x_shock
    end if
    
    if (T0 > 0.0d0) then
        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 Stagnation Temp (T01) = ", T0
        write(*, '(A, F14.4)') "Exit Stagnation Temp (T02) = ", T02
    end if
    
    if (P0 > 0.0d0) then
        write(*, '(A, F14.4)') "Inlet Pressure (P1) = ", P1
        write(*, '(A, F14.4)') "Exit Pressure (P2) = ", P2
        write(*, '(A, F14.4)') "Inlet Stagnation Pres (P01) = ", P0
        write(*, '(A, F14.4)') "Exit Stagnation Pres (P02) = ", P02
        if (T0 > 0.0d0) then
            write(*, '(A, F14.6)') "Inlet Density (rho1) = ", rho1
            write(*, '(A, F14.6)') "Exit Density (rho2) = ", rho2
        end if
    end if

contains

    ! Function to compute fL*/D for a given Mach number
    double precision function fanno_param(M, g)
        double precision, intent(in) :: M, g
        double precision :: Term1, Term2
        if (M <= 0.0d0) then
            fanno_param = 1.0d10
        else
            Term1 = (1.0d0 - M**2) / (g * M**2)
            Term2 = (g + 1.0d0) / (2.0d0 * g) * log(((g + 1.0d0) * M**2) / (2.0d0 + (g - 1.0d0) * M**2))
            fanno_param = Term1 + Term2
        end if
    end function fanno_param
    
    ! Function derivative of Fanno parameter dfL*/D / dM
    double precision function fanno_df(M, g)
        double precision, intent(in) :: M, g
        fanno_df = -2.0d0 * (1.0d0 - M**2) / (g * M**3 * (1.0d0 + (g - 1.0d0)/2.0d0 * M**2))
    end function fanno_df
    
    ! Newton-Raphson solver to find Mach number from a target Fanno parameter value
    double precision function solve_mach_from_fanno(target_val, g, lower_bound, upper_bound)
        double precision, intent(in) :: target_val, g, lower_bound, upper_bound
        double precision :: M, diff, f_val, df_val
        integer :: i
        
        if (target_val <= 0.0d0) then
            solve_mach_from_fanno = 1.0d0
            return
        end if
        
        ! Initial guess
        M = 0.5d0 * (lower_bound + upper_bound)
        if (lower_bound < 1.0d0 .and. upper_bound <= 1.0d0) then
            ! Subsonic guess
            M = 0.5d0
        else if (lower_bound >= 1.0d0) then
            ! Supersonic guess
            M = 2.0d0
        end if
        
        do i = 1, 100
            f_val = fanno_param(M, g) - target_val
            df_val = fanno_df(M, g)
            diff = f_val / df_val
            M = M - diff
            
            ! Keep M in bounds
            if (M < lower_bound) M = lower_bound + 1.0d-5
            if (M > upper_bound) M = upper_bound - 1.0d-5
            
            if (abs(diff) < 1.0d-12) exit
        end do
        
        solve_mach_from_fanno = M
    end function solve_mach_from_fanno
    
    ! Subroutine to solve for shock Mach number and position
    subroutine solve_shock_position(M1, f_fac, D_val, L_val, g, Ms1, Ms2, xs)
        double precision, intent(in) :: M1, f_fac, D_val, L_val, g
        double precision, intent(out) :: Ms1, Ms2, xs
        double precision :: low_M, high_M, mid_M, fL_D, test_val
        integer :: iter
        
        fL_D = f_fac * L_val / D_val
        
        ! We use bisection solver to find Ms1 in [1.0, M1]
        low_M = 1.0d0
        high_M = M1
        
        do iter = 1, 100
            mid_M = 0.5d0 * (low_M + high_M)
            
            ! Post-shock Mach number of mid_M
            Ms2 = sqrt((2.0d0 + (g - 1.0d0) * mid_M**2) / (2.0d0 * g * mid_M**2 - (g - 1.0d0)))
            
            ! Evaluate constraint: L_shock_Lstar = L*(M1) - L*(Ms1) + L*(Ms2)
            test_val = fanno_param(M1, g) - fanno_param(mid_M, g) + fanno_param(Ms2, g)
            
            if (test_val < fL_D) then
                ! Need higher mid_M (less shock distance)
                low_M = mid_M
            else
                ! Need lower mid_M
                high_M = mid_M
            end if
            
            if (abs(high_M - low_M) < 1.0d-12) exit
        end do
        
        Ms1 = mid_M
        Ms2 = sqrt((2.0d0 + (g - 1.0d0) * Ms1**2) / (2.0d0 * g * Ms1**2 - (g - 1.0d0)))
        xs = (fanno_param(M1, g) - fanno_param(Ms1, g)) * D_val / f_fac
    end subroutine solve_shock_position

end program fanno_flow
