Compressible Fanno Flow with Friction

program compressible_fanno_flow_duct
    implicit none
    integer :: iostat_val, iter
    double precision :: M1_mach, D_duct_mm, L_duct_m, f_fanning, P01_bar, T01_K, gamma_ratio
    double precision :: D_m, f4LD_1, Lmax1_m, f4LD_2, M2_mach, M_low, M_high, M_mid, f4LD_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
    logical :: is_choked

    ! Read inputs
    read(*,*,iostat=iostat_val) M1_mach       ! Inlet Mach Number (e.g. 0.35)
    read(*,*,iostat=iostat_val) D_duct_mm     ! Duct Diameter [mm] (e.g. 100.0)
    read(*,*,iostat=iostat_val) L_duct_m      ! Duct Length [m] (e.g. 25.0)
    read(*,*,iostat=iostat_val) f_fanning     ! Fanning Friction Factor (e.g. 0.005)
    read(*,*,iostat=iostat_val) P01_bar       ! Total Pressure [bar] (e.g. 6.0)
    read(*,*,iostat=iostat_val) T01_K         ! Total Temperature [K] (e.g. 300.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 Fanno flow calculation.'
        stop
    end if

    if (M1_mach <= 0.01d0 .or. D_duct_mm <= 0.0d0 .or. L_duct_m <= 0.0d0 .or. f_fanning <= 0.0d0) then
        write(*,*) 'ERROR: Mach, diameter, length, and friction factor must be positive.'
        stop
    end if

    D_m = D_duct_mm * 1.0d-3

    ! Inlet Fanno Parameter (4fL*/D)1
    f4LD_1 = (1.0d0 - (M1_mach**2)) / (gamma_ratio * (M1_mach**2)) + &
             ((gamma_ratio + 1.0d0) / (2.0d0 * gamma_ratio)) * &
             log(((gamma_ratio + 1.0d0) * (M1_mach**2)) / (2.0d0 + (gamma_ratio - 1.0d0) * (M1_mach**2)))

    Lmax1_m = f4LD_1 * D_m / (4.0d0 * f_fanning)

    if (L_duct_m >= Lmax1_m) then
        is_choked = .true.
        M2_mach = 1.0d0
        f4LD_2 = 0.0d0
    else
        is_choked = .false.
        f4LD_2 = f4LD_1 - (4.0d0 * f_fanning * L_duct_m) / D_m

        ! 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)
            f4LD_mid = (1.0d0 - (M_mid**2)) / (gamma_ratio * (M_mid**2)) + &
                       ((gamma_ratio + 1.0d0) / (2.0d0 * gamma_ratio)) * &
                       log(((gamma_ratio + 1.0d0) * (M_mid**2)) / (2.0d0 + (gamma_ratio - 1.0d0) * (M_mid**2)))
            if (M1_mach < 1.0d0) then
                if (f4LD_mid > f4LD_2) then
                    M_low = M_mid
                else
                    M_high = M_mid
                end if
            else
                if (f4LD_mid < f4LD_2) then
                    M_low = M_mid
                else
                    M_high = M_mid
                end if
            end if
        end do
        M2_mach = M_mid
    end if

    ! Sonic Reference Ratios
    T_Tstar_1 = (gamma_ratio + 1.0d0) / (2.0d0 + (gamma_ratio - 1.0d0) * (M1_mach**2))
    T_Tstar_2 = (gamma_ratio + 1.0d0) / (2.0d0 + (gamma_ratio - 1.0d0) * (M2_mach**2))

    P_Pstar_1 = (1.0d0 / M1_mach) * sqrt(T_Tstar_1)
    P_Pstar_2 = (1.0d0 / M2_mach) * sqrt(T_Tstar_2)

    P0_P0star_1 = (1.0d0 / M1_mach) * ((2.0d0 + (gamma_ratio - 1.0d0) * (M1_mach**2)) / &
                  (gamma_ratio + 1.0d0))**((gamma_ratio + 1.0d0) / (2.0d0 * (gamma_ratio - 1.0d0)))
    P0_P0star_2 = (1.0d0 / M2_mach) * ((2.0d0 + (gamma_ratio - 1.0d0) * (M2_mach**2)) / &
                  (gamma_ratio + 1.0d0))**((gamma_ratio + 1.0d0) / (2.0d0 * (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)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — COMPRESSIBLE FANNO FLOW ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Duct Diam / Length (f)    = ', D_duct_mm, ' mm / ', L_duct_m, ' m (f = ', f_fanning, ')'
    write(*,'(A,F10.2,A,F10.2,A)') 'Inlet Mach / Stagnation P0= ', M1_mach, ' / ', P01_bar, ' bar'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4)')   'EXIT MACH NUMBER (M2)     = ', M2_mach
    write(*,'(A,F10.2,A)') 'Sonic Choking Length Lmax = ', Lmax1_m, ' meters'
    if (is_choked) then
        write(*,'(A)') 'Flow Status               = CHOKED (Sonic Throat M2 = 1.0)'
    else
        write(*,'(A)') 'Flow Status               = UNCHOKED (Friction limited)'
    end if
    write(*,'(A,F10.3,A,F10.3,A)') 'Static Pressure P1 -> P2  = ', P1_bar, ' bar -> ', P2_bar, ' bar'
    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)') '============================================================'

end program compressible_fanno_flow_duct