#!/usr/bin/env python3
"""Reproducible illustrative petrophysical LAS; no measured well data.

Requires numpy and scipy. Run in any directory; outputs beside this script.
The model intentionally separates facies, fluids, invasion and tool smoothing.
It is an educational response model, not a calibrated tool forward model.
"""
from pathlib import Path
import json
import numpy as np
from scipy.ndimage import gaussian_filter1d

OUT = Path(__file__).resolve().parent
SEED = 20260929
rng = np.random.default_rng(SEED)
z = np.arange(500.0, 1500.0001, 0.25)
n = z.size
# top, base, facies, description, clean-fraction porosity, shale volume, Sw, calcite in matrix
zones = [
    (500, 620, 1, "Shale", .10, .87, 1.00, .00),
    (620, 800, 2, "Gas sandstone", .265, .025, .18, .00),
    (800, 820, 1, "Thin shale", .10, .91, 1.00, .00),
    (820, 960, 3, "Oil sandstone upper", .245, .035, .24, .00),
    (960, 968, 4, "Tight calcite-cemented sandstone", .028, .012, 1.00, .45),
    (968, 1110, 3, "Oil sandstone lower", .225, .045, .28, .00),
    (1110, 1210, 1, "Shale", .095, .86, 1.00, .00),
    (1210, 1410, 5, "Brine sandstone", .255, .025, 1.00, .00),
    (1410, 1500, 1, "Basal shale", .085, .89, 1.00, .00),
]

def band(sigma):
    x = gaussian_filter1d(rng.normal(size=n), sigma=sigma, mode="nearest")
    return (x-x.mean())/x.std()

broad, bed, fine = band(30), band(5), band(1)
facies = np.zeros(n, int)
zone_number = np.zeros(n, int)
phi = np.zeros(n)
vsh = np.zeros(n)
sw = np.ones(n)
calc = np.zeros(n)
for k, (top, base, fac, name, p, sh, sat, cc) in enumerate(zones):
    mask = (z >= top) & ((z < base) | ((k == len(zones)-1) & (z <= base)))
    frac = (z[mask]-top)/(base-top)
    facies[mask] = fac
    zone_number[mask] = k+1
    amp = .0025 if fac == 4 else (.006 if fac == 1 else .017)
    trend = .008*np.sin(2*np.pi*frac) if fac in (2,3,5) else 0
    phi[mask] = p+amp*(.7*broad[mask]+.35*bed[mask])+trend
    vsh[mask] = np.clip(sh+(.038 if fac == 1 else .010)*bed[mask], .004, .98)
    sw[mask] = np.clip(sat+.023*broad[mask]+.009*bed[mask], .12, .45) if fac in (2,3) else 1.0
    calc[mask] = cc
phi = np.clip(phi,.015,.34)
gas = facies == 2
oil = facies == 3
tight = facies == 4
shale = facies == 1

# Nuclear tools see partial invasion. Resistivity deep channel sees the virgin zone.
sx = np.ones(n)
sx[gas] = np.clip(sw[gas]+.26, .38, .62)
sx[oil] = np.clip(sw[oil]+.20, .35, .65)
rho_hc = np.ones(n)*1.025
rho_hc[gas] = .12
rho_hc[oil] = .80
hi_hc = np.ones(n)
hi_hc[gas] = .20
hi_hc[oil] = .95
rho_fl = sx*1.025+(1-sx)*rho_hc
hi_fl = sx+(1-sx)*hi_hc
rho_ma = (1-calc)*2.65+calc*2.71
rho_sh = 2.46+.000065*(z-500)+.010*broad
nphi_sh = .39-.000025*(z-500)+.011*bed
dt_sh = 115-.015*(z-500)+3.0*broad

gr_raw = 18+120*vsh+2.8*bed+1.3*fine
rhob_raw = (1-vsh)*((1-phi)*rho_ma+phi*rho_fl)+vsh*rho_sh
rhob_raw += .0035*fine
nphi_raw = (1-vsh)*(phi*hi_fl+.012*calc)+vsh*nphi_sh+.0025*fine
# Wyllie-style liquid reference; gas delay is an explicit illustrative 8-13 us/ft.
# It is not Wyllie applied to a free-gas velocity and is not Gassmann inversion.
dt_ma = (1-calc)*55.5+calc*47.5
dt_sand = dt_ma+phi*(189-dt_ma)
dt_sand += np.where(gas, 10.5+1.3*broad, np.where(oil, 1.5, 0))
dt_raw = (1-vsh)*dt_sand+vsh*dt_sh+.65*fine

# Archie response in clean sands, plus a small parallel shale conductivity term.
# Rw is fixed at formation conditions to isolate the fluid/facies effects.
rw = .12
m = np.where(tight, 2.15, 2.0)
r_clean = rw/(phi**m*sw**2)
r_shale = 1.7*np.exp(.14*broad)
rt = 1/((1-vsh)/r_clean+.018*vsh/r_shale)
rt[shale] = r_shale[shale]
rs_clean = rw/(phi**m*sx**2)
rs = 1/((1-vsh)/rs_clean+.018*vsh/r_shale)
rs[shale] = rt[shale]*np.exp(.025*bed[shale])
rt *= np.exp(.025*fine)
rs *= np.exp(.018*fine)

def smooth(x, metres):
    return gaussian_filter1d(x, metres/.25, mode="nearest")

gr = smooth(gr_raw,.35)
rhob = smooth(rhob_raw,.25)
nphi = smooth(nphi_raw,.45)
dt = smooth(dt_raw,.35)
lld = np.exp(smooth(np.log(rt),.70))
lls = np.exp(smooth(np.log(rs),.45))
cali = smooth(8.5+np.where(shale,.7+.28*np.maximum(bed,0),-.04)+.045*fine,.35)
drho = smooth(np.where(shale,.035,.008)+.004*bed,.25)
phie = (1-vsh)*phi
curves = [
    ("DEPT","M","Measured depth",z),
    ("GR","GAPI","Gamma ray",gr),
    ("CALI","IN","Caliper; bit size 8.5 in",cali),
    ("RHOB","G/CC","Bulk density",rhob),
    ("NPHI","V/V","Neutron porosity sandstone units",nphi),
    ("DT","US/FT","Compressional sonic slowness",dt),
    ("LLD","OHMM","Deep resistivity",lld),
    ("LLS","OHMM","Shallow resistivity partial invasion",lls),
    ("DRHO","G/CC","Density correction indicator",drho),
    ("PHIE","V/V","SYNTHETIC TRUTH effective porosity",phie),
    ("SW","V/V","SYNTHETIC TRUTH virgin water saturation",sw),
    ("VSH","V/V","SYNTHETIC TRUTH shale volume",vsh),
    ("FACIES","CODE","TRUTH 1 shale 2 gas 3 oil 4 tight calcite 5 brine",facies),
    ("ZONE","CODE","TRUTH sequential layer number 1 to 9",zone_number),
]
data = np.column_stack([c[3] for c in curves])
las = OUT/"Petrophysics_SYNTH_500-1500m.las"
header = [
    "~Version Information",
    "VERS. 2.0 : CWLS LAS version",
    "WRAP. NO : One depth sample per line",
    "~Well Information",
    "STRT.M 500.000 : Start depth",
    "STOP.M 1500.000 : Stop depth",
    "STEP.M 0.250 : Depth increment",
    "NULL. -999.25 : Missing value",
    "COMP. PETROPHYSICS.IO : Company",
    "WELL. SYNTH-500-1500 : Synthetic well name",
    "FLD. SYNTHETIC NINE LAYER MODEL : Field",
    "LOC. SYNTHETIC - NO REAL LOCATION : Location",
    "DATE. 2026-09-29 : Generation date",
    "UWI. SYNTHETIC-20260929 : Synthetic identifier",
    "~Parameter Information",
    "BS.IN 8.5 : Bit size",
    "RW.OHMM 0.12 : Assumed brine resistivity at formation conditions",
    "NREF. SANDSTONE : Neutron lithology calibration",
    "SEED. 20260929 : Reproducibility seed",
    "~Curve Information",
]
header += [f"{name}.{unit} : {desc}" for name,unit,desc,_ in curves]
header += ["~Other Information", "SYNTHETIC EDUCATIONAL DATA. No real well data.",
           "Shale means clay-rich shale/mudstone, not pure mineral clay.",
           "Log responses include correlated layering noise and 0.25-0.70 m Gaussian smoothing.",
           "Partial invasion in HC sands; a fixed formation-brine Rw is assumed.",
           "SW PHIE VSH FACIES ZONE are known model truth, not interpreted log results.",
           "Cemented interval is deliberately water-filled: high resistivity is caused by low porosity."]
header += [f"ZONE {i+1}: {a:.0f}-{b:.0f} M; {name}" for i,(a,b,_,name,*_) in enumerate(zones)]
header += ["Reference principles: SLB Defining Porosity; SLB Basic Well Log Interpretation;",
           "SLB sandstone-compatible scale; EPA Acoustic Logging; Bunch (2020) doi:10.1016/j.marpetgeo.2020.104424.",
           "~ASCII " + " ".join(c[0] for c in curves)]
with las.open("w",encoding="ascii",newline="\n") as f:
    f.write("\n".join(header)+"\n")
    np.savetxt(f,data,fmt=["%.3f"]+["%.5f"]*(data.shape[1]-1))

# Independent readback and geological plausibility checks.
text = las.read_text()
lines = text.splitlines()
ai = next(i for i,s in enumerate(lines) if s.upper().startswith("~ASCII"))
readback = np.loadtxt(lines[ai+1:])
assert readback.shape == (4001,14)
assert readback[0,0] == 500 and readback[-1,0] == 1500
assert np.allclose(np.diff(readback[:,0]),.25)
assert np.isfinite(readback).all() and np.all(readback[:,6:8] > 0)
assert np.allclose(readback,data,atol=.000051)
interior = np.ones(n,bool)
for a,b,*_ in zones:
    interior &= (abs(z-a)>2)&(abs(z-b)>2)
median = lambda x,fac: float(np.median(x[(facies==fac)&interior]))
dphi=(2.65-rhob)/(2.65-1.025)
assert median(gr,1)>100 and all(median(gr,f)<45 for f in [2,3,4,5])
assert median(dphi-nphi,2)>.12
assert abs(median(dphi-nphi,5))<.025
assert median(nphi-dphi,1)>.17
assert median(lld,2)>10*median(lld,5)
assert median(lld,3)>7*median(lld,5)
assert median(rhob,4)>2.58 and median(dt,4)<65 and median(nphi,4)<.065
stats = []
for i,(a,b,fac,name,*_) in enumerate(zones):
    ix = (zone_number==i+1)&interior
    stats.append(dict(top=a,base=b,facies=fac,name=name,**{nm:round(float(np.median(arr[ix])),4) for nm,_,_,arr in curves[1:9]}))
(OUT/"qc_summary.json").write_text(json.dumps(dict(samples=n,curves=len(curves),zones=stats),indent=2))
print(json.dumps(dict(file=str(las),samples=n,curves=len(curves),zones=stats),indent=2))