program ejector_eductor
    implicit none
    integer :: i, iostat_val, n_points
    double precision :: Pm_bar, Ps_bar, Pout_bar, Pm_Pa, Ps_Pa, Pout_Pa
    double precision :: dm_mm, dmix_mm, dm_m, dmix_m, Am, Amix, As, area_ratio
    double precision :: rho_m, rho_s, eta_nozzle, eta_diff
    double precision :: deltaP_motive, deltaP_lift, V_motive, mdot_m, Q_m
    double precision :: N_lift, Rm_ideal_max, Rm_calc, Rm, mdot_s, Q_s, Rv
    double precision :: mdot_out, rho_mix, V_mix, power_fluid_in, power_fluid_lift, efficiency
    double precision :: pm_var, dp_var, n_var, rm_var, rm_curr, n_curr
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Pm_bar
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid motive pressure.'
        stop
    end if
    read(*,*,iostat=iostat_val) Ps_bar
    read(*,*,iostat=iostat_val) Pout_bar
    read(*,*,iostat=iostat_val) dm_mm
    read(*,*,iostat=iostat_val) dmix_mm
    read(*,*,iostat=iostat_val) rho_m
    read(*,*,iostat=iostat_val) rho_s
    read(*,*,iostat=iostat_val) eta_nozzle
    read(*,*,iostat=iostat_val) eta_diff

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all ejector inputs.'
        stop
    end if

    ! Validation
    if (Pm_bar <= Ps_bar) then
        write(*,*) 'ERROR: Motive pressure must exceed suction pressure.'
        stop
    end if
    if (Pout_bar < Ps_bar .or. Pout_bar >= Pm_bar) then
        write(*,*) 'ERROR: Discharge pressure must be between suction and motive.'
        stop
    end if
    if (dm_mm <= 0.0d0 .or. dmix_mm <= dm_mm) then
        write(*,*) 'ERROR: Nozzle diameter must be positive and < mixing diameter.'
        stop
    end if

    ! Conversions
    Pm_Pa = Pm_bar * 1.0d5
    Ps_Pa = Ps_bar * 1.0d5
    Pout_Pa = Pout_bar * 1.0d5
    dm_m = dm_mm * 1.0d-3
    dmix_m = dmix_mm * 1.0d-3

    Am = (PI / 4.0d0) * (dm_m**2)
    Amix = (PI / 4.0d0) * (dmix_m**2)
    area_ratio = Am / Amix

    deltaP_motive = Pm_Pa - Ps_Pa
    deltaP_lift = Pout_Pa - Ps_Pa

    ! Primary Motive Jet
    V_motive = eta_nozzle * sqrt(2.0d0 * deltaP_motive / rho_m)
    mdot_m = rho_m * Am * V_motive
    Q_m = mdot_m / rho_m

    ! 1D Momentum Balance (Cunningham / ESDU 85032)
    N_lift = deltaP_lift / deltaP_motive
    Rm_ideal_max = ((1.0d0 / area_ratio) - 1.0d0) * sqrt(rho_s / rho_m)
    Rm_calc = Rm_ideal_max * max(0.0d0, (1.0d0 - (N_lift / max(0.05d0, 2.0d0 * area_ratio * (1.0d0 - area_ratio)))))
    Rm = max(0.0d0, min(Rm_ideal_max, Rm_calc * 0.90d0))

    mdot_s = Rm * mdot_m
    Q_s = mdot_s / rho_s
    if (Q_m > 0.0d0) then
        Rv = Q_s / Q_m
    else
        Rv = 0.0d0
    end if

    mdot_out = mdot_m + mdot_s
    if (Q_m + Q_s > 0.0d0) then
        rho_mix = mdot_out / (Q_m + Q_s)
    else
        rho_mix = rho_m
    end if
    V_mix = mdot_out / (rho_mix * Amix)

    power_fluid_in = Q_m * (Pm_Pa - Pout_Pa)
    power_fluid_lift = Q_s * (Pout_Pa - Ps_Pa)
    if (power_fluid_in > 0.0d0) then
        efficiency = min(65.0d0, max(0.0d0, (power_fluid_lift / power_fluid_in) * 100.0d0))
    else
        efficiency = 0.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — EJECTOR & EDUCTOR JET PUMP ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,ES12.4,A)') 'Motive Pressure Pm        = ', Pm_bar, ' bar'
    write(*,'(A,ES12.4,A)') 'Suction Pressure Ps       = ', Ps_bar, ' bar'
    write(*,'(A,ES12.4,A)') 'Discharge Pressure Pout   = ', Pout_bar, ' bar'
    write(*,'(A,ES12.4,A)') 'Motive Nozzle Diameter    = ', dm_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Mixing Chamber Diameter   = ', dmix_mm, ' mm'
    write(*,'(A,ES12.4)')   'Area Ratio Am/Amix        = ', area_ratio
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Motive Jet Velocity       = ', V_motive, ' m/s'
    write(*,'(A,ES12.4,A)') 'Motive Mass Flow (mdot_m) = ', mdot_m, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Suction Mass Flow (mdot_s)= ', mdot_s, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Total Discharge Mass Flow = ', mdot_out, ' kg/s'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4)')   'Entrainment Mass Ratio Rm = ', Rm
    write(*,'(A,ES12.4)')   'Volumetric Ratio Rv       = ', Rv
    write(*,'(A,ES12.4)')   'Pressure Lift Ratio N     = ', N_lift
    write(*,'(A,ES12.4,A)') 'Hydraulic Efficiency      = ', efficiency, ' %'
    write(*,'(A,ES12.4,A)') 'Throat Mixing Velocity    = ', V_mix, ' m/s'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- ENTRAINMENT RATIO RM VS MOTIVE PRESSURE PM ---'
    write(*,'(A)') '  Pm_bar       Rm(kg/kg)'
    n_points = 15
    do i = 1, n_points
        pm_var = (Ps_bar * 1.2d0) + (Pm_bar * 2.5d0 - Ps_bar * 1.2d0) * (dble(i - 1) / dble(n_points - 1))
        dp_var = (pm_var - Ps_bar) * 1.0d5
        n_var = (Pout_Pa - Ps_Pa) / max(1.0d0, dp_var)
        rm_var = Rm_ideal_max * max(0.0d0, (1.0d0 - (n_var / max(0.05d0, 2.0d0 * area_ratio * (1.0d0 - area_ratio))))) * 0.90d0
        write(*,'(2X,F10.2,2X,F10.4)') pm_var, max(0.0d0, rm_var)
    end do
    write(*,*)
    write(*,'(A)') '--- CHARACTERISTIC PUMP CURVE: PRESSURE LIFT N VS RM ---'
    write(*,'(A)') '  Rm(kg/kg)    Pressure_Lift_N'
    do i = 0, n_points
        rm_curr = (Rm_ideal_max * 1.1d0) * (dble(i) / dble(n_points))
        n_curr = max(0.0d0, 2.0d0 * area_ratio * (1.0d0 - area_ratio) * (1.0d0 - (rm_curr / max(1.0d-4, Rm_ideal_max))))
        write(*,'(2X,F10.4,2X,F10.5)') rm_curr, n_curr
    end do

end program ejector_eductor
