program pipe_networks
    implicit none
    
    integer, parameter :: max_pipes = 8
    double precision :: rho, mu, nu, Q_total, dZ
    double precision :: L(max_pipes), D(max_pipes), eps(max_pipes)
    double precision :: Q(max_pipes), V(max_pipes), Re(max_pipes), f(max_pipes)
    double precision :: dP(max_pipes), hf(max_pipes)
    double precision :: rel_rough(max_pipes), Ac(max_pipes)
    integer :: n_pipes, net_type  ! 1=series, 2=parallel
    double precision :: pi, g
    
    ! Hardy-Cross variables
    double precision :: dQ, sum_num, sum_den
    double precision :: f_old, Q_old(max_pipes)
    integer :: i, j, iter, max_iter, converged
    double precision :: tol
    
    ! Total results
    double precision :: dP_total, hf_total, pump_power
    
    pi = 3.14159265358979d0
    g = 9.81d0
    max_iter = 500
    tol = 1.0d-8
    
    ! Read inputs
    read(*,*) net_type    ! 1=series, 2=parallel
    read(*,*) n_pipes     ! Number of pipes
    read(*,*) rho         ! Density [kg/m3]
    read(*,*) mu          ! Dynamic viscosity [Pa.s]
    read(*,*) Q_total     ! Total flow rate [m3/s]
    read(*,*) dZ          ! Elevation difference [m]
    
    nu = mu / rho
    
    do i = 1, n_pipes
        read(*,*) L(i), D(i), eps(i)
        ! D and eps are in meters already (converted by PHP)
        Ac(i) = pi * D(i)**2 / 4.0d0
        rel_rough(i) = eps(i) / D(i)
    end do
    
    ! ============================================
    ! HEADER
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   PIPE NETWORK ANALYSIS'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ---'
    if (net_type == 1) then
        write(*,'(A)') '  Network Type            = SERIES'
    else
        write(*,'(A)') '  Network Type            = PARALLEL'
    end if
    write(*,'(A,I2)')       '  Number of pipes         = ', n_pipes
    write(*,'(A,F12.2,A)')  '  Fluid density ρ         = ', rho, ' kg/m³'
    write(*,'(A,ES12.4,A)') '  Dynamic viscosity μ     = ', mu, ' Pa·s'
    write(*,'(A,ES12.4,A)') '  Kinematic viscosity ν   = ', nu, ' m²/s'
    write(*,'(A,ES12.4,A)') '  Total flow rate Q       = ', Q_total, ' m³/s'
    write(*,'(A,F12.4,A)')  '  Total flow rate         = ', Q_total*1000.0d0*60.0d0, ' L/min'
    write(*,'(A,F12.4,A)')  '  Elevation difference    = ', dZ, ' m'
    write(*,*)
    
    write(*,'(A)') '--- PIPE DATA ---'
    write(*,'(A)') '  Pipe   L [m]     D [mm]    ε [mm]     ε/D        A [m²]'
    write(*,'(A)') '  ----   ------    ------    ------     ------     ------'
    do i = 1, n_pipes
        write(*,'(I5,2X,F8.2,2X,F8.2,2X,F8.4,3X,ES10.3,2X,ES10.4)') &
            i, L(i), D(i)*1000, eps(i)*1000, rel_rough(i), Ac(i)
    end do
    write(*,*)
    
    ! ============================================
    ! SERIES NETWORK
    ! ============================================
    if (net_type == 1) then
        write(*,'(A)') '--- SERIES ANALYSIS ---'
        write(*,'(A)') '  Same flow rate Q in all pipes.'
        write(*,*)
        
        dP_total = 0.0d0
        hf_total = 0.0d0
        
        do i = 1, n_pipes
            Q(i) = Q_total
            V(i) = Q(i) / Ac(i)
            Re(i) = rho * V(i) * D(i) / mu
            
            ! Friction factor
            call compute_friction(Re(i), rel_rough(i), f(i))
            
            ! Pressure drop (Darcy-Weisbach)
            dP(i) = f(i) * (L(i) / D(i)) * rho * V(i)**2 / 2.0d0
            hf(i) = dP(i) / (rho * g)
            
            dP_total = dP_total + dP(i)
            hf_total = hf_total + hf(i)
        end do
        
        ! Add elevation head
        dP_total = dP_total + rho * g * dZ
        hf_total = hf_total + dZ
        
        write(*,'(A)') '--- PER-PIPE RESULTS ---'
        write(*,'(A)') '  Pipe   Q [L/min]  V [m/s]    Re          f          ΔP [Pa]     h_f [m]'
        write(*,'(A)') '  ----   --------   ------     ----------  --------   ---------   -------'
        do i = 1, n_pipes
            write(*,'(I5,2X,F10.2,2X,F8.4,3X,ES10.4,2X,ES10.4,2X,F10.2,3X,F8.4)') &
                i, Q(i)*1000*60, V(i), Re(i), f(i), dP(i), hf(i)
        end do
        write(*,*)
        
        pump_power = dP_total * Q_total
        
        write(*,'(A)') '--- TOTAL SERIES RESULTS ---'
        write(*,'(A,F12.2,A)')  '  Total pressure drop     = ', dP_total, ' Pa'
        write(*,'(A,F12.4,A)')  '  Total pressure drop     = ', dP_total/1000, ' kPa'
        write(*,'(A,F12.4,A)')  '  Total head loss         = ', hf_total, ' m'
        write(*,'(A,F12.4,A)')  '  Elevation component     = ', rho*g*dZ, ' Pa'
        write(*,'(A,F12.4,A)')  '  Required pump power     = ', pump_power, ' W'
        write(*,*)
        
    ! ============================================
    ! PARALLEL NETWORK — HARDY-CROSS METHOD
    ! ============================================
    else
        write(*,'(A)') '--- PARALLEL ANALYSIS (HARDY-CROSS) ---'
        write(*,'(A)') '  Same pressure drop across all branches.'
        write(*,'(A)') '  Iterating to find flow distribution...'
        write(*,*)
        
        ! Initial guess: distribute proportional to D^2.5
        sum_num = 0.0d0
        do i = 1, n_pipes
            sum_num = sum_num + D(i)**2.5d0 / L(i)**0.5d0
        end do
        do i = 1, n_pipes
            Q(i) = Q_total * (D(i)**2.5d0 / L(i)**0.5d0) / sum_num
        end do
        
        write(*,'(A)') '--- CONVERGENCE HISTORY ---'
        write(*,'(A)') '  Iter    ΔQ [m³/s]      Max |Q change|   Σ Q [m³/s]'
        write(*,'(A)') '  ----    ----------      ---------------  ----------'
        
        converged = 0
        
        do iter = 1, max_iter
            ! Save old values
            do i = 1, n_pipes
                Q_old(i) = Q(i)
            end do
            
            ! Compute friction factors and head losses
            do i = 1, n_pipes
                V(i) = Q(i) / Ac(i)
                Re(i) = rho * abs(V(i)) * D(i) / mu
                if (Re(i) < 1.0d0) Re(i) = 1.0d0
                call compute_friction(Re(i), rel_rough(i), f(i))
            end do
            
            ! Hardy-Cross correction for parallel network
            ! For equal head loss: h_f1 = h_f2 = ... = h_fn
            ! Use head loss formulation: h_f = (8 f L Q|Q|) / (π² g D⁵)
            ! Correction: ΔQ = -Σ(R_i Q_i |Q_i|) / Σ(2 R_i |Q_i|)
            ! where R_i = 8 f_i L_i / (π² g D_i⁵)
            
            ! Actually for parallel, we equalize head losses
            ! Approach: adjust Q's so all have same head loss
            ! Using the Newton-like correction:
            
            ! Compute head losses
            do i = 1, n_pipes
                hf(i) = f(i) * L(i) * V(i) * abs(V(i)) / (2.0d0 * g * D(i))
            end do
            
            ! Target: average head loss
            sum_num = 0.0d0
            do i = 1, n_pipes
                sum_num = sum_num + hf(i)
            end do
            sum_num = sum_num / dble(n_pipes)  ! average h_f
            
            ! Adjust each pipe's flow to equalize head loss
            ! h_f = f L V² / (2gD) = f L Q² / (2g D A²)
            ! dh_f/dQ = f L 2Q / (2g D A²) = 2 h_f / Q (approximately)
            ! ΔQ_i = -(h_f_i - h_f_avg) / (2 h_f_i / Q_i)
            
            do i = 1, n_pipes
                if (abs(hf(i)) > 1.0d-15 .and. abs(Q(i)) > 1.0d-15) then
                    dQ = -(hf(i) - sum_num) / (2.0d0 * hf(i) / Q(i))
                    Q(i) = Q(i) + dQ * 0.5d0  ! relaxation factor
                end if
            end do
            
            ! Enforce mass conservation: Σ Q_i = Q_total
            sum_den = 0.0d0
            do i = 1, n_pipes
                sum_den = sum_den + Q(i)
            end do
            do i = 1, n_pipes
                Q(i) = Q(i) * Q_total / sum_den
            end do
            
            ! Check convergence
            dQ = 0.0d0
            do i = 1, n_pipes
                if (abs(Q(i) - Q_old(i)) > dQ) dQ = abs(Q(i) - Q_old(i))
            end do
            
            if (mod(iter, 10) == 1 .or. iter <= 5 .or. dQ < tol) then
                sum_den = 0.0d0
                do i = 1, n_pipes
                    sum_den = sum_den + Q(i)
                end do
                write(*,'(I5,3X,ES14.6,3X,ES14.6,3X,ES12.4)') iter, dQ, dQ, sum_den
            end if
            
            if (dQ < tol) then
                converged = 1
                write(*,*)
                write(*,'(A,I4,A)') '  ✓ Converged after ', iter, ' iterations'
                exit
            end if
        end do
        
        if (converged == 0) then
            write(*,*)
            write(*,'(A)') '  ⚠ WARNING: Did not converge within max iterations'
        end if
        write(*,*)
        
        ! Final results
        do i = 1, n_pipes
            V(i) = Q(i) / Ac(i)
            Re(i) = rho * abs(V(i)) * D(i) / mu
            call compute_friction(Re(i), rel_rough(i), f(i))
            dP(i) = f(i) * (L(i) / D(i)) * rho * V(i)**2 / 2.0d0
            hf(i) = dP(i) / (rho * g)
        end do
        
        write(*,'(A)') '--- PER-PIPE RESULTS ---'
        write(*,'(A)') '  Pipe   Q [L/min]  V [m/s]    Re          f          ΔP [Pa]     h_f [m]     Q/Q_t [%]'
        write(*,'(A)') '  ----   --------   ------     ----------  --------   ---------   -------     ---------'
        do i = 1, n_pipes
            write(*,'(I5,2X,F10.2,2X,F8.4,3X,ES10.4,2X,ES10.4,2X,F10.2,3X,F8.4,3X,F8.2)') &
                i, Q(i)*1000*60, V(i), Re(i), f(i), dP(i), hf(i), (Q(i)/Q_total)*100
        end do
        write(*,*)
        
        ! Use first pipe's head loss as reference (all should be equal)
        dP_total = dP(1)
        hf_total = hf(1)
        pump_power = dP_total * Q_total
        
        write(*,'(A)') '--- TOTAL PARALLEL RESULTS ---'
        write(*,'(A,F12.2,A)')  '  Common pressure drop    = ', dP_total, ' Pa'
        write(*,'(A,F12.4,A)')  '  Common pressure drop    = ', dP_total/1000, ' kPa'
        write(*,'(A,F12.4,A)')  '  Common head loss        = ', hf_total, ' m'
        write(*,'(A,F12.4,A)')  '  Required pump power     = ', pump_power, ' W'
        write(*,'(A,ES12.4,A)') '  Total flow rate check   = ', Q_total, ' m³/s'
        write(*,*)
    end if
    
    ! ============================================
    ! CHART DATA
    ! ============================================
    write(*,'(A)') '--- CHART_DATA ---'
    do i = 1, n_pipes
        write(*,'(I2,A,F12.6,A,F12.4,A,F12.2,A,F12.6)') &
            i, ',', Q(i)*1000*60, ',', V(i), ',', dP(i), ',', hf(i)
    end do
    write(*,'(A)') '--- END_CHART_DATA ---'
    
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  Darcy-Weisbach: ΔP = f(L/D)(ρV²/2)'
    write(*,'(A)') '  Colebrook-White: 1/√f = -2log₁₀(ε/D/3.7 + 2.51/(Re√f))'
    write(*,'(A)') '  Series: Q₁ = Q₂ = ... = Q_n; ΔP_total = Σ ΔP_i'
    write(*,'(A)') '  Parallel: ΔP₁ = ΔP₂ = ... = ΔP_n; Q_total = Σ Q_i'
    write(*,'(A)') '  Hardy-Cross: ΔQ = -Σ(RQ|Q|) / Σ(2R|Q|)'
    
contains
    
    subroutine compute_friction(Re_val, eD, f_out)
        double precision, intent(in) :: Re_val, eD
        double precision, intent(out) :: f_out
        double precision :: f_old_loc
        integer :: k
        
        if (Re_val < 2300.0d0) then
            f_out = 64.0d0 / Re_val
        else
            ! Swamee-Jain initial guess
            f_out = 0.25d0 / (log10(eD/3.7d0 + 5.74d0/Re_val**0.9d0))**2
            ! Colebrook-White iteration
            do k = 1, 100
                f_old_loc = f_out
                f_out = 1.0d0 / (-2.0d0 * log10(eD/3.7d0 + 2.51d0/(Re_val*sqrt(f_old_loc))))**2
                if (abs(f_out - f_old_loc) / f_out < 1.0d-10) exit
            end do
        end if
    end subroutine
    
end program pipe_networks
