! ============================================================================
! ThermoFluidCalc — Water Hammer (Coup de bélier) Solver
! Reference: Wylie & Streeter, Fluid Transients in Systems; Joukowsky (1898)
! ============================================================================
program water_hammer
    implicit none
    
    ! Inputs
    double precision :: L ! Pipe length [m]
    double precision :: D ! Pipe inner diameter [mm]
    double precision :: e_thick ! Pipe wall thickness [mm]
    double precision :: E_pipe ! Young's Modulus of pipe material [GPa]
    double precision :: V0 ! Initial flow velocity [m/s]
    double precision :: Kf ! Bulk modulus of fluid [GPa]
    double precision :: rho ! Fluid density [kg/m³]
    double precision :: tc ! Valve closure time [s]
    integer :: profile_option ! 1 = Instantaneous, 2 = Linear, 3 = Parabolic
    double precision :: P0 ! Static pressure [kPa]
    double precision :: Pv ! Vapor pressure [kPa]
    
    ! Constants
    double precision, parameter :: g = 9.81d0
    double precision, parameter :: pi = 3.141592653589793d0
    
    ! Outputs / Derived values
    double precision :: a ! Wave speed [m/s]
    double precision :: t_crit ! Critical time 2L/a [s]
    double precision :: DP_inst ! Theoretical instantaneous pressure surge [kPa]
    double precision :: DP_grad ! Theoretical gradual pressure surge (rigid column) [kPa]
    double precision :: DP_max ! Absolute maximum pressure surge from simulation [kPa]
    double precision :: P_max ! Maximum absolute pressure [kPa]
    double precision :: P_min ! Minimum absolute pressure [kPa]
    double precision :: F_support ! Force on pipe supports [kN]
    double precision :: hoop_stress ! Stress in pipe wall [MPa]
    logical :: cavitation_occurred
    
    ! MOC Solver variables
    integer, parameter :: N = 20 ! Number of spatial reaches
    double precision :: dx ! Spatial step [m]
    double precision :: dt ! Time step [s]
    integer :: N_steps ! Total time steps
    double precision :: t_max ! Simulation duration [s]
    double precision :: f_fric ! Darcy friction factor
    double precision :: Re ! Reynolds number
    double precision :: mu_visc ! Dynamic viscosity [Pa-s]
    double precision :: A_area ! Pipe cross-sectional area [m²]
    double precision :: Q0 ! Initial volumetric flow rate [m³/s]
    double precision :: B ! Impedance parameter
    double precision :: R_fric ! Friction resistance factor
    
    ! Grid arrays
    double precision :: H(0:N) ! Hydraulic head at current time step [m]
    double precision :: Q(0:N) ! Volumetric flow rate at current time step [m³/s]
    double precision :: H_new(0:N)
    double precision :: Q_new(0:N)
    
    ! Time series tracking at the valve (node N)
    integer, parameter :: MAX_TS_POINTS = 500
    double precision :: ts_time(MAX_TS_POINTS)
    double precision :: ts_pres(MAX_TS_POINTS)
    integer :: ts_count, skip_step
    
    ! Temporary variables
    integer :: i, nt
    double precision :: t_val, tau, CP, CM, Cv, Q_N_temp, H_N_temp, Hv, H0_head, term, eps_over_D
    character(len=20) :: closure_type
    
    ! Read inputs from stdin
    read(*,*) L
    read(*,*) D
    read(*,*) e_thick
    read(*,*) E_pipe
    read(*,*) V0
    read(*,*) Kf
    read(*,*) rho
    read(*,*) tc
    read(*,*) profile_option
    read(*,*) P0
    read(*,*) Pv
    
    ! Defaults & safety bounds
    if (L <= 0.0d0) L = 100.0d0
    if (D <= 0.0d0) D = 100.0d0
    if (e_thick <= 0.0d0) e_thick = 5.0d0
    if (E_pipe <= 0.0d0) E_pipe = 200.0d0
    if (Kf <= 0.0d0) Kf = 2.2d0
    if (rho <= 0.0d0) rho = 1000.0d0
    if (P0 <= 0.0d0) P0 = 200.0d0
    if (Pv < 0.0d0) Pv = 2.34d0
    if (tc < 0.0d0) tc = 0.0d0
    
    ! ── 1. WAVE SPEED CALCULATIONS (Joukowsky) ─────────────────
    ! bulk modulus and Young's modulus converted to Pa (1e9)
    ! a = sqrt(Kf/rho) / sqrt(1 + (Kf * D)/(E * e))
    a = sqrt((Kf * 1.0d9) / rho) / sqrt(1.0d0 + (Kf * D) / (E_pipe * e_thick))
    
    ! Critical time t_crit = 2L/a
    t_crit = 2.0d0 * L / a
    
    ! Closure type
    if (tc < t_crit) then
        closure_type = "Instantaneous"
    else
        closure_type = "Gradual"
    end if
    
    ! Theoretical surges
    DP_inst = (rho * a * V0) / 1000.0d0 ! in kPa
    if (tc > 0.0d0) then
        DP_grad = (rho * L * V0 / tc) / 1000.0d0 ! rigid column in kPa
    else
        DP_grad = DP_inst
    end if
    
    ! ── 2. PREPARE METHOD OF CHARACTERISTICS (MOC) ─────────────
    dx = L / dble(N)
    dt = dx / a
    
    ! Simulate for 5 wave cycles: T = 4L/a = 2 * t_crit
    t_max = 5.0d0 * (4.0d0 * L / a)
    N_steps = nint(t_max / dt)
    if (N_steps < 100) N_steps = 100
    
    ! Pipe area
    A_area = pi * (D / 1000.0d0)**2 / 4.0d0
    Q0 = V0 * A_area
    
    ! Friction factor estimation (Haaland)
    mu_visc = 0.001d0 ! water viscosity fallback
    Re = rho * abs(V0) * (D / 1000.0d0) / mu_visc
    if (Re < 2300.0d0) then
        if (Re > 0.0d0) then
            f_fric = 64.0d0 / Re
        else
            f_fric = 0.02d0
        end if
    else
        eps_over_D = (0.05d0 / 1000.0d0) / (D / 1000.0d0)
        term = (eps_over_D / 3.7d0)**1.11d0 + 6.9d0 / Re
        f_fric = 1.0d0 / (-1.8d0 * log10(term))**2
    end if
    
    ! MOC Parameters
    B = a / (g * A_area)
    R_fric = f_fric * dx / (2.0d0 * g * (D / 1000.0d0) * A_area**2)
    
    ! Initial Steady State Grade Line
    ! Reservoir at node 0 (fixed static pressure P0)
    H0_head = P0 / (rho * g / 1000.0d0)
    Hv = Pv / (rho * g / 1000.0d0)
    
    do i = 0, N
        H(i) = H0_head - dble(i) * R_fric * Q0**2
        Q(i) = Q0
    end do
    
    ! Initialize tracking
    DP_max = 0.0d0
    P_max = P0
    P_min = P0
    cavitation_occurred = .false.
    
    ts_count = 0
    skip_step = max(1, N_steps / MAX_TS_POINTS)
    
    ! Store initial point
    ts_count = ts_count + 1
    ts_time(ts_count) = 0.0d0
    ts_pres(ts_count) = H(N) * (rho * g / 1000.0d0)
    
    ! ── 3. MOC TRANSIENT LOOP ──────────────────────────────────
    do nt = 1, N_steps
        t_val = dble(nt) * dt
        
        ! Determine valve opening tau
        if (profile_option == 1) then
            ! Instantaneous
            tau = 0.0d0
        elseif (profile_option == 2) then
            ! Linear
            if (t_val < tc) then
                tau = 1.0d0 - (t_val / tc)
            else
                tau = 0.0d0
            end if
        else
            ! Parabolic
            if (t_val < tc) then
                tau = (1.0d0 - (t_val / tc))**2
            else
                tau = 0.0d0
            end if
        end if
        
        ! Interior nodes i = 1 ... N-1
        do i = 1, N-1
            CP = H(i-1) + B * Q(i-1) - R_fric * Q(i-1) * abs(Q(i-1))
            CM = H(i+1) - B * Q(i+1) + R_fric * Q(i+1) * abs(Q(i+1))
            
            H_new(i) = (CP + CM) / 2.0d0
            Q_new(i) = (CP - CM) / (2.0d0 * B)
        end do
        
        ! Reservoir Node i = 0 (Constant head)
        H_new(0) = H0_head
        CM = H(1) - B * Q(1) + R_fric * Q(1) * abs(Q(1))
        Q_new(0) = (H_new(0) - CM) / B
        
        ! Valve Node i = N
        CP = H(N-1) + B * Q(N-1) - R_fric * Q(N-1) * abs(Q(N-1))
        Cv = (tau * Q0)**2 / H0_head
        
        if (Cv > 0.0d0) then
            ! Quadratic formula for valve flow Q_N
            Q_N_temp = (-Cv * B + sqrt((Cv * B)**2 + 4.0d0 * Cv * CP)) / 2.0d0
            H_N_temp = CP - B * Q_N_temp
        else
            Q_N_temp = 0.0d0
            H_N_temp = CP
        end if
        
        ! Cavitation check at valve
        if (H_N_temp < Hv) then
            H_N_temp = Hv
            if (tau > 0.0d0) then
                Q_N_temp = (CP - Hv) / B
            else
                Q_N_temp = 0.0d0
            end if
            cavitation_occurred = .true.
        end if
        
        H_new(N) = H_N_temp
        Q_new(N) = Q_N_temp
        
        ! Update grid values
        H = H_new
        Q = Q_new
        
        ! Track absolute pressure at the valve
        H_N_temp = H(N)
        H_N_temp = max(H_N_temp, Hv) ! Physical bound
        
        ! Convert to pressure
        t_val = H_N_temp * (rho * g / 1000.0d0)
        
        if (t_val > P_max) P_max = t_val
        if (t_val < P_min) P_min = t_val
        
        ! Time series store
        if (mod(nt, skip_step) == 0 .and. ts_count < MAX_TS_POINTS) then
            ts_count = ts_count + 1
            ts_time(ts_count) = dble(nt) * dt
            ts_pres(ts_count) = t_val
        end if
    end do
    
    ! Ensure last point is captured
    if (ts_count < MAX_TS_POINTS) then
        ts_count = ts_count + 1
        ts_time(ts_count) = dble(N_steps) * dt
        ts_pres(ts_count) = H(N) * (rho * g / 1000.0d0)
    end if
    
    ! Calculate pressure surge surge
    DP_max = P_max - P0
    if (DP_max < 0.0d0) DP_max = 0.0d0
    
    ! Force on pipe supports (F = A * DP_max)
    ! DP_max is in kPa, A_area in m² -> F in kN
    F_support = A_area * DP_max
    
    ! Stress in pipe wall (hoop stress)
    ! sigma = P * D / (2 * e_thick)
    ! P_max is in kPa -> P_max / 1000 is in MPa
    hoop_stress = (P_max / 1000.0d0) * D / (2.0d0 * e_thick)
    
    ! ── 4. OUTPUT RESULTS IN KEY-VALUE FORMAT ─────────────────
    write(*, '(A, F14.2)') "Wave Speed = ", a
    write(*, '(A, F14.4)') "Critical Time = ", t_crit
    write(*, '(A, A)') "Closure Type = ", trim(closure_type)
    write(*, '(A, F14.2)') "DP Inst = ", DP_inst
    write(*, '(A, F14.2)') "DP Grad = ", DP_grad
    write(*, '(A, F14.2)') "Pressure Surge = ", DP_max
    write(*, '(A, F14.2)') "Max Pressure = ", P_max
    write(*, '(A, F14.2)') "Min Pressure = ", P_min
    write(*, '(A, A)') "Cavitation Risk = ", merge("Yes", "No ", cavitation_occurred)
    write(*, '(A, F14.3)') "Support Force = ", F_support
    write(*, '(A, F14.2)') "Hoop Stress = ", hoop_stress
    
    ! Timeline data points for PHP parsing
    write(*, '(A)') "--- TIMELINE DATA ---"
    do i = 1, ts_count
        write(*, '(F10.5, A, F12.2)') ts_time(i), ",", ts_pres(i)
    end do
    
end program water_hammer