program nozzle_flow
    implicit none
    
    double precision :: gamma, T0, P0, A_throat, A_exit
    double precision :: P_back, AR, R_gas
    double precision :: pi
    
    ! Isentropic properties
    double precision :: M_sub, M_sup, P_exit_sub, P_exit_sup
    double precision :: T_exit, rho_exit, V_exit, snd_exit
    double precision :: mdot, mdot_choked
    double precision :: P_e, T_e, M_e
    
    ! Shock variables
    double precision :: M_shock_up, M_shock_down
    double precision :: P2_P1, T2_T1, P02_P01
    double precision :: A_shock, AR_shock
    double precision :: P_exit_shock
    
    ! Regime detection
    integer :: regime  ! 1-6
    character(len=60) :: regime_name
    double precision :: P_crit1, P_crit2, P_crit3
    
    ! Chart data
    integer :: i, n_pts
    double precision :: x, AR_x, M_x, P_x, T_x
    double precision :: x_throat, x_shock
    
    ! Bisection variables
    double precision :: M_lo, M_hi, M_mid, AR_mid
    integer :: j
    
    pi = 3.14159265358979d0
    R_gas = 287.0d0  ! J/(kg·K) for air
    n_pts = 60
    
    ! Read inputs
    read(*,*) gamma      ! Specific heat ratio
    read(*,*) T0         ! Stagnation temperature [K]
    read(*,*) P0         ! Stagnation pressure [kPa]
    read(*,*) A_throat   ! Throat area [cm2]
    read(*,*) A_exit     ! Exit area [cm2]
    read(*,*) P_back     ! Back pressure [kPa]
    
    AR = A_exit / A_throat
    
    ! ============================================
    ! FIND SUBSONIC AND SUPERSONIC EXIT MACH
    ! ============================================
    ! Subsonic solution (bisection: M from 0.01 to 1.0)
    M_lo = 0.001d0; M_hi = 1.0d0
    do j = 1, 200
        M_mid = (M_lo + M_hi) / 2.0d0
        AR_mid = area_mach_ratio(M_mid, gamma)
        if (AR_mid > AR) then
            M_lo = M_mid
        else
            M_hi = M_mid
        end if
    end do
    M_sub = M_mid
    
    ! Supersonic solution (bisection: M from 1.0 to 20)
    M_lo = 1.0d0; M_hi = 20.0d0
    do j = 1, 200
        M_mid = (M_lo + M_hi) / 2.0d0
        AR_mid = area_mach_ratio(M_mid, gamma)
        if (AR_mid < AR) then
            M_lo = M_mid
        else
            M_hi = M_mid
        end if
    end do
    M_sup = M_mid
    
    ! Exit pressures for subsonic and supersonic isentropic cases
    P_exit_sub = P0 / isen_P_ratio(M_sub, gamma)
    P_exit_sup = P0 / isen_P_ratio(M_sup, gamma)
    
    ! Normal shock at exit: compute post-shock pressure
    M_shock_down = sqrt(((gamma-1.0d0)*M_sup**2 + 2.0d0) / &
                        (2.0d0*gamma*M_sup**2 - (gamma-1.0d0)))
    P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M_sup**2 - 1.0d0)
    P_exit_shock = P_exit_sup * P2_P1
    
    ! Critical back pressures for regime boundaries
    P_crit1 = P_exit_sub      ! Below: choked flow begins
    P_crit2 = P_exit_shock    ! Below: shock moves into nozzle → exits
    P_crit3 = P_exit_sup      ! Below: overexpanded → design → underexpanded
    
    ! ============================================
    ! DETERMINE OPERATING REGIME
    ! ============================================
    if (P_back >= P0) then
        regime = 1
        regime_name = 'No flow (Pb >= P0)'
        M_e = 0.0d0
        P_e = P_back
        T_e = T0
    else if (P_back >= P_crit1) then
        regime = 2
        regime_name = 'Subsonic throughout (not choked)'
        ! Find actual exit Mach from P_back
        M_e = mach_from_pressure(P_back, P0, gamma)
        P_e = P_back
        T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
    else if (P_back >= P_crit2) then
        regime = 3
        regime_name = 'Normal shock in diverging section'
        ! Need to find shock location
        M_e = M_sub  ! Will be adjusted
        P_e = P_back
        T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
    else if (P_back > P_crit3) then
        regime = 4
        regime_name = 'Overexpanded (oblique shocks at exit)'
        M_e = M_sup
        P_e = P_exit_sup
        T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
    else if (abs(P_back - P_crit3) < 0.001d0) then
        regime = 5
        regime_name = 'Design condition (perfectly expanded)'
        M_e = M_sup
        P_e = P_exit_sup
        T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
    else
        regime = 6
        regime_name = 'Underexpanded (expansion fans at exit)'
        M_e = M_sup
        P_e = P_exit_sup
        T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
    end if
    
    ! Mass flow rate
    if (P_back < P_crit1) then
        ! Choked
        mdot_choked = P0*1000.0d0 * (A_throat*1.0d-4) * &
            sqrt(gamma/(R_gas*T0)) * (2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
        mdot = mdot_choked
    else
        ! Not choked — subsonic throughout
        mdot = P0*1000.0d0 * (A_throat*1.0d-4) * M_e * &
            sqrt(gamma/(R_gas*T0)) / isen_P_ratio(M_e, gamma) * &
            sqrt(1.0d0 + (gamma-1.0d0)/2.0d0*M_e**2)
        mdot_choked = P0*1000.0d0 * (A_throat*1.0d-4) * &
            sqrt(gamma/(R_gas*T0)) * (2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
    end if
    
    ! ============================================
    ! OUTPUT
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CONVERGING-DIVERGING NOZZLE FLOW ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ---'
    write(*,'(A,F12.4)')    '  Specific heat ratio γ   = ', gamma
    write(*,'(A,F12.2,A)')  '  Stagnation temp T₀      = ', T0, ' K'
    write(*,'(A,F12.4,A)')  '  Stagnation pressure P₀  = ', P0, ' kPa'
    write(*,'(A,F12.4,A)')  '  Throat area A*          = ', A_throat, ' cm²'
    write(*,'(A,F12.4,A)')  '  Exit area A_e           = ', A_exit, ' cm²'
    write(*,'(A,F12.4)')    '  Area ratio A_e/A*       = ', AR
    write(*,'(A,F12.4,A)')  '  Back pressure P_b       = ', P_back, ' kPa'
    write(*,*)
    
    write(*,'(A)') '--- OPERATING REGIME ---'
    write(*,'(A,I1,A,A)')   '  Regime ', regime, ': ', trim(regime_name)
    write(*,*)
    
    write(*,'(A)') '--- CRITICAL BACK PRESSURES ---'
    write(*,'(A,F12.4,A)')  '  P_sub (choke onset)     = ', P_crit1, ' kPa'
    write(*,'(A,F12.4,A)')  '  P_shock (shock at exit) = ', P_crit2, ' kPa'
    write(*,'(A,F12.4,A)')  '  P_design (isentropic)   = ', P_crit3, ' kPa'
    write(*,*)
    
    write(*,'(A)') '--- ISENTROPIC EXIT SOLUTIONS ---'
    write(*,'(A,F12.6)')    '  M_exit (subsonic)       = ', M_sub
    write(*,'(A,F12.6)')    '  M_exit (supersonic)     = ', M_sup
    write(*,'(A,F12.4,A)')  '  P_exit (subsonic)       = ', P_exit_sub, ' kPa'
    write(*,'(A,F12.4,A)')  '  P_exit (supersonic)     = ', P_exit_sup, ' kPa'
    write(*,*)
    
    write(*,'(A)') '--- ACTUAL EXIT CONDITIONS ---'
    write(*,'(A,F12.6)')    '  Exit Mach number        = ', M_e
    write(*,'(A,F12.4,A)')  '  Exit pressure           = ', P_e, ' kPa'
    write(*,'(A,F12.2,A)')  '  Exit temperature        = ', T_e, ' K'
    if (M_e > 0.0d0) then
        snd_exit = sqrt(gamma * R_gas * T_e)
        V_exit = M_e * snd_exit
        write(*,'(A,F12.2,A)')  '  Exit velocity           = ', V_exit, ' m/s'
        write(*,'(A,F12.2,A)')  '  Speed of sound (exit)   = ', snd_exit, ' m/s'
    end if
    write(*,*)
    
    write(*,'(A)') '--- MASS FLOW RATE ---'
    write(*,'(A,F12.6,A)')  '  ṁ (actual)              = ', mdot, ' kg/s'
    write(*,'(A,F12.6,A)')  '  ṁ (choked max)          = ', mdot_choked, ' kg/s'
    if (P_back < P_crit1) then
        write(*,'(A)')      '  Status: CHOKED (ṁ = ṁ_max)'
    else
        write(*,'(A,F8.2,A)') '  Status: NOT CHOKED (', (mdot/mdot_choked)*100.0d0, '% of max)'
    end if
    write(*,*)
    
    ! Normal shock info if applicable
    if (regime == 3) then
        write(*,'(A)') '--- NORMAL SHOCK IN NOZZLE ---'
        ! Find shock location by bisection
        call find_shock_location(AR, P_back, P0, gamma, A_throat, A_shock, AR_shock, &
                                  M_shock_up, M_shock_down)
        
        P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M_shock_up**2 - 1.0d0)
        T2_T1 = (1.0d0 + 2.0d0*gamma/(gamma+1.0d0)*(M_shock_up**2-1.0d0)) * &
                ((2.0d0+(gamma-1.0d0)*M_shock_up**2)/((gamma+1.0d0)*M_shock_up**2))
        P02_P01 = (((gamma+1.0d0)*M_shock_up**2)/((gamma-1.0d0)*M_shock_up**2+2.0d0)) &
                  **(gamma/(gamma-1.0d0)) * &
                  (1.0d0/(1.0d0+2.0d0*gamma/(gamma+1.0d0)*(M_shock_up**2-1.0d0))) &
                  **(1.0d0/(gamma-1.0d0))
        
        write(*,'(A,F12.4)')    '  Shock at A/A*           = ', AR_shock
        write(*,'(A,F12.6)')    '  M upstream of shock     = ', M_shock_up
        write(*,'(A,F12.6)')    '  M downstream of shock   = ', M_shock_down
        write(*,'(A,F12.4)')    '  P2/P1 across shock      = ', P2_P1
        write(*,'(A,F12.4)')    '  T2/T1 across shock      = ', T2_T1
        write(*,'(A,F12.6)')    '  P02/P01 (stag. loss)    = ', P02_P01
        write(*,'(A,F8.2,A)')   '  Stag. pressure loss     = ', (1.0d0-P02_P01)*100.0d0, ' %'
        write(*,*)
    end if
    
    ! ============================================
    ! NOZZLE PROFILE TABLE
    ! ============================================
    write(*,'(A)') '--- NOZZLE FLOW PROFILE ---'
    write(*,'(A)') '  x/L       A/A*       M          P/P0       T/T0'
    write(*,'(A)') '  ------    ------     ------     ------     ------'
    
    x_throat = 0.5d0  ! Throat at mid-nozzle
    
    do i = 1, n_pts
        x = dble(i) / dble(n_pts)
        
        ! Area ratio: linear converging then diverging
        if (x <= x_throat) then
            AR_x = AR + (1.0d0 - AR) * (x / x_throat)
        else
            AR_x = 1.0d0 + (AR - 1.0d0) * ((x - x_throat) / (1.0d0 - x_throat))
        end if
        
        if (AR_x < 1.0d0) AR_x = 1.0d0
        
        ! Find Mach based on regime
        if (regime <= 2) then
            M_x = mach_from_area_sub(AR_x, gamma)
        else if (x <= x_throat) then
            M_x = mach_from_area_sub(AR_x, gamma)
        else
            if (regime >= 4) then
                M_x = mach_from_area_sup(AR_x, gamma)
            else
                M_x = mach_from_area_sup(AR_x, gamma)
            end if
        end if
        
        P_x = 1.0d0 / isen_P_ratio(M_x, gamma)
        T_x = 1.0d0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_x**2)
        
        write(*,'(F7.4,3X,F8.4,3X,F8.4,3X,F8.4,3X,F8.4)') x, AR_x, M_x, P_x, T_x
    end do
    write(*,*)
    
    ! ============================================
    ! CHART DATA
    ! ============================================
    write(*,'(A)') '--- CHART_DATA ---'
    do i = 1, n_pts
        x = dble(i) / dble(n_pts)
        
        if (x <= x_throat) then
            AR_x = AR + (1.0d0 - AR) * (x / x_throat)
        else
            AR_x = 1.0d0 + (AR - 1.0d0) * ((x - x_throat) / (1.0d0 - x_throat))
        end if
        if (AR_x < 1.0d0) AR_x = 1.0d0
        
        if (regime <= 2) then
            M_x = mach_from_area_sub(AR_x, gamma)
        else if (x <= x_throat) then
            M_x = mach_from_area_sub(AR_x, gamma)
        else
            if (regime >= 4) then
                M_x = mach_from_area_sup(AR_x, gamma)
            else
                M_x = mach_from_area_sup(AR_x, gamma)
            end if
        end if
        
        P_x = 1.0d0 / isen_P_ratio(M_x, gamma)
        T_x = 1.0d0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_x**2)
        
        write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,F10.6)') &
            x, ',', AR_x, ',', M_x, ',', P_x, ',', T_x
    end do
    write(*,'(A)') '--- END_CHART_DATA ---'
    
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  A/A* = (1/M)[(2/(g+1))(1+(g-1)/2 M^2)]^((g+1)/(2(g-1)))'
    write(*,'(A)') '  T0/T = 1 + (g-1)/2 x M^2'
    write(*,'(A)') '  P0/P = (T0/T)^(g/(g-1))'
    write(*,'(A)') '  mdot_choked = P0*A*sqrt(g/(R*T0)) x (2/(g+1))^((g+1)/(2(g-1)))'
    write(*,'(A)') '  Normal shock: M2^2 = [(g-1)M1^2+2]/[2g*M1^2-(g-1)]'
    
contains

    ! Area-Mach ratio function
    function area_mach_ratio(M, g) result(ar)
        double precision, intent(in) :: M, g
        double precision :: ar
        ar = (1.0d0/M) * ((2.0d0/(g+1.0d0)) * &
             (1.0d0 + (g-1.0d0)/2.0d0 * M**2))**((g+1.0d0)/(2.0d0*(g-1.0d0)))
    end function
    
    ! Isentropic pressure ratio P0/P
    function isen_P_ratio(M, g) result(pr)
        double precision, intent(in) :: M, g
        double precision :: pr
        pr = (1.0d0 + (g-1.0d0)/2.0d0 * M**2)**(g/(g-1.0d0))
    end function
    
    ! Find Mach from pressure ratio (subsonic)
    function mach_from_pressure(Pb, P0_val, g) result(M_out)
        double precision, intent(in) :: Pb, P0_val, g
        double precision :: M_out
        M_out = sqrt(2.0d0/(g-1.0d0) * ((P0_val/Pb)**((g-1.0d0)/g) - 1.0d0))
    end function
    
    ! Find subsonic Mach from area ratio (bisection)
    function mach_from_area_sub(ar_target, g) result(M_out)
        double precision, intent(in) :: ar_target, g
        double precision :: M_out, ml, mh, mm, arm
        integer :: jj
        ml = 0.001d0; mh = 1.0d0
        do jj = 1, 200
            mm = (ml + mh) / 2.0d0
            arm = area_mach_ratio(mm, g)
            if (arm > ar_target) then
                ml = mm
            else
                mh = mm
            end if
        end do
        M_out = mm
    end function
    
    ! Find supersonic Mach from area ratio (bisection)
    function mach_from_area_sup(ar_target, g) result(M_out)
        double precision, intent(in) :: ar_target, g
        double precision :: M_out, ml, mh, mm, arm
        integer :: jj
        ml = 1.0d0; mh = 20.0d0
        do jj = 1, 200
            mm = (ml + mh) / 2.0d0
            arm = area_mach_ratio(mm, g)
            if (arm < ar_target) then
                ml = mm
            else
                mh = mm
            end if
        end do
        M_out = mm
    end function
    
    ! Find shock location in diverging section
    subroutine find_shock_location(AR_exit, Pb, P0_val, g, A_t_val, A_sh, AR_sh, M_up, M_dn)
        double precision, intent(in) :: AR_exit, Pb, P0_val, g, A_t_val
        double precision, intent(out) :: A_sh, AR_sh, M_up, M_dn
        double precision :: ar_lo, ar_hi, ar_mid
        double precision :: M1, M2, P02_01, P0_2, AR_downstream
        double precision :: P_exit_test, M_exit_test
        integer :: kk
        
        ar_lo = 1.0d0
        ar_hi = AR_exit
        
        do kk = 1, 200
            ar_mid = (ar_lo + ar_hi) / 2.0d0
            
            ! Supersonic Mach at shock location
            M1 = mach_from_area_sup(ar_mid, g)
            
            ! Post-shock Mach
            M2 = sqrt(((g-1.0d0)*M1**2 + 2.0d0) / (2.0d0*g*M1**2 - (g-1.0d0)))
            
            ! Stagnation pressure ratio across shock
            P02_01 = (((g+1.0d0)*M1**2)/((g-1.0d0)*M1**2+2.0d0))**(g/(g-1.0d0)) * &
                      (1.0d0/(1.0d0+2.0d0*g/(g+1.0d0)*(M1**2-1.0d0)))**(1.0d0/(g-1.0d0))
            
            P0_2 = P0_val * P02_01
            
            ! Post-shock A/A* for M2
            AR_downstream = area_mach_ratio(M2, g)
            
            ! Exit ratio: A_exit/A*_new = (AR_exit/ar_mid) * AR_downstream
            M_exit_test = mach_from_area_sub((AR_exit/ar_mid) * AR_downstream, g)
            P_exit_test = P0_2 / isen_P_ratio(M_exit_test, g)
            
            if (P_exit_test > Pb) then
                ar_lo = ar_mid
            else
                ar_hi = ar_mid
            end if
        end do
        
        AR_sh = ar_mid
        A_sh = ar_mid * A_t_val
        M_up = mach_from_area_sup(ar_mid, g)
        M_dn = sqrt(((g-1.0d0)*M_up**2 + 2.0d0) / (2.0d0*g*M_up**2 - (g-1.0d0)))
    end subroutine

end program nozzle_flow