program duct_sizing_regain
    implicit none
    integer :: method_mode, shape_mode, iostat_val, iter
    double precision :: Q_m3h, dP_L_target_Pam, v_target_ms, rect_width_mm
    double precision :: Q_m3s, rho_air, mu_air, eps_rough_m
    double precision :: D_circ_m, D_circ_mm, v_actual_ms, Re_air, f_darcy
    double precision :: dP_L_actual, Pv_Pa, b_height_mm, b_height_m, De_calc_m, aspect_ratio
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) method_mode      ! 1=Equal Friction, 2=Velocity Method
    read(*,*,iostat=iostat_val) shape_mode       ! 1=Circular, 2=Rectangular
    read(*,*,iostat=iostat_val) Q_m3h            ! Air Volume Flow Rate [m3/h]
    read(*,*,iostat=iostat_val) dP_L_target_Pam  ! Target Friction Rate [Pa/m] (e.g. 1.0)
    read(*,*,iostat=iostat_val) v_target_ms      ! Target Air Velocity [m/s] (e.g. 6.0)
    read(*,*,iostat=iostat_val) rect_width_mm    ! Constrained Rectangular Width a [mm]

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

    Q_m3s = Q_m3h / 3600.0d0
    rho_air = 1.204d0     ! Standard air at 20 C [kg/m3]
    mu_air = 1.81d-5      ! Air viscosity [Pa.s]
    eps_rough_m = 0.09d-3 ! Galvanized steel duct roughness [m]

    ! 1. Sizing Circular Duct
    if (method_mode == 1) then
        ! Equal Friction sizing (iterate on D to match dP/L)
        D_circ_m = 0.30d0
        do iter = 1, 40
            v_actual_ms = Q_m3s / ((PI / 4.0d0) * (D_circ_m**2))
            Re_air = (rho_air * v_actual_ms * D_circ_m) / mu_air
            f_darcy = 0.25d0 / (log10((eps_rough_m / (3.7d0 * D_circ_m)) + (5.74d0 / (max(1000.0d0, Re_air)**0.9d0)))**2)
            dP_L_actual = f_darcy * (rho_air * (v_actual_ms**2)) / (2.0d0 * D_circ_m)
            D_circ_m = D_circ_m * ((dP_L_actual / max(0.01d0, dP_L_target_Pam))**0.20d0)
        end do
    else
        ! Velocity Sizing
        D_circ_m = sqrt((4.0d0 * Q_m3s) / (PI * max(0.5d0, v_target_ms)))
    end if

    D_circ_mm = D_circ_m * 1000.0d0
    v_actual_ms = Q_m3s / ((PI / 4.0d0) * (D_circ_m**2))
    Re_air = (rho_air * v_actual_ms * D_circ_m) / mu_air
    f_darcy = 0.25d0 / (log10((eps_rough_m / (3.7d0 * D_circ_m)) + (5.74d0 / (max(1000.0d0, Re_air)**0.9d0)))**2)
    dP_L_actual = f_darcy * (rho_air * (v_actual_ms**2)) / (2.0d0 * D_circ_m)
    Pv_Pa = 0.5d0 * rho_air * (v_actual_ms**2)

    ! 2. Rectangular Duct Equivalence (Huebscher formula)
    ! De = 1.30 * (a * b)^0.625 / (a + b)^0.25
    b_height_mm = 200.0d0
    if (shape_mode == 2 .and. rect_width_mm > 10.0d0) then
        do iter = 1, 50
            De_calc_m = 1.30d0 * (((rect_width_mm * 1.0d-3 * b_height_mm * 1.0d-3)**0.625d0) / &
                        ((rect_width_mm * 1.0d-3 + b_height_mm * 1.0d-3)**0.25d0))
            b_height_mm = b_height_mm * (D_circ_m / max(0.01d0, De_calc_m))
        end do
        aspect_ratio = max(rect_width_mm, b_height_mm) / min(rect_width_mm, b_height_mm)
    else
        b_height_mm = D_circ_mm
        aspect_ratio = 1.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ASHRAE / SMACNA DUCT SIZING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A)')  'Air Flow Rate Q           = ', Q_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  'Equivalent Circular Diam  = ', D_circ_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Air Velocity in Duct      = ', v_actual_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Velocity Pressure Pv      = ', Pv_Pa, ' Pa'
    write(*,'(A,F10.3,A)')  'Friction Loss Rate dP/L   = ', dP_L_actual, ' Pa/m'
    write(*,'(A,ES12.4)')   'Reynolds Number Re        = ', Re_air
    if (shape_mode == 2) then
        write(*,'(A)') '------------------------------------------------------------'
        write(*,'(A,F10.1,A,F10.1,A)') 'Rectangular Duct Size     = ', rect_width_mm, ' x ', b_height_mm, ' mm'
        write(*,'(A,F10.2,A)')  'Duct Aspect Ratio (a/b)   = ', aspect_ratio, ':1'
    end if
    write(*,'(A)') '============================================================'

end program duct_sizing_regain
