program venturi_nozzle
    implicit none
    integer :: i, iter, iostat_val
    double precision :: Dpipe, Dthroat, rho, mu, deltaP, gamma_exp, Cd_user
    double precision :: beta, A2, E_approach, Cd, Q, mdot, V2, Re_D, Re_throat, nu
    double precision :: Cd_prev, Q_prev, dP_i, Q_i, Re_i, Cd_i
    double precision, parameter :: PI = 3.141592653589793d0

    read(*,*,iostat=iostat_val) Dpipe
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid pipe diameter input.'
        stop
    end if
    read(*,*,iostat=iostat_val) Dthroat
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) deltaP
    read(*,*,iostat=iostat_val) gamma_exp
    read(*,*,iostat=iostat_val) Cd_user
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all venturi inputs.'
        stop
    end if
    if (Dpipe <= 0.0d0 .or. Dthroat <= 0.0d0) then
        write(*,*) 'ERROR: Pipe and throat diameters must be positive.'
        stop
    end if
    if (Dthroat >= Dpipe) then
        write(*,*) 'ERROR: Throat diameter must be smaller than pipe diameter.'
        stop
    end if
    if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. deltaP < 0.0d0) then
        write(*,*) 'ERROR: Density, viscosity must be positive; deltaP non-negative.'
        stop
    end if
    if (gamma_exp <= 0.0d0 .or. gamma_exp > 1.0d0) gamma_exp = 1.0d0

    beta = Dthroat / Dpipe
    A2 = PI / 4.0d0 * Dthroat**2
    E_approach = 1.0d0 / sqrt(1.0d0 - beta**4)
    nu = mu / rho

    ! Iterate for Cd and Q
    if (Cd_user > 0.01d0) then
        Cd = Cd_user
    else
        Cd = 0.995d0
    end if

    do iter = 1, 100
        Cd_prev = Cd
        Q = Cd * gamma_exp * E_approach * A2 * sqrt(2.0d0 * deltaP / rho)
        V2 = Q / A2
        Re_D = rho * (Q / (PI/4.0d0*Dpipe**2)) * Dpipe / mu
        Re_throat = rho * V2 * Dthroat / mu
        if (Cd_user > 0.01d0) then
            Cd = Cd_user
        else
            Cd = venturi_cd(Re_D, beta)
        end if
        if (abs(Cd - Cd_prev) < 1.0d-10) exit
    end do

    Q = Cd * gamma_exp * E_approach * A2 * sqrt(2.0d0 * deltaP / rho)
    mdot = rho * Q
    V2 = Q / A2
    Re_D = rho * (Q / (PI/4.0d0*Dpipe**2)) * Dpipe / mu
    Re_throat = rho * V2 * Dthroat / mu

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   VENTURI & NOZZLE FLOW METER ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- GEOMETRY ------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Pipe Diameter            = ', Dpipe, ' m'
    write(*,'(A,ES12.4,A)') '  Throat Diameter          = ', Dthroat, ' m'
    write(*,'(A,ES12.4)')   '  Beta Ratio               = ', beta
    write(*,'(A,ES12.4,A)') '  Throat Area              = ', A2, ' m2'
    write(*,'(A,ES12.4)')   '  Approach Velocity Factor = ', E_approach
    write(*,*)
    write(*,'(A)') '--- FLOW RESULTS --------------------------------------------'
    write(*,'(A,ES12.4)')   '  Discharge Coefficient Cd = ', Cd
    write(*,'(A,ES12.4)')   '  Expansion Factor         = ', gamma_exp
    write(*,'(A,ES12.4,A)') '  Volume Flow Q            = ', Q, ' m3/s'
    write(*,'(A,ES12.4,A)') '  Mass Flow                = ', mdot, ' kg/s'
    write(*,'(A,ES12.4,A)') '  Throat Velocity          = ', V2, ' m/s'
    write(*,'(A,ES12.4)')   '  Re (pipe)                = ', Re_D
    write(*,'(A,ES12.4)')   '  Re (throat)              = ', Re_throat
    write(*,'(A,ES12.4,A)') '  Differential Pressure    = ', deltaP, ' Pa'
    write(*,*)

    write(*,'(A)') '--- CD VS RE SWEEP ------------------------------------------'
    write(*,'(A)') '  Re_D          Cd'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 50
        Re_i = 1.0d3 * (1.0d7/1.0d3)**(dble(i-1)/49.0d0)
        Cd_i = venturi_cd(Re_i, beta)
        write(*,'(ES12.4,2X,ES12.4)') Re_i, Cd_i
    end do
    write(*,*)

    write(*,'(A)') '--- Q VS DELTAP SWEEP ---------------------------------------'
    write(*,'(A)') '  deltaP[Pa]    Q[m3/s]'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 40
        dP_i = max(deltaP, 100.0d0) * 0.05d0 * dble(i)
        Q_i = Cd * gamma_exp * E_approach * A2 * sqrt(2.0d0*dP_i/rho)
        write(*,'(ES12.4,2X,ES12.4)') dP_i, Q_i
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Q = Cd * gamma * E * A2 * sqrt(2 deltaP / rho).'
    write(*,'(A)') '  E = 1/sqrt(1 - beta^4); beta = D_throat/D_pipe.'
    write(*,'(A)') '  Cd = simplified Reader-Harris/Gallagher for venturi.'

contains
    double precision function venturi_cd(Re_in, beta_in)
        implicit none
        double precision, intent(in) :: Re_in, beta_in
        ! Simplified ISO 5167 venturi Cd correlation
        double precision :: base
        base = 0.9858d0 - 0.196d0 * beta_in**4.5d0
        if (Re_in > 2.0d5) then
            venturi_cd = base
        else if (Re_in > 1.0d4) then
            venturi_cd = base - 0.005d0 * (2.0d5 - Re_in) / 1.9d5
        else
            venturi_cd = base - 0.01d0
        end if
        if (venturi_cd < 0.90d0) venturi_cd = 0.90d0
        if (venturi_cd > 1.00d0) venturi_cd = 1.00d0
    end function venturi_cd
end program venturi_nozzle