program supersonic_converging_diverging_nozzle
    implicit none
    integer :: iostat_val, iter
    double precision :: Dt_throat_mm, Ae_At_ratio, L_div_mm, P0_bar, T0_K, Pb_bar, gamma_ratio
    double precision :: At_m2, Ae_m2, R_gas, mdot_choked_kgs, Me_design, Pe_design_bar
    double precision :: M_sub_exit, P_crit1_bar, P_crit2_bar, Me_actual, Pe_actual_bar
    double precision :: x_shock_pct, Ms_pre, M_low, M_high, M_mid, ar_mid
    double precision :: Ve_actual_ms, thrust_N, CF_thrust
    character(len=40) :: nozzle_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Dt_throat_mm    ! Throat Diameter [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) Ae_At_ratio     ! Exit Area Ratio Ae/At (e.g. 2.50)
    read(*,*,iostat=iostat_val) L_div_mm        ! Diverging Length [mm] (e.g. 180.0)
    read(*,*,iostat=iostat_val) P0_bar          ! Stagnation Pressure [bar] (e.g. 10.0)
    read(*,*,iostat=iostat_val) T0_K            ! Stagnation Temp [K] (e.g. 500.0)
    read(*,*,iostat=iostat_val) Pb_bar          ! Back Pressure [bar] (e.g. 1.013)
    read(*,*,iostat=iostat_val) gamma_ratio     ! Specific Heat Ratio (e.g. 1.40)

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

    if (Dt_throat_mm <= 0.0d0 .or. Ae_At_ratio <= 1.0d0 .or. P0_bar <= 0.0d0 .or. Pb_bar <= 0.0d0) then
        write(*,*) 'ERROR: All geometric and pressure parameters must be positive, Ae/At > 1.'
        stop
    end if

    R_gas = 287.05d0
    At_m2 = (PI / 4.0d0) * ((Dt_throat_mm * 1.0d-3)**2)
    Ae_m2 = At_m2 * Ae_At_ratio

    ! Choked Mass Flow Rate [kg/s]
    mdot_choked_kgs = (P0_bar * 1.0d5 * At_m2 / sqrt(T0_K)) * sqrt(gamma_ratio / R_gas) * &
                      ((2.0d0 / (gamma_ratio + 1.0d0))**((gamma_ratio + 1.0d0) / (2.0d0 * (gamma_ratio - 1.0d0))))

    ! Solve for Design Supersonic Exit Mach Me_design via Bisection
    M_low = 1.001d0
    M_high = 8.0d0
    do iter = 1, 80
        M_mid = 0.5d0 * (M_low + M_high)
        ar_mid = (1.0d0 / M_mid) * (((2.0d0 + (gamma_ratio - 1.0d0) * (M_mid**2)) / &
                 (gamma_ratio + 1.0d0))**((gamma_ratio + 1.0d0) / (2.0d0 * (gamma_ratio - 1.0d0))))
        if (ar_mid < Ae_At_ratio) then
            M_low = M_mid
        else
            M_high = M_mid
        end if
    end do
    Me_design = M_mid
    Pe_design_bar = P0_bar * (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (Me_design**2))**(-gamma_ratio / (gamma_ratio - 1.0d0))

    ! Solve for Subsonic Exit Mach M_sub_exit
    M_low = 0.001d0
    M_high = 0.999d0
    do iter = 1, 80
        M_mid = 0.5d0 * (M_low + M_high)
        ar_mid = (1.0d0 / M_mid) * (((2.0d0 + (gamma_ratio - 1.0d0) * (M_mid**2)) / &
                 (gamma_ratio + 1.0d0))**((gamma_ratio + 1.0d0) / (2.0d0 * (gamma_ratio - 1.0d0))))
        if (ar_mid > Ae_At_ratio) then
            M_low = M_mid
        else
            M_high = M_mid
        end if
    end do
    M_sub_exit = M_mid
    P_crit1_bar = P0_bar * (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M_sub_exit**2))**(-gamma_ratio / (gamma_ratio - 1.0d0))

    ! Normal Shock at Exit Plane Pressure P_crit2
    P_crit2_bar = Pe_design_bar * (2.0d0 * gamma_ratio * (Me_design**2) - (gamma_ratio - 1.0d0)) / (gamma_ratio + 1.0d0)

    ! Regime Classification & Shock Position Determination
    if (Pb_bar >= P_crit1_bar) then
        nozzle_regime = 'Subsonic Venturi Flow (Unchoked Throat)'
        x_shock_pct = 0.0d0
        Me_actual = M_sub_exit
        Pe_actual_bar = Pb_bar
    else if (Pb_bar > P_crit2_bar) then
        nozzle_regime = 'Normal Shock in Diverging Section'
        ! Shock position interpolated between throat (0%) and exit (100%)
        x_shock_pct = ((P_crit1_bar - Pb_bar) / (P_crit1_bar - P_crit2_bar)) * 100.0d0
        Me_actual = 0.65d0
        Pe_actual_bar = Pb_bar
    else if (abs(Pb_bar - P_crit2_bar) < 0.02d0 * P_crit2_bar) then
        nozzle_regime = 'Normal Shock at Exit Plane'
        x_shock_pct = 100.0d0
        Me_actual = sqrt((1.0d0 + 0.5d0*(gamma_ratio-1.0d0)*(Me_design**2)) / (gamma_ratio*(Me_design**2) - 0.5d0*(gamma_ratio-1.0d0)))
        Pe_actual_bar = Pb_bar
    else if (Pb_bar > Pe_design_bar) then
        nozzle_regime = 'Over-Expanded (Oblique Exit Shocks)'
        x_shock_pct = 100.0d0
        Me_actual = Me_design
        Pe_actual_bar = Pe_design_bar
    else if (abs(Pb_bar - Pe_design_bar) < 0.02d0 * Pe_design_bar) then
        nozzle_regime = 'Perfectly Expanded (Design Condition)'
        x_shock_pct = 100.0d0
        Me_actual = Me_design
        Pe_actual_bar = Pe_design_bar
    else
        nozzle_regime = 'Under-Expanded (Exit Expansion Fan)'
        x_shock_pct = 100.0d0
        Me_actual = Me_design
        Pe_actual_bar = Pe_design_bar
    end if

    ! Exit velocity & Thrust
    Ve_actual_ms = Me_actual * sqrt(gamma_ratio * R_gas * (T0_K / (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (Me_actual**2))))
    thrust_N = mdot_choked_kgs * Ve_actual_ms + (Pe_actual_bar - Pb_bar) * 1.0d5 * Ae_m2
    CF_thrust = thrust_N / (P0_bar * 1.0d5 * At_m2)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERSONIC DE LAVAL NOZZLE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2)') 'Throat Diam / Area Ratio = ', Dt_throat_mm, ' mm / ', Ae_At_ratio
    write(*,'(A,F10.2,A,F10.3,A)')'Stagnation P0 / Back Pb  = ', P0_bar, ' bar / ', Pb_bar, ' bar'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,A)')       'Flow Regime               = ', trim(nozzle_regime)
    write(*,'(A,F10.3)')   'Design Supersonic Mach Me = ', Me_design
    write(*,'(A,F10.3,A)') 'Design Exit Pressure Pe   = ', Pe_design_bar, ' bar'
    if (x_shock_pct > 0.0d0 .and. x_shock_pct < 100.0d0) then
        write(*,'(A,F10.1,A,F10.1,A)') 'Normal Shock Location     = ', x_shock_pct, ' % (x = ', x_shock_pct*L_div_mm/100.0d0, ' mm)'
    end if
    write(*,'(A,F10.3,A)') 'Mass Flow Rate (Choked)   = ', mdot_choked_kgs, ' kg/s'
    write(*,'(A,F10.1,A,F10.3)') 'Thrust Force / Coeff CF   = ', thrust_N, ' N / ', CF_thrust
    write(*,'(A)') '============================================================'

end program supersonic_converging_diverging_nozzle