program droplet_evaporation
    implicit none
    integer :: i, iostat_val, n_points
    double precision :: d0_um, d0_m, T_inf_C, T_inf_K, Ts_C, Ts_K
    double precision :: U_rel, rho_L, Lv_kJ, Lv_J, kg, Cpg, rho_g, mug, Dab
    double precision :: Pr, Sc, BT, Re0, F_conv, K_evap, tau_life, mdot0
    double precision :: t_curr, d_curr, m_curr, X_stop
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) d0_um
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid droplet diameter.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_inf_C
    read(*,*,iostat=iostat_val) Ts_C
    read(*,*,iostat=iostat_val) U_rel
    read(*,*,iostat=iostat_val) rho_L
    read(*,*,iostat=iostat_val) Lv_kJ
    read(*,*,iostat=iostat_val) kg
    read(*,*,iostat=iostat_val) Cpg
    read(*,*,iostat=iostat_val) rho_g
    read(*,*,iostat=iostat_val) mug
    read(*,*,iostat=iostat_val) Dab

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

    ! Validation
    if (d0_um <= 0.0d0 .or. rho_L <= 0.0d0 .or. Lv_kJ <= 0.0d0) then
        write(*,*) 'ERROR: Diameter, density, and latent heat must be positive.'
        stop
    end if
    if (T_inf_C <= Ts_C) then
        write(*,*) 'ERROR: Ambient gas temperature must exceed droplet surface temperature.'
        stop
    end if
    if (kg <= 0.0d0 .or. Cpg <= 0.0d0 .or. rho_g <= 0.0d0 .or. mug <= 0.0d0) then
        write(*,*) 'ERROR: Gas transport properties must be positive.'
        stop
    end if

    ! Conversions
    d0_m = d0_um * 1.0d-6
    T_inf_K = T_inf_C + 273.15d0
    Ts_K = Ts_C + 273.15d0
    Lv_J = Lv_kJ * 1.0d3

    ! Dimensionless Numbers
    Pr = (Cpg * mug) / kg
    if (Dab > 0.0d0) then
        Sc = mug / (rho_g * Dab)
    else
        Sc = Pr
    end if
    Re0 = (rho_g * U_rel * d0_m) / mug

    ! Spalding Transfer Number
    BT = (Cpg * (T_inf_K - Ts_K)) / Lv_J

    ! Ranz-Marshall Convective Factor
    F_conv = 1.0d0 + 0.3d0 * sqrt(max(0.0d0, Re0)) * (Pr**(1.0d0 / 3.0d0))

    ! Evaporation Constant K (d2-law) in m2/s
    K_evap = (8.0d0 * kg / (rho_L * Cpg)) * log(1.0d0 + BT) * F_conv

    ! Lifetime & Initial Evaporation Rate
    if (K_evap > 1.0d-30) then
        tau_life = (d0_m**2) / K_evap ! seconds
    else
        tau_life = 1.0d30
    end if
    mdot0 = (PI / 2.0d0) * rho_L * d0_m * K_evap ! kg/s

    ! Stopping Distance
    X_stop = (rho_L * (d0_m**2) * U_rel) / (18.0d0 * mug * (1.0d0 + 0.15d0 * (max(0.01d0, Re0)**0.687d0)))

    ! Output Formatted Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — DROPLET EVAPORATION & D2-LAW ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,ES12.4,A)') 'Initial Diameter d0       = ', d0_um, ' um'
    write(*,'(A,ES12.4,A)') 'Ambient Gas Temp T_inf    = ', T_inf_C, ' C'
    write(*,'(A,ES12.4,A)') 'Droplet Surface Temp Ts   = ', Ts_C, ' C'
    write(*,'(A,ES12.4,A)') 'Relative Velocity U_rel   = ', U_rel, ' m/s'
    write(*,'(A,ES12.4,A)') 'Liquid Density            = ', rho_L, ' kg/m3'
    write(*,'(A,ES12.4,A)') 'Latent Heat of Vapor. Lv  = ', Lv_kJ, ' kJ/kg'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4)')   'Prandtl Number Pr         = ', Pr
    write(*,'(A,ES12.4)')   'Schmidt Number Sc         = ', Sc
    write(*,'(A,ES12.4)')   'Initial Reynolds Re0      = ', Re0
    write(*,'(A,ES12.4)')   'Spalding Transfer Num BT  = ', BT
    write(*,'(A,ES12.4)')   'Ranz-Marshall Factor      = ', F_conv
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,ES12.4,A)') 'Evaporation Constant K    = ', K_evap * 1.0d6, ' mm2/s'
    write(*,'(A,ES12.4,A)') 'Total Droplet Lifetime    = ', tau_life * 1.0d3, ' ms'
    write(*,'(A,ES12.4,A)') 'Initial Mass Evap. Rate   = ', mdot0, ' kg/s'
    write(*,'(A,ES12.4,A)') 'Aerodynamic Stop Distance = ', X_stop * 1.0d3, ' mm'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- DROPLET DIAMETER AND MASS VS TIME ---'
    write(*,'(A)') '  Time(ms)     d(um)        d2/d02       Mass(ug)'
    n_points = 20
    do i = 0, n_points
        t_curr = tau_life * (dble(i) / dble(n_points))
        if (t_curr >= tau_life) then
            d_curr = 0.0d0
        else
            d_curr = sqrt(max(0.0d0, (d0_m**2) - K_evap * t_curr))
        end if
        m_curr = (PI / 6.0d0) * rho_L * (d_curr**3)
        write(*,'(2X,F10.2,2X,F10.2,2X,F10.4,2X,ES12.4)') t_curr * 1.0d3, d_curr * 1.0d6, &
            (d_curr / d0_m)**2, m_curr * 1.0d9
    end do

end program droplet_evaporation