! ============================================================================
! ThermoFluidCalc — Bernoulli & Extended Energy Equation Solver
! ============================================================================
program bernoulli_energy
    implicit none
    
    ! Input variables
    integer :: mode ! 1 = Simple Bernoulli, 2 = Extended Energy
    integer :: solve_for ! 1=P2, 2=V2, 3=z2, 4=Q, 5=h_p, 6=h_L
    double precision :: rho ! Fluid density [kg/m³]
    double precision :: P1, V1, z1 ! Point 1: Pressure [Pa], Velocity [m/s], Elevation [m]
    double precision :: P2, V2, z2 ! Point 2: Pressure [Pa], Velocity [m/s], Elevation [m]
    double precision :: D1, D2 ! Pipe diameters [m]
    double precision :: Q ! Volumetric flow rate [m³/s]
    double precision :: h_p, h_t, h_L ! Pump head [m], Turbine head [m], Head loss [m]
    
    ! Constants
    double precision, parameter :: g = 9.81d0
    double precision, parameter :: pi = 3.141592653589793d0
    
    ! Intermediate and output variables
    double precision :: hp1, hv1, hz1, H1
    double precision :: hp2, hv2, hz2, H2
    double precision :: C1, C2, RHS, Q2, V2_sq
    double precision :: W_p, W_t ! Power [W]
    character(len=50) :: solved_var_name
    
    ! Read input from stdin
    read(*,*) mode
    read(*,*) solve_for
    read(*,*) rho
    read(*,*) P1
    read(*,*) V1
    read(*,*) z1
    read(*,*) P2
    read(*,*) V2
    read(*,*) z2
    read(*,*) D1
    read(*,*) D2
    read(*,*) Q
    read(*,*) h_p
    read(*,*) h_t
    read(*,*) h_L
    
    ! Override simple Bernoulli head parameters to zero if mode=1
    if (mode == 1) then
        h_p = 0.0d0
        h_t = 0.0d0
        h_L = 0.0d0
    end if
    
    ! Calculate velocity from flow rate and diameters if available
    ! Point 1
    if (Q > 0.0d0 .and. D1 > 0.0d0) then
        V1 = 4.0d0 * Q / (pi * D1**2)
    else if (V1 > 0.0d0 .and. D1 > 0.0d0 .and. Q <= 0.0d0) then
        Q = V1 * pi * D1**2 / 4.0d0
    end if
    
    ! Point 2 (only if we are not solving for V2 or Q)
    if (solve_for /= 2 .and. solve_for /= 4) then
        if (Q > 0.0d0 .and. D2 > 0.0d0) then
            V2 = 4.0d0 * Q / (pi * D2**2)
        else if (V2 > 0.0d0 .and. D2 > 0.0d0 .and. Q <= 0.0d0) then
            Q = V2 * pi * D2**2 / 4.0d0
        end if
    end if
    
    ! Perform solving based on target variable
    select case (solve_for)
        case (1) ! Solve for P2
            solved_var_name = "Downstream Pressure (P2)"
            P2 = P1 + rho * g * ((V1**2 - V2**2) / (2.0d0 * g) + (z1 - z2) + h_p - h_t - h_L)
            
        case (2) ! Solve for V2
            solved_var_name = "Downstream Velocity (V2)"
            V2_sq = V1**2 + 2.0d0 * g * ((P1 - P2) / (rho * g) + (z1 - z2) + h_p - h_t - h_L)
            if (V2_sq < 0.0d0) then
                write(*,*) "ERROR: No physical solution for V2 (energy state is impossible)"
                stop
            end if
            V2 = sqrt(V2_sq)
            if (D2 > 0.0d0) then
                Q = V2 * pi * D2**2 / 4.0d0
            end if
            
        case (3) ! Solve for z2
            solved_var_name = "Downstream Elevation (z2)"
            z2 = z1 + (P1 - P2) / (rho * g) + (V1**2 - V2**2) / (2.0d0 * g) + h_p - h_t - h_L
            
        case (4) ! Solve for Q
            solved_var_name = "System Flow Rate (Q)"
            
            ! Coefficients for Q^2 relation
            C1 = 0.0d0
            C2 = 0.0d0
            if (D1 > 0.0d0) C1 = 8.0d0 / (g * pi**2 * D1**4)
            if (D2 > 0.0d0) C2 = 8.0d0 / (g * pi**2 * D2**4)
            
            if (D1 > 0.0d0 .and. D2 > 0.0d0) then
                if (abs(D1 - D2) < 1.0d-7) then
                    write(*,*) "ERROR: Cannot solve for Q when diameters are equal (V1 = V2)"
                    stop
                end if
                RHS = (P2 - P1) / (rho * g) + (z2 - z1) + h_t + h_L - h_p
                Q2 = RHS / (C1 - C2)
                if (Q2 < 0.0d0) then
                    write(*,*) "ERROR: No physical solution for Q (check energy terms direction)"
                    stop
                end if
                Q = sqrt(Q2)
                V1 = 4.0d0 * Q / (pi * D1**2)
                V2 = 4.0d0 * Q / (pi * D2**2)
            else if (D1 > 0.0d0 .and. D2 <= 0.0d0) then
                ! V2 is fixed, solve for Q
                RHS = (P2 - P1) / (rho * g) + V2**2 / (2.0d0 * g) + (z2 - z1) + h_t + h_L - h_p
                Q2 = RHS / C1
                if (Q2 < 0.0d0) then
                    write(*,*) "ERROR: No physical solution for Q"
                    stop
                end if
                Q = sqrt(Q2)
                V1 = 4.0d0 * Q / (pi * D1**2)
            else if (D1 <= 0.0d0 .and. D2 > 0.0d0) then
                ! V1 is fixed, solve for Q
                RHS = (P1 - P2) / (rho * g) + V1**2 / (2.0d0 * g) + (z1 - z2) + h_p - h_t - h_L
                Q2 = RHS / C2
                if (Q2 < 0.0d0) then
                    write(*,*) "ERROR: No physical solution for Q"
                    stop
                end if
                Q = sqrt(Q2)
                V2 = 4.0d0 * Q / (pi * D2**2)
            else
                write(*,*) "ERROR: Cannot solve for Q without at least one pipe diameter"
                stop
            end if
            
        case (5) ! Solve for h_p
            solved_var_name = "Required Pump Head (h_p)"
            h_p = (P2 - P1) / (rho * g) + (V2**2 - V1**2) / (2.0d0 * g) + (z2 - z1) + h_t + h_L
            
        case (6) ! Solve for h_L
            solved_var_name = "Head Loss (h_L)"
            h_L = (P1 - P2) / (rho * g) + (V1**2 - V2**2) / (2.0d0 * g) + (z1 - z2) + h_p - h_t
            if (h_L < 0.0d0) then
                write(*,*) "ERROR: Calculated Head Loss is negative (violates thermodynamics, check inputs)"
                stop
            end if
            
        case default
            write(*,*) "ERROR: Invalid solve target"
            stop
    end select
    
    ! Calculate final component heads
    hp1 = P1 / (rho * g)
    hv1 = V1**2 / (2.0d0 * g)
    hz1 = z1
    H1 = hp1 + hv1 + hz1
    
    hp2 = P2 / (rho * g)
    hv2 = V2**2 / (2.0d0 * g)
    hz2 = z2
    H2 = hp2 + hv2 + hz2
    
    ! Calculate power values
    W_p = 0.0d0
    W_t = 0.0d0
    if (Q > 0.0d0) then
        if (h_p > 0.0d0) W_p = rho * g * Q * h_p
        if (h_t > 0.0d0) W_t = rho * g * Q * h_t
    end if
    
    ! Output details
    write(*, '(A, A)') "Solved Variable = ", trim(solved_var_name)
    write(*, '(A, F14.4)') "Fluid Density (rho) = ", rho
    write(*, '(A, F14.4)') "Upstream Pressure (P1) = ", P1 / 1000.0d0 ! Pa to kPa
    write(*, '(A, F14.4)') "Upstream Velocity (V1) = ", V1
    write(*, '(A, F14.4)') "Upstream Elevation (z1) = ", z1
    write(*, '(A, F14.4)') "Downstream Pressure (P2) = ", P2 / 1000.0d0 ! Pa to kPa
    write(*, '(A, F14.4)') "Downstream Velocity (V2) = ", V2
    write(*, '(A, F14.4)') "Downstream Elevation (z2) = ", z2
    write(*, '(A, F14.4)') "Flow Rate (Q_lmin) = ", Q * 1000.0d0 * 60.0d0 ! m3/s to L/min
    write(*, '(A, F14.4)') "Flow Rate (Q_m3h) = ", Q * 3600.0d0 ! m3/s to m3/h
    write(*, '(A, F14.4)') "Pipe 1 Diameter (D1) = ", D1 * 1000.0d0 ! m to mm
    write(*, '(A, F14.4)') "Pipe 2 Diameter (D2) = ", D2 * 1000.0d0 ! m to mm
    write(*, '(A, F14.4)') "Pump Head Added (h_p) = ", h_p
    write(*, '(A, F14.4)') "Turbine Head Extracted (h_t) = ", h_t
    write(*, '(A, F14.4)') "Head Loss (h_L) = ", h_L
    write(*, '(A, F14.4)') "Pressure Head 1 (hp1) = ", hp1
    write(*, '(A, F14.4)') "Velocity Head 1 (hv1) = ", hv1
    write(*, '(A, F14.4)') "Elevation Head 1 (hz1) = ", hz1
    write(*, '(A, F14.4)') "Total Head 1 (H1) = ", H1
    write(*, '(A, F14.4)') "Pressure Head 2 (hp2) = ", hp2
    write(*, '(A, F14.4)') "Velocity Head 2 (hv2) = ", hv2
    write(*, '(A, F14.4)') "Elevation Head 2 (hz2) = ", hz2
    write(*, '(A, F14.4)') "Total Head 2 (H2) = ", H2
    write(*, '(A, F14.4)') "Pump Power (W_p) = ", W_p / 1000.0d0 ! W to kW
    write(*, '(A, F14.4)') "Pump Power (HP) = ", W_p / 745.7d0 ! W to HP
    write(*, '(A, F14.4)') "Turbine Power (W_t) = ", W_t / 1000.0d0 ! W to kW
    write(*, '(A, F14.4)') "Turbine Power (HP) = ", W_t / 745.7d0 ! W to HP
    
end program bernoulli_energy
