"""FIELD3D ``.flp`` parameter file reader (port of Matlab ``read_flp3d.m``)."""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from at_py.readwrite.flp import _mlimit_from_line, _opt_from_line, _title_from_line
from at_py.readwrite.vector import format_readvector_lines, parse_readvector_lines
_MAX_MODES_FIELD3D = 4000
def _normalize_flp3d_opt(raw: str) -> str:
"""Default third option char ``O`` if ``len(Opt) <= 2`` (no 4th ``C``; ``read_flp3d.m``)."""
if len(raw) <= 2:
chars = list(raw)
while len(chars) < 3:
chars.append(" ")
chars[2] = "O"
return "".join(chars[:3])
return raw
def _fscanf_floats_slash_stops(
lines: list[str],
start: int,
n_want: int,
) -> tuple[list[float], int]:
"""Emulate ``fscanf('%f')``: ``/`` stops conversion (``readRcvrBearings.m``)."""
out: list[float] = []
j = start
while j < len(lines) and len(out) < n_want:
part = lines[j].split("!", 1)[0]
for tok in part.replace(",", " ").split():
tok = tok.strip()
if not tok:
continue
if tok == "/":
return out, j + 1
try:
out.append(float(tok))
except ValueError as e:
raise ValueError(f"invalid float token {tok!r} in receiver bearings") from e
if len(out) >= n_want:
return out[:n_want], j + 1
j += 1
return out, j
def _parse_rcvr_bearings(lines: list[str], i: int) -> tuple[np.ndarray, int]:
"""Port of ``readRcvrBearings.m`` (after count line at ``lines[i]``)."""
if i >= len(lines):
raise ValueError("receiver bearings: missing count line")
n_theta_req = _mlimit_from_line(lines[i])
i += 1
if i >= len(lines):
raise ValueError("receiver bearings: missing data")
floats, j = _fscanf_floats_slash_stops(lines, i, n_theta_req)
if n_theta_req > 2:
if len(floats) < 2:
raise ValueError("receiver bearings: need at least two values when Ntheta > 2")
theta = np.linspace(floats[0], floats[1], n_theta_req, dtype=np.float64)
else:
if len(floats) < n_theta_req:
raise ValueError(
f"receiver bearings: expected {n_theta_req} value(s), got {len(floats)}"
)
theta = np.array(floats[:n_theta_req], dtype=np.float64)
if len(theta) >= 2 and np.isclose(theta[-1], theta[0] + 360.0):
theta = theta[:-1]
return theta, j
def _parse_node_line(line: str) -> tuple[float, float, str]:
"""Parse one mesh node: ``x_km``, ``y_km``, quoted mode file name."""
s = line.split("!", 1)[0]
parts = s.split("'")
if len(parts) < 3:
raise ValueError(f"field3d .flp node line missing quoted mode file name: {line!r}")
name = parts[1].strip()
head = parts[0].split()
if len(head) < 2:
raise ValueError(f"field3d .flp node line needs x y before quoted name: {line!r}")
return float(head[0]), float(head[1]), name
def _parse_elt_line(line: str) -> tuple[int, int, int]:
"""Parse one triangle: three 1-based node indices."""
s = line.split("!", 1)[0]
toks = s.replace(",", " ").split()
if len(toks) < 3:
raise ValueError(f"field3d .flp element line needs three node indices: {line!r}")
return int(toks[0]), int(toks[1]), int(toks[2])
[docs]
@dataclass(frozen=True)
class FieldFlp3DResult:
"""In-memory ``field3d`` driver parameters from a ``.flp`` file (``read_flp3d.m``)."""
title: str
opt: str
component: str
m_limit: int
source_x_m: np.ndarray
source_y_m: np.ndarray
n_sx: int
n_sy: int
source_z_m: np.ndarray
receiver_z_m: np.ndarray
receiver_range_m: np.ndarray
r_min_m: float
r_max_m: float
theta_deg: np.ndarray
node_x_m: np.ndarray
node_y_m: np.ndarray
mode_file_names: tuple[str, ...]
elt_nodes: np.ndarray
"""Shape ``(3, n_elts)``, 1-based node indices as in the file."""
[docs]
def parse_field_flp3d(text: str) -> FieldFlp3DResult:
"""Parse FIELD3D ``.flp`` text (Matlab ``read_flp3d.m``).
Horizontal coordinates and receiver ranges in the file are **km**; results use **meters**.
Trailing lines after the triangular mesh block are ignored (extended FIELD3D inputs).
"""
lines = text.splitlines()
i = 0
if i >= len(lines):
raise ValueError("empty .flp file")
title = _title_from_line(lines[i])
i += 1
if i >= len(lines):
raise ValueError(".flp: missing Opt line")
opt_raw = _opt_from_line(lines[i])
i += 1
opt = _normalize_flp3d_opt(opt_raw)
comp = opt[2] if len(opt) >= 3 else "P"
if i >= len(lines):
raise ValueError(".flp: missing MLimit line")
m_limit = min(_MAX_MODES_FIELD3D, _mlimit_from_line(lines[i]))
i += 1
sx_km, n_sx, i = parse_readvector_lines(lines, i)
sy_km, n_sy, i = parse_readvector_lines(lines, i)
source_x_m = (sx_km * 1000.0).astype(np.float64)
source_y_m = (sy_km * 1000.0).astype(np.float64)
sz, _nsz, i = parse_readvector_lines(lines, i)
rz, _nrz, i = parse_readvector_lines(lines, i)
rr_km, _nrr, i = parse_readvector_lines(lines, i)
receiver_range_m = (rr_km * 1000.0).astype(np.float64)
r_min_m = float(receiver_range_m[0])
r_max_m = float(receiver_range_m[-1])
theta_deg, i = _parse_rcvr_bearings(lines, i)
if i >= len(lines):
raise ValueError("field3d .flp: missing NNodes line")
n_nodes = _mlimit_from_line(lines[i])
i += 1
xs_km: list[float] = []
ys_km: list[float] = []
names: list[str] = []
for _ in range(n_nodes):
if i >= len(lines):
raise ValueError("field3d .flp: incomplete node block")
xk, yk, nm = _parse_node_line(lines[i])
xs_km.append(xk)
ys_km.append(yk)
names.append(nm)
i += 1
if i >= len(lines):
raise ValueError("field3d .flp: missing NElts line")
n_elts = _mlimit_from_line(lines[i])
i += 1
elt = np.zeros((3, n_elts), dtype=np.int32)
for e in range(n_elts):
if i >= len(lines):
raise ValueError("field3d .flp: incomplete element block")
a, b, c = _parse_elt_line(lines[i])
elt[0, e] = a
elt[1, e] = b
elt[2, e] = c
i += 1
return FieldFlp3DResult(
title=title,
opt=opt,
component=comp,
m_limit=m_limit,
source_x_m=source_x_m,
source_y_m=source_y_m,
n_sx=n_sx,
n_sy=n_sy,
source_z_m=sz.astype(np.float64),
receiver_z_m=rz.astype(np.float64),
receiver_range_m=receiver_range_m,
r_min_m=r_min_m,
r_max_m=r_max_m,
theta_deg=theta_deg,
node_x_m=(np.array(xs_km, dtype=np.float64) * 1000.0),
node_y_m=(np.array(ys_km, dtype=np.float64) * 1000.0),
mode_file_names=tuple(names),
elt_nodes=elt,
)
[docs]
def parse_field_flp3d_bytes(data: bytes, *, encoding: str = "utf-8") -> FieldFlp3DResult:
"""Like :func:`parse_field_flp3d` after decoding ``data`` with ``encoding``."""
return parse_field_flp3d(data.decode(encoding, errors="replace"))
def format_rcvr_bearings_block(theta_deg: np.ndarray) -> str:
"""Emit receiver bearings block (inverse of :func:`_parse_rcvr_bearings`)."""
t = np.asarray(theta_deg, dtype=np.float64).reshape(-1)
n = int(t.size)
if n < 1:
raise ValueError("receiver bearings: empty array")
if n > 2:
lin = np.linspace(float(t[0]), float(t[-1]), n, dtype=np.float64)
if np.allclose(t, lin):
return f"{n}\n{t[0]:.17g} {t[-1]:.17g} /\n"
return f"{n}\n" + " ".join(f"{float(x):.17g}" for x in t) + "\n"
return f"{n}\n" + " ".join(f"{float(x):.17g}" for x in t) + "\n"