#!/usr/bin/env python3
"""
pseudoskin.py — Topic 3.4 reference implementation (Pseudo-Skin Components Sc, Sc″).

Well Productivity Programme · Course 01 · Module 03 · Topic 3.4
Implements the Brons–Marting partial-completion skin, the Cinco-Ley deviation
skin, the Jones–Watts h/hp damage multiplier, the combined completion skin,
and a Joshi (1988) horizontal-well equivalent skin used for the SYNTHETIC
Dunbar DB-17 vignette.

All formulas and constants follow TECH-BRIEF-m03.md (Topic 3.4 section), which
is extracted from the authoritative topic page C01_M03_T4. Known source-page
quirks (design-table Sc = +5.0 vs worked-example +3.08; Cinco-Ley table −2.1
vs worked −1.26; Jones–Watts J = 0.214 chain) are reproduced/documented here —
see README.md and the __main__ demo.

Run:  python3 pseudoskin.py
"""

import math

# ---------------------------------------------------------------- constants
LN_TERM_GK22 = 7.71          # ln(0.472 * 1650 / 0.35) — full-log GK-22 value
LN7 = 7.0                    # the "≈7" classroom approximation
GK22 = dict(k=85.0, h=42.0, mu=1.8, B=1.32, pbar=4200.0, pwf=2500.0,
            rw=0.35, re=1650.0, q=782.0, S_prime=14.0)


# ----------------------------------------------------- dimensionless groups
def penetration_ratio(hp, h):
    """b = hp / h (perforated fraction of net pay)."""
    return hp / h


def h_dimensionless(h, rw, kv_kh=1.0):
    """hD = (h/rw) * sqrt(kh/kv).  GK-22: 120 isotropic, 169.7 at kv/kh = 0.5."""
    return (h / rw) * math.sqrt(1.0 / kv_kh)


# ------------------------------------------------------------ Brons–Marting
def g_brons_marting(b, hD):
    """G(b) auxiliary function as printed on the topic page:
    G(b) = 1.35 / [ b^0.825 * (1 + b^3.5/hD) ].
    Reproduces the worked example G(0.5) ≈ 2.40 at hD = 120.
    NOTE: this closed form degrades below b ≈ 0.4 (Sc becomes non-monotonic);
    the page's design table at low b matches the hD>20 approximation instead.
    """
    return 1.35 / (b ** 0.825 * (1.0 + b ** 3.5 / hD))


def brons_marting(b, hD):
    """Partial-completion skin, Brons & Marting (1961), as on the topic page:
    Sc = (1/b − 1) [ ln(hD/b) − G(b) ].
    Worked example: b=0.5, hD=120 → Sc ≈ +3.08 (page table prints +5.0 — quirk)."""
    if b >= 1.0:
        return 0.0
    return (1.0 / b - 1.0) * (math.log(hD / b) - g_brons_marting(b, hD))


def brons_marting_approx(b, h, rw):
    """hD > 20 approximation: Sc ≈ ((1−b)/b) [ ln(h/(b·rw)) − 0.75 ].
    This is the form the page's published design table tracks (e.g. +4.7 ≈ +5.0
    at b = 0.5, +12.2 ≈ +11.5 at b = 0.3)."""
    if b >= 1.0:
        return 0.0
    return ((1.0 - b) / b) * (math.log(h / (b * rw)) - 0.75)


# ---------------------------------------------------------------- Cinco-Ley
def theta_prime(theta_deg, kv_kh=1.0):
    """Anisotropy-corrected deviation angle: θ′ = arctan[ sqrt(kv/kh) · tanθ ]."""
    return math.degrees(math.atan(math.sqrt(kv_kh) * math.tan(math.radians(theta_deg))))


def cinco_ley(theta_deg, hD, kv_kh=1.0):
    """Deviation pseudo-skin, Cinco-Ley et al. (1975), BP NWP Eq 2.30:
    Sc″ = −(θ′/41)^2.06 − (θ′/56)^1.865 · log10(hD/100),   valid θ ≤ 75°.
    Worked example: θ=45°, isotropic, hD=120 → Sc″ ≈ −1.26 (page table −2.1 — quirk)."""
    tp = theta_prime(theta_deg, kv_kh)
    if tp > 75.0:
        raise ValueError("Cinco-Ley valid to 75 deg; use Joshi (1988) beyond")
    return -((tp / 41.0) ** 2.06) - ((tp / 56.0) ** 1.865) * math.log10(hD / 100.0)


# ------------------------------------------------- combining & Jones–Watts
def combined_completion_skin(Sc, Sc2):
    """Combined completion skin (BP NWP Eq 2.31): Sc,total = Sc + Sc″."""
    return Sc + Sc2


def jones_watts(h, hp, Sd, Sc=0.0):
    """Jones–Watts (Eq 4b.1): S = (h/hp)·Sd + Sc.
    Partial completion AMPLIFIES true damage — flow squeezes through hp of
    damaged rock, so each unit of Sd costs h/hp times more."""
    return (h / hp) * Sd + Sc


def j_index(k, h, mu, B, ln_term, S):
    """Productivity index J = 0.00708 k h / [ μ B (ln_term + S) ]  (stb/d/psi)."""
    return 0.00708 * k * h / (mu * B * (ln_term + S))


def j_ratio_ln7(S1, S2):
    """J2/J1 = (7 + S1)/(7 + S2) — classroom ln≈7 comparison form."""
    return (LN7 + S1) / (LN7 + S2)


# ------------------------------------------------ Joshi horizontal (DB-17)
def joshi_horizontal(kh, kv, h, L, rw, reh, mu, B):
    """Joshi (1988) horizontal-well PI with Economides anisotropy correction.
    Returns (Jh, X, Seq, Jv_ideal, ratio) where Seq is the equivalent skin a
    vertical well would need to match Jh:  X = ln(0.472 reh/rw) + Seq.
    Used ONLY for the SYNTHETIC Dunbar DB-17 training vignette."""
    beta = math.sqrt(kh / kv)
    a = (L / 2.0) * math.sqrt(0.5 + math.sqrt(0.25 + (2.0 * reh / L) ** 4))
    term1 = math.log((a + math.sqrt(a ** 2 - (L / 2.0) ** 2)) / (L / 2.0))
    term2 = (beta * h / L) * math.log(beta * h / (2.0 * rw))
    X = term1 + term2
    ln_v = math.log(0.472 * reh / rw)
    Jh = 0.00708 * kh * h / (mu * B * X)
    Jv = 0.00708 * kh * h / (mu * B * ln_v)
    return Jh, X, X - ln_v, Jv, Jh / Jv


# ---------------------------------------------------------------- __main__
if __name__ == '__main__':
    g = GK22
    print('=' * 74)
    print('Topic 3.4 — pseudo-skin reference calculations (TECH-BRIEF-m03 canon)')
    print('=' * 74)

    # 1. Brons–Marting worked example — 50%-perforated GK-22 hypothetical
    hD_iso = h_dimensionless(g['h'], g['rw'])                      # 120.0
    Sc_05 = brons_marting(0.5, hD_iso)
    print('\n[1] Brons–Marting, b = 0.5, hD = %.0f (GK-22 isotropic):' % hD_iso)
    print('    G(0.5)        = %.3f   (page: 2.401)' % g_brons_marting(0.5, hD_iso))
    print('    Sc            = +%.2f  (page worked example: +3.08; design table'
          ' prints +5.0 — see README)' % Sc_05)
    print('    approximation = +%.2f  (matches the published table, not the'
          ' worked chain)' % brons_marting_approx(0.5, g['h'], g['rw']))

    # 2. Cinco-Ley worked example — 45 deg deviated GK-22, isotropic
    Sc2_45 = cinco_ley(45.0, hD_iso)
    print('\n[2] Cinco-Ley, theta = 45 deg, isotropic, hD = 120:')
    print('    Sc"           = %.2f   (page worked example: -1.26; table prints'
          ' -2.1 — see README)' % Sc2_45)

    # 3. GK-22 actual audit — scenario (a)
    hD_aniso = h_dimensionless(g['h'], g['rw'], kv_kh=0.5)         # 169.7
    Sc_a = brons_marting(1.0, hD_aniso)
    Sc2_a = cinco_ley(4.0, hD_aniso, kv_kh=0.5)
    print('\n[3] GK-22 actual audit (b = 1.0, theta = 4 deg, kv/kh = 0.5,'
          ' hD = %.1f):' % hD_aniso)
    print('    theta-prime   = %.2f deg' % theta_prime(4.0, 0.5))
    print('    Sc            = %.3f   (exactly zero: full penetration)' % Sc_a)
    print('    Sc"           = %.3f  (canon: -0.005)' % Sc2_a)
    print('    with Dq = +0.001  ->  Sd = S\' = +14.0 CONFIRMED (all true damage)')

    # 4. Scenario (b): 45-deg deviated + 50% perforated (brief worked values)
    comb_b = combined_completion_skin(Sc_05, Sc2_45)
    print('\n[4] Scenario (b): b = 0.5 AND theta = 45 deg (isotropic, hD = 120):')
    print('    Sc,total      = +%.2f + (%.2f) = +%.2f' % (Sc_05, Sc2_45, comb_b))

    # 5. Jones–Watts worked row — 50%-perforated GK-22 with Sd = +14
    S_jw = jones_watts(g['h'], 21.0, 14.0, Sc_05)
    J_pub = j_index(g['k'], 21.0, g['mu'], g['B'], LN_TERM_GK22, 14.0 + Sc_05)
    J_strict = j_index(g['k'], g['h'], g['mu'], g['B'], LN_TERM_GK22, S_jw)
    dd = g['pbar'] - g['pwf']
    print('\n[5] Jones–Watts, GK-22 hypothetical (h 42, hp 21, Sd +14, Sc +3.08):')
    print('    S_total       = (42/21)x14 + 3.08 = +%.1f  (page: +31.1)' % S_jw)
    print('    J (published basis, k*hp) = %.3f stb/d/psi -> q = %.0f stb/d'
          % (J_pub, J_pub * dd))
    print('      (page prints J = 0.214, q = 364, extra loss 418 stb/d)')
    print('    J (strict k*h basis)      = %.3f stb/d/psi -> q = %.0f stb/d'
          % (J_strict, J_strict * dd))
    print('      (the two formulations bracket the answer — see README quirk note)')

    # 6. Acid-vs-perforation trade (b = 0.5, Sd = +10, ln≈7 basis)
    print('\n[6] Acid-vs-perfs trade (b = 0.5, Sd = +10, ln~7 basis):')
    S_now = jones_watts(2.0, 1.0, 10.0, 3.1)        # h/hp = 2 -> +23.1
    for label, S_new in [('do nothing', S_now), ('acid only (Sd->0)', 3.1),
                         ('perforate only (b->1)', 10.0), ('both', 0.0)]:
        frac = LN7 / (LN7 + S_new)
        print('    %-22s S\' = %+5.1f   J/Jideal = %.1f%%   gain = %.2fx'
              % (label, S_new, 100 * frac, j_ratio_ln7(S_now, S_new)))

    # 7. SYNTHETIC Dunbar DB-17 horizontal vignette (Joshi 1988)
    print('\n[7] Dunbar DB-17 horizontal vignette — SYNTHETIC, TRAINING USE ONLY:')
    Jh, X, Seq, Jv, ratio = joshi_horizontal(kh=45.0, kv=13.5, h=80.0, L=1500.0,
                                             rw=0.35, reh=2980.0, mu=0.9, B=1.25)
    print('    kh 45 md, kv/kh 0.3, h 80 ft, L 1,500 ft, reh 2,980 ft')
    print('    Joshi denominator X = %.3f   vs vertical ln term %.3f'
          % (X, math.log(0.472 * 2980.0 / 0.35)))
    print('    equivalent skin Seq = %.2f   (geometry acting as stimulation)' % Seq)
    print('    Jh = %.2f stb/d/psi vs Jv(S=0) = %.2f  ->  %.2fx vertical'
          % (Jh, Jv, ratio))
    print('    at 800 psi drawdown: qh = %.0f vs qv = %.0f stb/d'
          % (Jh * 800.0, Jv * 800.0))

    print('\nAll values match TECH-BRIEF-m03.md Topic 3.4 canon (worked-example'
          ' basis).')