"""radial_inflow.py — Topic 2.1 reference implementation
Course 01 · Module 02 · Topic 2.1 — Darcy's Law for Radial Flow
Well Productivity Programme (Abalt Solutions)

Single-phase (oil) inflow calculations in oilfield units:
  - linear Darcy flow (the sand-column experiment)
  - radial steady-state (SS) inflow
  - radial pseudo-steady-state (PSS) inflow
  - the logarithmic pressure-traverse profile P(r)

The 0.00708 field constant is rebuilt from first principles in
field_constant() so it can never be mis-remembered. The radial Darcy law in
pure cgs is q[cc/s] = 2*pi*k[D]*h[cm]*dP[atm]/(mu[cp]*ln(re/rw)); converting
the field units mD, ft, psi and the rate cc/s -> bbl/day gives 0.007082
(= 1/141.2). See README for the term-by-term chain and a note on the brief's
shorthand derivation.

__main__ reproduces the brief's North Sea / KRM-4 worked example:
    SS  Q = 1,472 STB/d   (ln re/rw = 8.22, denom 14.39)
    PSS Q = 1,620 STB/d today, 1,019 @ 4,200, 371 @ 3,500 (denom 13.07)
Run:  python3 radial_inflow.py
"""

import math

# ---------------------------------------------------------------------------
# Field constant
# ---------------------------------------------------------------------------

def field_constant():
    """Rebuild the 0.00708 radial-inflow constant from unit conversions.

        Q[STB/d] = const * k[mD] * h[ft] * dP[psi] / (mu[cp]*B[RB/STB]*denom)

    Radial Darcy law in cgs:  q[cc/s] = 2*pi*k[D]*h[cm]*dP[atm]/(mu[cp]*denom)
    Substitute one field unit of each input (k=1 mD, h=1 ft, dP=1 psi, mu=1 cp)
    and convert the resulting cc/s rate to bbl/day:
        k : 1 mD            = 1/1000 D
        h : 1 ft            = 30.48 cm
        dP: 1 psi           = 1/14.696 atm
        q : cc/s -> bbl/day = * 86,400 s/day / 158,987 cc/bbl
    """
    k_D = 1.0 / 1000.0
    h_cm = 30.48
    dP_atm = 1.0 / 14.696
    q_ccs = 2.0 * math.pi * k_D * h_cm * dP_atm
    return q_ccs * 86_400.0 / 158_987.0


# Canonical constant used on the WPP decks/brief (= 1/141.2).
C_FIELD = 0.00708


# ---------------------------------------------------------------------------
# Linear Darcy flow (the sand column)
# ---------------------------------------------------------------------------

def linear_darcy_rate(k_md, area_ft2, dP_psi, mu_cp, length_ft, B=1.0):
    """Linear (1-D) Darcy flow through a constant-area sand column.

        Q[STB/d] = 0.00708 * k * A * dP / (mu * L * B)

    Note: for a linear element the geometric "denominator" is simply the
    length L; A is the cross-sectional area. With B = 1.0 the result is a
    reservoir-volume rate; pass B to convert to stock-tank.
    Returned in STB/d (reservoir bbl/d if B = 1.0). Uses the same 0.00708
    bookkeeping constant as the radial form.
    """
    # The radial constant 0.00708 already folds 2*pi; remove it for the
    # plane-linear analogue by dividing by 2*pi to recover the pure
    # Darcy unit conversion 1.127e-3.
    c_lin = C_FIELD / (2.0 * math.pi)   # = 1.127e-3
    return c_lin * k_md * area_ft2 * dP_psi / (mu_cp * length_ft * B)


def darcy_velocity(k_md, mu_cp, dPdx_psi_ft):
    """Darcy (superficial) velocity u = -(k/mu) dP/dx, ft/day.

    Convenience helper, returns magnitude in ft/day for a given gradient.
    """
    # 1.127e-3 converts mD, psi, cp, ft -> bbl/(d.ft^2); *5.615 ft^3/bbl
    c_lin = C_FIELD / (2.0 * math.pi)
    return c_lin * k_md / mu_cp * dPdx_psi_ft * 5.615


# ---------------------------------------------------------------------------
# Radial inflow
# ---------------------------------------------------------------------------

def ln_re_rw(re_ft, rw_ft):
    """Natural log of the radius ratio (the geometric heart of radial flow)."""
    return math.log(re_ft / rw_ft)


def radial_ss_rate(k_md, h_ft, P_e, Pwf, mu_cp, B, re_ft, rw_ft, S=0.0):
    """Steady-state (constant-pressure boundary) radial oil inflow.

        Q = 0.00708 k h (Pe - Pwf) / [mu B (ln(re/rw) + S)]

    Reference pressure is the boundary value Pe (aquifer / injection support).
    Returns Q in STB/d.
    """
    denom = mu_cp * B * (ln_re_rw(re_ft, rw_ft) + S)
    return C_FIELD * k_md * h_ft * (P_e - Pwf) / denom


def radial_pss_rate(k_md, h_ft, P_bar, Pwf, mu_cp, B, re_ft, rw_ft, S=0.0):
    """Pseudo-steady-state (no-flow boundary) radial oil inflow.

        Q = 0.00708 k h (P_bar - Pwf) / [mu B (ln(re/rw) - 0.75 + S)]

    Reference pressure is the volumetric average P_bar. The -0.75 term comes
    from averaging the logarithmic profile over the drainage area
    (ln(re/rw) - 0.75 = ln(0.472 re/rw)). Returns Q in STB/d.
    """
    denom = mu_cp * B * (ln_re_rw(re_ft, rw_ft) - 0.75 + S)
    return C_FIELD * k_md * h_ft * (P_bar - Pwf) / denom


def pressure_profile(r_ft, Pwf, P_ref, rw_ft, re_ft):
    """Logarithmic radial pressure traverse P(r).

        P(r) = Pwf + (P_ref - Pwf) * ln(r/rw) / ln(re/rw)

    This is the integrated shell equation P(r) = Pwf + [q mu /(2 pi k h)] ln(r/rw)
    written with the drawdown (P_ref - Pwf) substituted for the leading group,
    so the endpoints are exact: P(rw) = Pwf and P(re) = P_ref.
    Pass P_ref = Pe for the SS profile (or P_bar for an approximate PSS funnel).
    """
    return Pwf + (P_ref - Pwf) * math.log(r_ft / rw_ft) / math.log(re_ft / rw_ft)


def drawdown_fraction_within(r_ft, rw_ft, re_ft):
    """Fraction of total drawdown consumed between rw and r_ft.

        f(r) = ln(r/rw) / ln(re/rw)
    """
    return math.log(r_ft / rw_ft) / math.log(re_ft / rw_ft)


def shell_area(r_ft, h_ft):
    """Cylindrical shell area A(r) = 2 pi r h (ft^2)."""
    return 2.0 * math.pi * r_ft * h_ft


# ---------------------------------------------------------------------------
# __main__ — reproduce the brief's North Sea / KRM-4 worked example
# ---------------------------------------------------------------------------

if __name__ == "__main__":
    # KRM-4 canonical inputs (TECH-BRIEF-m02 §1.2 / Topic 2.1 worked example)
    k, h = 18.0, 95.0          # mD, ft
    mu, B = 1.4, 1.25          # cp, RB/STB  (use Bo(P_bar)=1.25, NOT Bob)
    re, rw = 1320.0, 0.354     # ft (160-acre spacing), ft
    Pwf = 3100.0               # psia (separator back-pressure)
    PR = 4850.0                # psia (boundary / current average)

    print("=== Topic 2.1 radial_inflow.py — KRM-4 worked example ===\n")

    print("Field constant 0.00708 rebuilt from units : %.6f  (target 0.00708)"
          % field_constant())
    print("Reciprocal 1/const                         : %.2f   (target 141.2)\n"
          % (1.0 / field_constant()))

    ratio = re / rw
    ln = ln_re_rw(re, rw)
    print("Step 1  radius ratio re/rw    = %.0f" % ratio)
    print("        ln(re/rw)             = %.3f   (deck 8.22)" % ln)

    num = C_FIELD * k * h * (PR - Pwf)
    print("Step 2  numerator 0.00708 k h dP = %.0f   (deck 21,187)" % num)

    denom_ss = mu * B * ln
    Q_ss = radial_ss_rate(k, h, PR, Pwf, mu, B, re, rw, S=0.0)
    print("Step 3  SS denom mu B ln      = %.2f   (deck 14.39)" % denom_ss)
    print("        SS  Q                 = %.0f STB/d   (deck 1,472)\n" % Q_ss)

    denom_pss = mu * B * (ln - 0.75)
    print("PSS denom mu B (ln-0.75)      = %.2f   (deck 13.07)" % denom_pss)
    for Pbar, target in ((4850, 1620), (4200, 1019), (3500, 371)):
        q = radial_pss_rate(k, h, Pbar, Pwf, mu, B, re, rw, S=0.0)
        print("PSS Q @ P_bar = %4d psia     = %.0f STB/d   (deck %d)"
              % (Pbar, q, target))

    print("\nShell areas (geometry):")
    print("  A(re) = 2 pi re h           = %.0f ft^2   (deck ~790,000)"
          % shell_area(re, h))
    print("  A(rw) = 2 pi rw h           = %.0f ft^2   (deck 211)"
          % shell_area(rw, h))

    print("\nPressure funnel (SS profile, P_ref = 4,850):")
    for r in (rw, 10.0, 22.0, re):
        P = pressure_profile(r, Pwf, PR, rw, re)
        f = drawdown_fraction_within(r, rw, re)
        print("  r = %7.3f ft : P = %6.1f psia   drawdown consumed = %5.1f%%"
              % (r, P, f * 100.0))
    print("  -> ~40%% of drawdown within 10 ft; ~50%% by 22 ft (deck)")

    print("\nLinear Darcy check (sand column, B=1): finite, positive ->",
          "%.1f bbl/d" % linear_darcy_rate(1000.0, 1.0, 100.0, 1.0, 1.0))