program restriction_orifice
    implicit none
    integer :: fluid_type, iostat_val, iter, n_holes
    double precision :: D_pipe_mm, W_kgh, P1_barA, dP_perm_bar, rho_kgm3, Pv_barA
    double precision :: D_pipe_m, Q_m3s, P2_barA, sigma_cav, do_mm, do_m
    double precision :: beta_ratio, Cd_val, dP_orifice_bar, t_plate_mm
    character(len=32) :: cav_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type      ! 1=Liquid, 2=Gas
    read(*,*,iostat=iostat_val) D_pipe_mm       ! Pipe Inside Diameter [mm]
    read(*,*,iostat=iostat_val) W_kgh           ! Mass Flow Rate [kg/h]
    read(*,*,iostat=iostat_val) P1_barA         ! Upstream Pressure P1 [bar a]
    read(*,*,iostat=iostat_val) dP_perm_bar     ! Required Permanent Delta P [bar]
    read(*,*,iostat=iostat_val) rho_kgm3        ! Fluid Density [kg/m3]
    read(*,*,iostat=iostat_val) Pv_barA         ! Vapor Pressure Pv [bar a] (for liquid)
    read(*,*,iostat=iostat_val) n_holes         ! Number of Orifice Holes (1=Single, >1=Multi-hole)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for restriction orifice sizing.'
        stop
    end if

    if (dP_perm_bar >= P1_barA .or. W_kgh <= 0.0d0 .or. D_pipe_mm <= 0.0d0) then
        write(*,*) 'ERROR: Delta P must be < P1 and flow rate must be positive.'
        stop
    end if

    D_pipe_m = D_pipe_mm * 1.0d-3
    Q_m3s = (W_kgh / 3600.0d0) / max(1.0d-1, rho_kgm3)
    P2_barA = P1_barA - dP_perm_bar

    ! Cavitation Index for liquids
    if (fluid_type == 1) then
        sigma_cav = (P2_barA - Pv_barA) / max(1.0d-2, dP_perm_bar)
        if (sigma_cav < 1.0d0) then
            cav_status = 'Severe Choked Cavitation / Flashing'
        else if (sigma_cav < 1.8d0) then
            cav_status = 'Incipient Cavitation Hazard'
        else
            cav_status = 'Safe / No Cavitation'
        end if
    else
        sigma_cav = 99.0d0
        cav_status = 'Compressible Gas Service'
    end if

    ! Iterative solution for beta ratio (Crane TP410)
    ! dP_perm = (1 - beta^1.9) * dP_orifice
    ! Q = Cd * A_o * sqrt(2 * dP_orifice / rho) / sqrt(1 - beta^4)
    beta_ratio = 0.45d0
    Cd_val = 0.62d0
    do iter = 1, 50
        dP_orifice_bar = dP_perm_bar / max(0.1d0, 1.0d0 - (beta_ratio**1.9d0))
        do_m = sqrt((4.0d0 * Q_m3s * sqrt(1.0d0 - beta_ratio**4)) / &
                    (PI * Cd_val * sqrt(2.0d0 * (dP_orifice_bar * 1.0d5) / rho_kgm3)))
        do_m = do_m / sqrt(dble(max(1, n_holes)))
        beta_ratio = (do_m * sqrt(dble(max(1, n_holes)))) / D_pipe_m
        if (beta_ratio > 0.85d0) beta_ratio = 0.85d0
        if (beta_ratio < 0.10d0) beta_ratio = 0.10d0
    end do
    do_mm = do_m * 1000.0d0

    ! Plate minimum thickness to resist bending stress (S_allow = 138 MPa)
    t_plate_mm = max(3.0d0, do_mm * sqrt((3.0d0 * (dP_perm_bar * 1.0d5)) / (4.0d0 * 138.0d6)) * 1000.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — THICK-PLATE RESTRICTION ORIFICE ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A)')  'Pipe Inside Diameter      = ', D_pipe_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Mass Flow Rate W          = ', W_kgh, ' kg/h'
    write(*,'(A,F10.2,A)')  'Upstream Pressure P1      = ', P1_barA, ' bar(a)'
    write(*,'(A,F10.2,A)')  'Permanent Pressure Loss   = ', dP_perm_bar, ' bar'
    write(*,'(A,I2)')       'Number of Holes           = ', n_holes
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Orifice Hole Diameter do  = ', do_mm, ' mm each'
    write(*,'(A,F10.4)')    'Diameter Ratio Beta (d/D) = ', beta_ratio
    write(*,'(A,F10.2)')    'Cavitation Index Sigma    = ', sigma_cav
    write(*,'(A,A)')        'Cavitation Risk Status    = ', trim(cav_status)
    write(*,'(A,F10.2,A)')  'Recommended Plate Thk     = ', t_plate_mm, ' mm'
    write(*,'(A)') '============================================================'

end program restriction_orifice
