program rocket_nozzle
    implicit none
    integer :: i, iostat_val
    double precision :: Pc, Tc, gamma, R_gas, M_mol, Pe, Pa
    double precision :: At, Ae, Dt, De, Me, Te, Ve, rho_e, mdot
    double precision :: CF, Isp, F_thrust, c_star, expansion_ratio
    double precision :: Pt_ratio, Tt_ratio, rho_ratio, a_ratio
    double precision :: M_i, AR_i, P_i, T_i, V_i, F_i, CF_i
    double precision :: R_specific, cp, a_throat, rho_c, V_throat
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: g0 = 9.80665d0

    read(*,*,iostat=iostat_val) Pc
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid chamber pressure input.'
        stop
    end if
    read(*,*,iostat=iostat_val) Tc
    read(*,*,iostat=iostat_val) gamma
    read(*,*,iostat=iostat_val) M_mol
    read(*,*,iostat=iostat_val) Pe
    read(*,*,iostat=iostat_val) Pa
    read(*,*,iostat=iostat_val) Dt
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all nozzle inputs.'
        stop
    end if
    if (Pc <= 0.0d0 .or. Tc <= 0.0d0 .or. gamma <= 1.0d0) then
        write(*,*) 'ERROR: Pc, Tc must be positive and gamma > 1.'
        stop
    end if
    if (M_mol <= 0.0d0) then
        write(*,*) 'ERROR: Molar mass must be positive.'
        stop
    end if
    if (Dt <= 0.0d0) then
        write(*,*) 'ERROR: Throat diameter must be positive.'
        stop
    end if
    if (Pe <= 0.0d0) Pe = 101325.0d0
    if (Pa < 0.0d0) Pa = 0.0d0

    R_specific = 8314.46d0 / M_mol   ! J/(kg K)
    cp = gamma * R_specific / (gamma - 1.0d0)
    At = PI / 4.0d0 * Dt**2

    ! Characteristic velocity c*
    c_star = sqrt(gamma * R_specific * Tc) / gamma / &
             sqrt((2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(gamma-1.0d0)))

    ! Mass flow rate (choked throat)
    rho_c = Pc / (R_specific * Tc)
    a_throat = sqrt(gamma * R_specific * Tc * &
               (2.0d0/(gamma+1.0d0)))
    mdot = Pc * At / c_star

    ! Exit Mach number from pressure ratio Pe/Pc
    Me = exit_mach(Pc, Pe, gamma)

    ! Expansion ratio Ae/At from exit Mach
    expansion_ratio = area_ratio(Me, gamma)
    Ae = At * expansion_ratio
    De = sqrt(4.0d0 * Ae / PI)

    ! Exit conditions
    Te = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * Me**2)
    Ve = Me * sqrt(gamma * R_specific * Te)
    rho_e = Pe / (R_specific * Te)

    ! Thrust
    F_thrust = mdot * Ve + (Pe - Pa) * Ae

    ! Thrust coefficient CF = F / (Pc At)
    CF = F_thrust / (Pc * At)

    ! Specific impulse
    Isp = F_thrust / (mdot * g0)

    ! Isentropic ratios at throat (M=1)
    Pt_ratio = (1.0d0 + (gamma-1.0d0)/2.0d0)**(-gamma/(gamma-1.0d0))
    Tt_ratio = 2.0d0 / (gamma + 1.0d0)
    V_throat = sqrt(gamma * R_specific * Tc * Tt_ratio)

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ROCKET NOZZLE DESIGN (DE LAVAL) ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Chamber Pressure Pc       = ', Pc, ' Pa'
    write(*,'(A,ES12.4,A)') '  Chamber Temperature Tc    = ', Tc, ' K'
    write(*,'(A,ES12.4)')   '  Specific Heat Ratio gamma = ', gamma
    write(*,'(A,ES12.4,A)') '  Molar Mass                = ', M_mol, ' kg/kmol'
    write(*,'(A,ES12.4,A)') '  R_specific                = ', R_specific, ' J/(kg.K)'
    write(*,'(A,ES12.4,A)') '  Exit Pressure Pe          = ', Pe, ' Pa'
    write(*,'(A,ES12.4,A)') '  Ambient Pressure Pa       = ', Pa, ' Pa'
    write(*,'(A,ES12.4,A)') '  Throat Diameter Dt        = ', Dt, ' m'
    write(*,*)
    write(*,'(A)') '--- NOZZLE GEOMETRY -----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Throat Area At            = ', At, ' m2'
    write(*,'(A,ES12.4,A)') '  Exit Area Ae              = ', Ae, ' m2'
    write(*,'(A,ES12.4,A)') '  Exit Diameter De          = ', De, ' m'
    write(*,'(A,ES12.4)')   '  Expansion Ratio Ae/At     = ', expansion_ratio
    write(*,*)
    write(*,'(A)') '--- EXIT CONDITIONS -----------------------------------------'
    write(*,'(A,ES12.4)')   '  Exit Mach Number Me       = ', Me
    write(*,'(A,ES12.4,A)') '  Exit Temperature Te       = ', Te, ' K'
    write(*,'(A,ES12.4,A)') '  Exit Velocity Ve          = ', Ve, ' m/s'
    write(*,'(A,ES12.4,A)') '  Exit Density              = ', rho_e, ' kg/m3'
    write(*,*)
    write(*,'(A)') '--- PERFORMANCE ---------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Mass Flow Rate            = ', mdot, ' kg/s'
    write(*,'(A,ES12.4,A)') '  Characteristic Velocity   = ', c_star, ' m/s'
    write(*,'(A,ES12.4,A)') '  Thrust                    = ', F_thrust, ' N'
    write(*,'(A,ES12.4,A)') '  Thrust                    = ', F_thrust/1000.0d0, ' kN'
    write(*,'(A,ES12.4)')   '  Thrust Coefficient CF     = ', CF
    write(*,'(A,ES12.4,A)') '  Specific Impulse Isp      = ', Isp, ' s'
    write(*,'(A,ES12.4,A)') '  Effective Exhaust Vel     = ', Isp*g0, ' m/s'
    write(*,*)
    write(*,'(A)') '--- THROAT CONDITIONS (M=1) ---------------------------------'
    write(*,'(A,ES12.4)')   '  P_throat / Pc             = ', Pt_ratio
    write(*,'(A,ES12.4)')   '  T_throat / Tc             = ', Tt_ratio
    write(*,'(A,ES12.4,A)') '  Throat Velocity           = ', V_throat, ' m/s'
    write(*,*)

    ! Nozzle profile: Mach, P, T, V vs A/At
    write(*,'(A)') '--- NOZZLE PROFILE ------------------------------------------'
    write(*,'(A)') '  A/At          Mach          P[Pa]         T[K]          V[m/s]'
    write(*,'(A)') '  ----------------------------------------------------------------'
    do i = 1, 80
        AR_i = 1.0d0 + (expansion_ratio - 1.0d0) * (dble(i-1)/79.0d0)**1.5
        M_i = mach_from_area_ratio(AR_i, gamma)
        T_i = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2)
        P_i = Pc * (T_i/Tc)**(gamma/(gamma-1.0d0))
        V_i = M_i * sqrt(gamma * R_specific * T_i)
        write(*,'(F10.4,2X,F10.5,2X,ES12.4,2X,F10.2,2X,F12.2)') &
            AR_i, M_i, P_i, T_i, V_i
    end do
    write(*,*)

    ! Thrust vs expansion ratio sweep
    write(*,'(A)') '--- THRUST VS EXPANSION RATIO SWEEP -------------------------'
    write(*,'(A)') '  Ae/At         Me            CF            F[kN]         Isp[s]'
    write(*,'(A)') '  ----------------------------------------------------------------'
    do i = 1, 50
        AR_i = 1.0d0 + 99.0d0 * dble(i-1) / 49.0d0
        M_i = mach_from_area_ratio(AR_i, gamma)
        T_i = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2)
        P_i = Pc * (T_i/Tc)**(gamma/(gamma-1.0d0))
        V_i = M_i * sqrt(gamma * R_specific * T_i)
        F_i = mdot * V_i + (P_i - Pa) * At * AR_i
        CF_i = F_i / (Pc * At)
        write(*,'(F10.2,2X,F10.5,2X,F10.5,2X,F12.4,2X,F10.2)') &
            AR_i, M_i, CF_i, F_i/1000.0d0, F_i/(mdot*g0)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Isentropic flow: T/Tc = 1/(1+(gamma-1)/2 M^2).'
    write(*,'(A)') '  Area-Mach: A/At = (1/M)((2+(gamma-1)M^2)/(gamma+1))^((gamma+1)/(2(gamma-1))).'
    write(*,'(A)') '  Thrust: F = mdot Ve + (Pe-Pa) Ae.'
    write(*,'(A)') '  c* = sqrt(gamma R Tc)/gamma / sqrt((2/(gamma+1))^((gamma+1)/(gamma-1))).'
    write(*,'(A)') '  CF = F/(Pc At); Isp = F/(mdot g0).'

contains

    double precision function exit_mach(Pcin, Pein, gam)
        implicit none
        double precision, intent(in) :: Pcin, Pein, gam
        double precision :: ratio, lo, hi, mid, P_mid
        integer :: it
        ratio = Pein / Pcin
        if (ratio >= 1.0d0) then
            exit_mach = 0.0d0
            return
        end if
        lo = 1.0d0
        hi = 50.0d0
        do it = 1, 200
            mid = 0.5d0*(lo+hi)
            P_mid = (1.0d0 + (gam-1.0d0)/2.0d0*mid**2)**(- gam/(gam-1.0d0))
            if (P_mid > ratio) then
                lo = mid
            else
                hi = mid
            end if
        end do
        exit_mach = 0.5d0*(lo+hi)
    end function exit_mach

    double precision function area_ratio(Min, gam)
        implicit none
        double precision, intent(in) :: Min, gam
        double precision :: term
        if (Min < 1.0d-10) then
            area_ratio = 1.0d30
            return
        end if
        term = (2.0d0 + (gam-1.0d0)*Min**2) / (gam+1.0d0)
        area_ratio = (1.0d0/Min) * term**((gam+1.0d0)/(2.0d0*(gam-1.0d0)))
    end function area_ratio

    double precision function mach_from_area_ratio(AR_in, gam)
        implicit none
        double precision, intent(in) :: AR_in, gam
        double precision :: lo, hi, mid, ar_mid
        integer :: it
        ! Supersonic branch
        lo = 1.0d0
        hi = 50.0d0
        do it = 1, 200
            mid = 0.5d0*(lo+hi)
            ar_mid = area_ratio(mid, gam)
            if (ar_mid < AR_in) then
                lo = mid
            else
                hi = mid
            end if
        end do
        mach_from_area_ratio = 0.5d0*(lo+hi)
    end function mach_from_area_ratio

end program rocket_nozzle