program steam_trap_sizing
    implicit none
    integer :: iostat_val, trap_type
    double precision :: P1_barG, P2_barG, mdot_c_kgh, safety_factor
    double precision :: P1_barA, P2_barA, hf1, hf2, hfg2, vg2
    double precision :: x_flash_frac, mdot_flash_kgh, Q_flash_m3h
    double precision :: D_pipe_mm, D_pipe_m, V_mix_ms, design_cap_kgh
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) P1_barG        ! Steam Supply Pressure [bar g]
    read(*,*,iostat=iostat_val) P2_barG        ! Condensate Return Pressure [bar g]
    read(*,*,iostat=iostat_val) mdot_c_kgh     ! Running Condensate Load [kg/h]
    read(*,*,iostat=iostat_val) safety_factor  ! Sizing Safety Factor (1.5 - 3.0)
    read(*,*,iostat=iostat_val) trap_type      ! 1=Float-Thermostatic, 2=Thermodynamic, 3=Inverted Bucket

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs for steam trap sizing.'
        stop
    end if

    if (P1_barG <= P2_barG) then
        write(*,*) 'ERROR: Steam inlet pressure must be greater than return pressure.'
        stop
    end if

    P1_barA = P1_barG + 1.01325d0
    P2_barA = P2_barG + 1.01325d0

    ! Empirical Steam Enthalpies [kJ/kg] & Specific Volume
    hf1 = 419.0d0 * (P1_barA**0.235d0)
    hf2 = 419.0d0 * (P2_barA**0.235d0)
    hfg2 = 2257.0d0 / (P2_barA**0.075d0)
    vg2 = 1.673d0 / (P2_barA**0.95d0) ! m3/kg

    ! Flash Steam Fraction
    x_flash_frac = max(0.0d0, (hf1 - hf2) / hfg2)
    mdot_flash_kgh = mdot_c_kgh * x_flash_frac
    Q_flash_m3h = mdot_flash_kgh * vg2

    ! Rated Trap Sizing Capacity
    design_cap_kgh = mdot_c_kgh * safety_factor

    ! Condensate Return Line Sizing (Target mixture velocity <= 15 m/s)
    ! Flow area A = (Q_flash_m3s) / V_target
    D_pipe_m = sqrt((4.0d0 * (Q_flash_m3h / 3600.0d0)) / (PI * 12.0d0))
    D_pipe_mm = max(20.0d0, D_pipe_m * 1000.0d0)
    V_mix_ms = (Q_flash_m3h / 3600.0d0) / (PI * (D_pipe_mm * 1.0d-3)**2 / 4.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — STEAM TRAP & FLASH STEAM RECOVERY ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Steam Supply Pressure P1  = ', P1_barG, ' bar(g)'
    write(*,'(A,F10.2,A)')  'Condensate Return P2      = ', P2_barG, ' bar(g)'
    write(*,'(A,F10.2,A)')  'Running Condensate Load   = ', mdot_c_kgh, ' kg/h'
    write(*,'(A,F10.2,A)')  'Trap Sizing Capacity (SF) = ', design_cap_kgh, ' kg/h'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Flash Steam Percentage    = ', x_flash_frac * 100.0d0, ' %'
    write(*,'(A,F10.2,A)')  'Flash Steam Generated     = ', mdot_flash_kgh, ' kg/h'
    write(*,'(A,F10.2,A)')  'Flash Vapor Volume Flow   = ', Q_flash_m3h, ' m3/h'
    write(*,'(A,F10.1,A)')  'Recommended Return Pipe ID= ', D_pipe_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Condensate Line Velocity  = ', V_mix_ms, ' m/s'
    write(*,'(A)') '============================================================'

end program steam_trap_sizing