"""
skin_basics.py — Topic 3.1 reference implementation (Course 01 · Module 03 · WPP).

Skin-factor fundamentals: the skin pressure drop, Hawkins' formula, effective
(apparent) wellbore radius, the 7-approximation productivity ratio / flow
efficiency, and ideal-vs-damaged radial pressure profiles.

All equations and canonical numbers follow the Module 03 technical brief
(TECH-BRIEF-m03.md) and the BP Completion Design Manual field-unit forms:

    q  = 0.00708 k h (pR - pwf) / [ mu B (ln(0.472 re/rw) + S') ]   [stb/d]
    dp_skin = 141.2 q mu B S / (k h)                                 [psi]
    Hawkins: S = (k/ks - 1) ln(rs/rw)
    rw' = rw exp(-S)
    FE  = 7 / (7 + S)        (ln(0.472 re/rw) ~ 7 screening approximation)

Running `python3 skin_basics.py` reproduces the GK-22 worked chain used in
Lecture 3.1 (74.18 psi per skin unit, dp_skin = 1,038 psi = 61% of drawdown,
J 1.380 -> 0.460 stb/d/psi, 2,346 vs 782 stb/d, FE = 0.333).

(c) 2026 Abalt Solutions Ltd - training reference code.
"""

import math

# ---------------------------------------------------------------- constants
C_FLOW = 0.00708          # field-unit radial-flow constant
C_DP = 141.2              # = 1 / 0.00708, skin pressure-drop constant
LN7 = 7.0                 # ln(0.472 re/rw) screening approximation


# ------------------------------------------------------------ core formulae
def ln_term(re_ft, rw_ft):
    """ln(0.472 re/rw) — pseudo-steady-state log term (GK-22: 7.709)."""
    return math.log(0.472 * re_ft / rw_ft)


def psi_per_skin_unit(q_stbd, mu_cp, B_rb_stb, k_md, h_ft):
    """Cost of ONE skin unit in psi: 141.2 q mu B / (k h)."""
    return C_DP * q_stbd * mu_cp * B_rb_stb / (k_md * h_ft)


def dp_skin(q_stbd, mu_cp, B_rb_stb, k_md, h_ft, S):
    """Skin pressure drop, psi: dp_skin = 141.2 q mu B S / (k h)."""
    return psi_per_skin_unit(q_stbd, mu_cp, B_rb_stb, k_md, h_ft) * S


def skin_from_dp(dp_psi, q_stbd, mu_cp, B_rb_stb, k_md, h_ft):
    """Inverse: S = 0.00708 k h dp_skin / (q mu B)."""
    return C_FLOW * k_md * h_ft * dp_psi / (q_stbd * mu_cp * B_rb_stb)


def hawkins(k_md, ks_md, rs_ft, rw_ft):
    """Hawkins (1956): S = (k/ks - 1) ln(rs/rw)."""
    return (k_md / ks_md - 1.0) * math.log(rs_ft / rw_ft)


def hawkins_rs(S, k_md, ks_md, rw_ft):
    """Inverse Hawkins: damage radius rs = rw exp[S / (k/ks - 1)]."""
    return rw_ft * math.exp(S / (k_md / ks_md - 1.0))


def apparent_rw(rw_ft, S):
    """Effective (apparent) wellbore radius rw' = rw exp(-S)."""
    return rw_ft * math.exp(-S)


def prats_skin(xf_ft, rw_ft):
    """Prats fracture skin: Ss = ln(2 rw / Xf)."""
    return math.log(2.0 * rw_ft / xf_ft)


def s_min(re_ft, rw_ft):
    """Minimum (floor) skin: Smin = ln(rw / 0.472 re)."""
    return math.log(rw_ft / (0.472 * re_ft))


def productivity_index(k_md, h_ft, mu_cp, B_rb_stb, re_ft, rw_ft, S=0.0):
    """J = 0.00708 k h / [mu B (ln(0.472 re/rw) + S)]  [stb/d/psi]."""
    return C_FLOW * k_md * h_ft / (mu_cp * B_rb_stb * (ln_term(re_ft, rw_ft) + S))


def j_ratio_7(S):
    """Screening productivity ratio Jdamaged/Jideal = 7/(7+S) (= Standing FE)."""
    return LN7 / (LN7 + S)


def flow_efficiency(S, ln_val=LN7):
    """Standing flow efficiency FE = ln/(ln+S); default uses the 7-approx."""
    return ln_val / (ln_val + S)


# ------------------------------------------------------- pressure profiles
def pressure_profiles(radii_ft, pwf_psi, q_stbd, mu_cp, B_rb_stb, k_md, h_ft,
                      ks_md, rs_ft, rw_ft):
    """
    Ideal and damaged radial pressure traverses p(r), anchored at the measured
    flowing pressure pwf (damaged-well sandface).

    Per-log-cycle slope m = 141.2 q mu B / (k h); inside the damaged annulus
    the slope steepens by k/ks. The ideal profile is offset upward by dp_skin
    so the two curves coincide beyond rs (far field identical).

    Returns list of (r, ln(r/rw), p_damaged, p_ideal) tuples.
    """
    m = psi_per_skin_unit(q_stbd, mu_cp, B_rb_stb, k_md, h_ft)
    ratio = k_md / ks_md
    step = m * (ratio - 1.0) * math.log(rs_ft / rw_ft)   # = dp_skin (Hawkins)
    rows = []
    for r in radii_ft:
        if r <= rs_ft:
            p_dam = pwf_psi + m * ratio * math.log(r / rw_ft)
        else:
            p_dam = (pwf_psi + m * ratio * math.log(rs_ft / rw_ft)
                     + m * math.log(r / rs_ft))
        p_ideal = pwf_psi + step + m * math.log(r / rw_ft)
        rows.append((r, math.log(r / rw_ft), p_dam, p_ideal))
    return rows


# ------------------------------------------------------------------- demo
if __name__ == '__main__':
    # GK-22 canonical dataset (locked, TECH-BRIEF-m03 §1a)
    k, h, mu, B = 85.0, 42.0, 1.8, 1.32
    pR, pwf = 4200.0, 2500.0
    rw, re = 0.35, 1650.0
    q, S = 782.0, 14.0
    ks_ratio, rs = 0.145, 3.77

    dd = pR - pwf
    m = psi_per_skin_unit(q, mu, B, k, h)
    dps = dp_skin(q, mu, B, k, h, S)
    j_id = productivity_index(k, h, mu, B, re, rw, S=0.0)
    j_7 = j_id * j_ratio_7(S)
    q_id = j_id * dd

    print('GK-22 skin pressure-drop chain (Lecture 3.1)')
    print('-' * 56)
    print(f'ln(0.472 re/rw)          = {ln_term(re, rw):6.3f}   (canon 7.709)')
    print(f'psi per skin unit        = {m:7.2f}  (canon 74.18)')
    print(f'dp_skin = 14 x 74.18     = {dps:7.0f}  psi (canon 1,038)')
    print(f'share of 1,700-psi dd    = {100 * dps / dd:6.1f}%  (canon 61%)')
    print(f'Darcy-work pressure      = {dd - dps:7.0f}  psi (canon 662)')
    print(f'J ideal (S=0)            = {j_id:7.3f}  stb/d/psi (canon 1.380)')
    print(f'J actual (7-approx)      = {j_7:7.3f}  stb/d/psi (canon 0.460)')
    print(f'q ideal at pwf 2,500     = {q_id:7.0f}  stb/d (canon 2,346)')
    print(f'q actual / lost          = {q:5.0f} / {q_id - q:5.0f} stb/d '
          f'({100 * (q_id - q) / q_id:.0f}% lost; canon 1,564 / 67%)')
    print(f'FE = 7/(7+S)             = {flow_efficiency(S):7.3f}  (canon 0.333)')
    print()
    print('Hawkins / effective-radius cross-checks')
    print('-' * 56)
    s_hk = hawkins(k, ks_ratio * k, rs, rw)
    print(f'Hawkins S (ks/k 0.145, rs 3.77 ft) = {s_hk:6.2f}  (canon +14.0)')
    print(f'inverse rs at S=+14                = {hawkins_rs(S, k, ks_ratio * k, rw):6.2f} ft (canon 3.77)')
    print(f"apparent rw' at S = +14            = {apparent_rw(rw, S):.2e} ft (strangled)")
    print(f"apparent rw' at S = -5             = {apparent_rw(rw, -5):6.1f} ft (canon 52)")
    print(f'Prats Ss (Xf 400 ft, rw 0.33)      = {prats_skin(400, 0.33):6.2f}  (Kessler-7 canon -6.41)')
    print(f'Smin (re 2,980 ft, rw 0.33)        = {s_min(2980, 0.33):6.2f}  (canon ~ -8.4 -> floor -7)')
    print()
    print('Pressure profile end-points (anchored at pwf 2,500)')
    print('-' * 56)
    radii = [0.35, 3.77, 1650.0]
    for r, lnx, p_dam, p_id in pressure_profiles(radii, pwf, q, mu, B, k, h,
                                                 ks_ratio * k, rs, rw):
        print(f'r = {r:7.2f} ft   ln(r/rw) = {lnx:5.3f}   '
              f'p_damaged = {p_dam:6.0f} psi   p_ideal = {p_id:6.0f} psi')
    gap = (pressure_profiles([rw], pwf, q, mu, B, k, h, ks_ratio * k, rs, rw)[0][3]
           - pwf)
    print(f'sandface gap p_ideal - p_damaged   = {gap:6.0f} psi (= dp_skin, canon 1,038)')