program spray_atomization
    implicit none
    integer :: i, iostat_val, n_points
    double precision :: d_nozzle_mm, d_m, deltaP_bar, deltaP_Pa, Cd
    double precision :: rho_L, mu_L, sigma, rho_g
    double precision :: A_orifice, V_ideal, V_inj, mdot_L, Q_L
    double precision :: Re_L, We_L, We_g, Oh
    double precision :: term1, term2, smd_m, D32, D10, Dv50
    double precision :: rho_ratio, theta_rad, spray_angle, Lb_m, breakup_length
    double precision :: q_rr, X_rr, d_val, cum_vol, pdf_val
    double precision :: p_var, p_var_Pa, m_var, smd_var
    double precision, parameter :: PI = 3.141592653589793d0
    character(len=64) :: regime_name

    ! Read inputs from standard input
    read(*,*,iostat=iostat_val) d_nozzle_mm
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid nozzle diameter input.'
        stop
    end if
    read(*,*,iostat=iostat_val) deltaP_bar
    read(*,*,iostat=iostat_val) Cd
    read(*,*,iostat=iostat_val) rho_L
    read(*,*,iostat=iostat_val) mu_L
    read(*,*,iostat=iostat_val) sigma
    read(*,*,iostat=iostat_val) rho_g

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

    ! Validation
    if (d_nozzle_mm <= 0.0d0 .or. deltaP_bar <= 0.0d0 .or. Cd <= 0.0d0) then
        write(*,*) 'ERROR: Geometric and pressure inputs must be positive.'
        stop
    end if
    if (rho_L <= 0.0d0 .or. mu_L <= 0.0d0 .or. sigma <= 0.0d0 .or. rho_g <= 0.0d0) then
        write(*,*) 'ERROR: Physical fluid properties must be positive.'
        stop
    end if

    ! Conversions
    d_m = d_nozzle_mm * 1.0d-3
    deltaP_Pa = deltaP_bar * 1.0d5
    A_orifice = (PI / 4.0d0) * (d_m**2)

    ! Hydrodynamics & Flow Rates
    V_ideal = sqrt(2.0d0 * deltaP_Pa / rho_L)
    V_inj = Cd * V_ideal
    mdot_L = Cd * A_orifice * sqrt(2.0d0 * rho_L * deltaP_Pa)
    Q_L = mdot_L / rho_L

    ! Non-dimensional numbers
    Re_L = (rho_L * V_inj * d_m) / mu_L
    We_L = (rho_L * (V_inj**2) * d_m) / sigma
    We_g = (rho_g * (V_inj**2) * d_m) / sigma
    Oh = mu_L / sqrt(rho_L * sigma * d_m)

    ! Breakup Regime Identification
    if (Oh < 0.1d0 .and. We_L < 10.0d0) then
        regime_name = 'Rayleigh Capillary Breakup'
    else if (We_g < 13.0d0) then
        regime_name = 'First Wind-Induced Regime'
    else if (We_g < 40.3d0) then
        regime_name = 'Second Wind-Induced Regime'
    else
        regime_name = 'Catastrophic Shear Atomization'
    end if

    ! Sauter Mean Diameter (SMD D32) — Lefebvre-Hiroyasu Hybrid Model
    term1 = 2.25d0 * (sigma**0.25d0) * (mu_L**0.16d0) * (rho_L**0.20d0)
    term2 = (mdot_L**0.22d0) * (deltaP_Pa**(-0.43d0))
    smd_m = term1 * term2

    if (We_g > 40.0d0) then
        smd_m = smd_m * ((40.0d0 / We_g)**0.12d0)
    end if

    D32 = max(1.0d0, smd_m * 1.0d6) ! in micrometers
    D10 = D32 * 0.65d0
    Dv50 = D32 * 1.18d0 ! Mass Median Diameter

    ! Spray Cone Angle (Arai/Hiroyasu)
    rho_ratio = rho_g / rho_L
    theta_rad = 0.057d0 * (rho_ratio**0.10d0) * (max(100.0d0, Re_L)**0.25d0)
    spray_angle = min(120.0d0, max(5.0d0, theta_rad * (180.0d0 / PI) * 1.8d0))

    ! Liquid Core Length
    Lb_m = 7.0d0 * d_m * sqrt(rho_L / rho_g) * (1.0d0 + 0.4d0 * Oh)
    breakup_length = Lb_m * 1.0d3 ! in mm

    ! Rosin-Rammler parameters
    q_rr = 2.85d0
    X_rr = Dv50 / (log(2.0d0)**(1.0d0 / q_rr))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SPRAY & DROPLET ATOMIZATION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,ES12.4,A)') 'Nozzle Diameter d0        = ', d_nozzle_mm, ' mm'
    write(*,'(A,ES12.4,A)') 'Injection Delta-P         = ', deltaP_bar, ' bar'
    write(*,'(A,ES12.4)')   'Discharge Coefficient Cd  = ', Cd
    write(*,'(A,ES12.4,A)') 'Liquid Density            = ', rho_L, ' kg/m3'
    write(*,'(A,ES12.4,A)') 'Liquid Viscosity          = ', mu_L, ' Pa.s'
    write(*,'(A,ES12.4,A)') 'Surface Tension           = ', sigma, ' N/m'
    write(*,'(A,ES12.4,A)') 'Gas Density               = ', rho_g, ' kg/m3'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Injection Velocity        = ', V_inj, ' m/s'
    write(*,'(A,ES12.4,A)') 'Mass Flow Rate            = ', mdot_L, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Volumetric Flow Rate      = ', Q_L, ' m3/s'
    write(*,'(A,ES12.4)')   'Liquid Reynolds Number    = ', Re_L
    write(*,'(A,ES12.4)')   'Gas Weber Number Weg      = ', We_g
    write(*,'(A,ES12.4)')   'Ohnesorge Number Oh       = ', Oh
    write(*,'(A,A)')        'Breakup Regime            = ', trim(regime_name)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Sauter Mean Diameter D32  = ', D32, ' um'
    write(*,'(A,ES12.4,A)') 'Arithmetic Mean D10       = ', D10, ' um'
    write(*,'(A,ES12.4,A)') 'Mass Median Diameter Dv50 = ', Dv50, ' um'
    write(*,'(A,ES12.4,A)') 'Spray Cone Angle          = ', spray_angle, ' deg'
    write(*,'(A,ES12.4,A)') 'Liquid Core Breakup Length= ', breakup_length, ' mm'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- ROSIN-RAMMLER DROPLET SIZE DISTRIBUTION ---'
    write(*,'(A)') '  d_um         Q_cum(%)     PDF'
    n_points = 25
    do i = 0, n_points
        d_val = (Dv50 * 2.8d0) * (dble(i) / dble(n_points))
        if (d_val <= 1.0d-6) then
            cum_vol = 0.0d0
            pdf_val = 0.0d0
        else
            cum_vol = 1.0d0 - exp(-((d_val / X_rr)**q_rr))
            pdf_val = (q_rr / X_rr) * ((d_val / X_rr)**(q_rr - 1.0d0)) * exp(-((d_val / X_rr)**q_rr))
        end if
        write(*,'(2X,F10.2,2X,F10.2,2X,ES12.4)') d_val, cum_vol * 100.0d0, pdf_val * 100.0d0
    end do
    write(*,*)
    write(*,'(A)') '--- SMD D32 VS INJECTION PRESSURE ---'
    write(*,'(A)') '  DeltaP_bar   D32_um'
    do i = 1, 15
        p_var = (deltaP_bar * 0.2d0) + (deltaP_bar * 2.8d0) * (dble(i - 1) / 14.0d0)
        p_var_Pa = p_var * 1.0d5
        m_var = Cd * A_orifice * sqrt(2.0d0 * rho_L * p_var_Pa)
        smd_var = (term1 * (m_var**0.22d0) * (p_var_Pa**(-0.43d0))) * 1.0d6
        write(*,'(2X,F10.2,2X,F10.2)') p_var, max(1.0d0, smd_var)
    end do

end program spray_atomization
