#!/usr/bin/env python3
"""
hawkins.py — Hawkins' formation-damage skin (Sd) reference implementation.

Well Productivity Programme · Course 01 · Module 03 · Topic 3.3 (Formation Damage Skin).
Running case: Gashaka GK-22 skin audit (Niger Delta, Agbada Fm).

All numbers reproduce TECH-BRIEF-m03 §1a and the Lecture 3.3a / 3.3b / 3.3c decks:

    single-zone inverse : ks/k = 0.145, Sd = +14   ->  rs = 3.77 ft
    two-zone composite  : solids 6.80 mD -> 0.60 ft  =>  +6.19
                          filtrate 21.25 mD -> 8.08 ft =>  +7.81   (sum ~= +14)
    post-acid (single)  : ks,post = 0.76 k           ->  Sd ~= +1.0

Hawkins (1956):  Sd = (k/ks - 1) * ln(rs/rw)
Multi-zone (series, outer ratios use formation k — see TECH-BRIEF inconsistency note):
    Sd_total = sum_i (k/ks_i - 1) * ln(r_outer_i / r_inner_i)

Run:  python3 hawkins.py
"""

import math

# --- GK-22 canonical inputs (TECH-BRIEF-m03 §1a, LOCKED) ---
K = 85.0      # formation permeability, mD
RW = 0.35     # wellbore radius, ft
PHI = 0.22    # porosity, fraction
H = 42.0      # net pay, ft
RE = 1650.0   # drainage radius, ft


# ---------------------------------------------------------------------------
# Core Hawkins relations
# ---------------------------------------------------------------------------
def sd_forward(ks_over_k, rs, rw=RW):
    """Forward Hawkins: anatomy (severity ks/k, depth rs) -> skin Sd."""
    return (1.0 / ks_over_k - 1.0) * math.log(rs / rw)


def rs_inverse(sd, ks_over_k, rw=RW):
    """Inverse for depth: skin Sd + severity ks/k -> damaged radius rs (ft)."""
    bracket = 1.0 / ks_over_k - 1.0          # = k/ks - 1
    return rw * math.exp(sd / bracket)


def ks_over_k_inverse(sd, rs, rw=RW):
    """Inverse for severity: skin Sd + depth rs -> ks/k (the ambiguity curve)."""
    return 1.0 / (1.0 + sd / math.log(rs / rw))


def sd_composite(zones, rw=RW, k=K):
    """
    Multi-zone composite skin. `zones` is an ordered list (inner -> outer) of
    dicts: {'ks': ks_mD, 'r_outer': ft}. Each zone's ratio uses formation k;
    the inner radius of zone i is the outer radius of zone i-1 (rw for the first).
    Returns (total_Sd, [per-zone Sd]).
    """
    terms = []
    r_inner = rw
    for z in zones:
        ks = z['ks']
        r_outer = z['r_outer']
        terms.append((k / ks - 1.0) * math.log(r_outer / r_inner))
        r_inner = r_outer
    return sum(terms), terms


def solve_zone2_radius(sd_zone2, ks2, r_inner, k=K):
    """Solve the outer radius of zone 2 so its Hawkins term equals sd_zone2."""
    bracket = k / ks2 - 1.0
    return r_inner * math.exp(sd_zone2 / bracket)


# ---------------------------------------------------------------------------
# Derived quantities
# ---------------------------------------------------------------------------
def damaged_volumes(rs, h=H, phi=PHI, rw=RW):
    """Rock volume and pore volume of the damaged annulus rw->rs."""
    rock_ft3 = math.pi * (rs ** 2 - rw ** 2) * h
    rock_bbl = rock_ft3 / 5.615
    pore_bbl = rock_bbl * phi
    return rock_ft3, rock_bbl, pore_bbl


def j_ratio(sd_before, sd_after, ln_term=7.0):
    """Productivity ratio J_after/J_before using the ln~=7 approximation (Topic 3.1)."""
    return (ln_term + sd_before) / (ln_term + sd_after)


def acid_scenario(restoration, sd_pre, ks_over_k_pre, rw=RW):
    """
    Single-zone acid response. `restoration` = fraction of formation k restored
    (e.g. 0.76). Keeps the damaged radius rs fixed at the pre-acid value and
    recomputes Sd at the new, milder severity.
    """
    rs = rs_inverse(sd_pre, ks_over_k_pre, rw)
    ks_over_k_post = restoration           # ks,post / k
    sd_post = sd_forward(ks_over_k_post, rs, rw)
    return sd_post, rs


# ---------------------------------------------------------------------------
# Demonstration / self-check
# ---------------------------------------------------------------------------
def _check(label, value, expected, tol):
    ok = abs(value - expected) <= tol
    flag = 'OK ' if ok else 'XX '
    print(f"  [{flag}] {label:<42} {value:8.3f}  (expect ~{expected})")
    return ok


def main():
    print("=" * 68)
    print("GK-22 HAWKINS FORMATION-DAMAGE AUDIT  (Topic 3.3)")
    print(f"  k = {K} mD · rw = {RW} ft · h = {H} ft · phi = {PHI} · re = {RE} ft")
    print("=" * 68)
    all_ok = True

    # 1 · Single-zone inverse: ks/k = 0.145, Sd = +14 -> rs = 3.77 ft
    print("\n1 · Single-zone inverse (measured severity ks/k = 0.145):")
    ks_over_k = 0.145
    rs = rs_inverse(14.0, ks_over_k)
    all_ok &= _check("rs from Sd=+14, ks/k=0.145 (ft)", rs, 3.77, 0.03)
    sd_chk = sd_forward(ks_over_k, rs)
    all_ok &= _check("forward check Sd", sd_chk, 14.0, 0.05)
    rock_ft3, rock_bbl, pore_bbl = damaged_volumes(rs)
    all_ok &= _check("damaged rock volume (ft3)", rock_ft3, 1860.0, 40.0)
    all_ok &= _check("damaged rock (bbl)", rock_bbl, 331.0, 8.0)
    all_ok &= _check("damaged pore volume (bbl)", pore_bbl, 73.0, 3.0)

    # 2 · Two-zone composite: +6.19 (solids) + 7.81 (filtrate) ~= +14
    print("\n2 · Two-zone composite (solids collar + filtrate halo):")
    # Zone 1: mud solids, ks1 = 0.08k = 6.80 mD, rs1 = 0.60 ft
    ks1 = 0.08 * K
    sd1 = sd_forward(ks1 / K, 0.60)
    all_ok &= _check("Zone 1 (6.80 mD -> 0.60 ft) Sd", sd1, 6.19, 0.05)
    # Zone 2: filtrate/clay, ks2 = 0.25k = 21.25 mD; size rs2 so term = 7.81
    ks2 = 0.25 * K
    rs2 = solve_zone2_radius(7.81, ks2, r_inner=0.60)
    all_ok &= _check("Zone 2 outer radius rs2 (ft)", rs2, 8.08, 0.05)
    total, terms = sd_composite(
        [{'ks': ks1, 'r_outer': 0.60}, {'ks': ks2, 'r_outer': rs2}]
    )
    all_ok &= _check("Zone 2 (21.25 mD -> 8.08 ft) Sd", terms[1], 7.81, 0.05)
    all_ok &= _check("composite total Sd", total, 14.0, 0.05)

    # 3 · Inverse ambiguity sanity: shallow&savage vs deep&mild both give +14
    print("\n3 · Inverse-ambiguity curve at fixed Sd = +14:")
    all_ok &= _check("ks/k for rs = 1.0 ft", ks_over_k_inverse(14.0, 1.0), 0.0697, 0.01)
    all_ok &= _check("ks/k for rs = 12.0 ft", ks_over_k_inverse(14.0, 12.0), 0.202, 0.01)

    # 4 · Post-acid single-zone: 76% restoration -> Sd ~= +1.0
    print("\n4 · Acid response (76% average permeability restoration):")
    sd_post, rs_post = acid_scenario(0.76, sd_pre=14.0, ks_over_k_pre=0.145)
    # Brief: Sd,post = (1/0.76 - 1)*ln(3.77/0.35) = 0.316*2.377 = 0.75, quoted ~= +1.0
    all_ok &= _check("post-acid Sd (0.75, rounds to +1.0)", sd_post, 0.75, 0.05)
    # Canon J ratio / q use the rounded post-acid Sd = +1 (TECH-BRIEF QUICK CANON CARD)
    jr = j_ratio(14.0, 1.0)
    all_ok &= _check("J ratio at Sd +14 -> +1 (21/8)", jr, 2.625, 0.01)
    q_post = 782.0 * jr
    all_ok &= _check("q post-acid at pwf 2,500 (stb/d)", q_post, 2053.0, 6.0)

    print("\n" + "=" * 68)
    print("ALL CHECKS PASSED" if all_ok else "*** SOME CHECKS FAILED ***")
    print("=" * 68)
    return 0 if all_ok else 1


if __name__ == "__main__":
    raise SystemExit(main())