""" rns_helpers.py ============== Thin Python wrapper around the RNS v1.1d C executable. Provides: * run_rns(...) : run a single RNS computation, return a dict of results * run_sequence(...) : run a sequence (with -l/-n) and return a DataFrame * Several convenience helpers for the common sequence types * parse_metric_grid(...) : parse -p 3 (full metric and pressure on the grid) * load_eos(file) : read an EOS table for plotting * Physical constants matching those built into rns.c The wrapper assumes the executables `rns`, `rns_high`, `rns_vh` and an `eos/` subdirectory exist in BUILD_DIR (default: ./rns_build). Author: built for the RNS tutorial notebook. """ from __future__ import annotations import os import re import subprocess from dataclasses import dataclass from pathlib import Path from typing import Optional, Sequence import numpy as np import pandas as pd # ------------------------------------------------------------------ # Physical constants (matching the values hard-coded in rns.c) # ------------------------------------------------------------------ C_LIGHT = 2.9979e10 # cm/s G_NEWT = 6.6732e-8 # cgs M_SUN = 1.987e33 # g M_BARYON = 1.66e-24 # g KM_PER_CM = 1e-5 GM_SUN_C2_KM = G_NEWT * M_SUN / (C_LIGHT**2) * KM_PER_CM # ~ 1.477 km # ------------------------------------------------------------------ # Default build directory (override via environment or set_build_dir) # ------------------------------------------------------------------ BUILD_DIR = Path(os.environ.get("RNS_BUILD_DIR", "rns_build")).resolve() def set_build_dir(path: str | os.PathLike) -> None: """Point the wrapper at a different build directory.""" global BUILD_DIR BUILD_DIR = Path(path).resolve() def _binary(grid: str = "std") -> Path: name = { "std": "rns", "high": "rns_high", "vh": "rns_vh", "xh": "rns_xh", "hng": "hng", }[grid] path = BUILD_DIR / name if not path.exists(): raise FileNotFoundError( f"RNS binary not found at {path}. " f"Compile with `make` inside {BUILD_DIR}." ) return path def _eos_path(eos: str) -> str: """Resolve `eos='eosC'` to a path; pass-through if already a path.""" if os.path.sep in eos or os.path.exists(eos): return str(eos) return str(BUILD_DIR / "eos" / eos) # ------------------------------------------------------------------ # Output parsing (-p 1, vertical) # ------------------------------------------------------------------ # Maps the human-readable label printed by rns.c (case 1) to a short key. _P1_LABELS = { "e_c": "e_c", # central energy density / (10^15 g/cm^3) "M": "M", # gravitational mass (M_sun) "M_0": "M_0", # rest mass (M_sun) "R_e": "R_e", # equatorial circumferential radius (km) "Omega": "Omega", # angular velocity (10^4 s^-1) "Omega_p": "Omega_p", # Kepler angular velocity of test particle (10^4 s^-1) "T/W": "T_W", # rotational/gravitational energy ratio "cJ/GM_sun^2": "cJ_GMsun2", # dimensionless angular momentum "I": "I", # moment of inertia (10^45 g cm^2) "Phi_2": "Phi_2", # quadrupole moment (10^42 g cm^2) "h+": "h_plus", # ISCO height, co-rotating (km) "h-": "h_minus", # ISCO height, counter-rotating (km) "Z_p": "Z_p", # polar redshift "Z_f": "Z_f", # forward equatorial redshift "Z_b": "Z_b", # backward equatorial redshift "omega_c/Omega": "omegac_Omega", # central frame-drag / Omega "r_e": "r_e_km", # coordinate equatorial radius (km, tab EOS only) "r_p/r_e": "r_ratio", # axes ratio } _P1_RE = re.compile( r"^\s*([-+]?(?:[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?|nan|inf)|---)" r"\s+(\S+(?:\s*\([^)]*\))?)\s*$" ) def _to_float(s: str) -> float: """Convert a printed numeric token (including '---', '-nan', 'nan') to float.""" if s == "---": return np.nan try: return float(s) # handles 'nan' and '-nan' natively except ValueError: return np.nan def _parse_p1(text: str) -> dict[str, float]: """Parse vertical (-p 1) output. Returns a dict (NaN for '---' or 'nan').""" out: dict[str, float] = {} for line in text.splitlines(): m = _P1_RE.match(line) if not m: continue val_str, label_with_unit = m.groups() # strip units in parens to get the bare label label = re.sub(r"\s*\([^)]*\)\s*$", "", label_with_unit).strip() if label not in _P1_LABELS: continue key = _P1_LABELS[label] out[key] = _to_float(val_str) return out def parse_p1_rows(text: str) -> list[tuple[str, float, str]]: """Parse vertical (-p 1) output into (quantity, value, unit) rows.""" rows: list[tuple[str, float, str]] = [] for line in text.splitlines(): m = _P1_RE.match(line) if not m: continue val_str, label_with_unit = m.groups() quantity = label_with_unit.split("(", 1)[0].strip() unit = "" if "(" in label_with_unit and ")" in label_with_unit: unit = label_with_unit[label_with_unit.find("(") + 1 : label_with_unit.rfind(")")].strip() rows.append((quantity, _to_float(val_str), unit)) return rows # ------------------------------------------------------------------ # Output parsing (-p 2, horizontal table for sequences) # ------------------------------------------------------------------ # Column order from rns.c print_header() (case 2): # rho_c M M_0 R Omega Omega_p T/W cJ/GM_s^2 I Phi_2 h_plus h_minus # Z_p Z_b Z_f omega_c/Omega r_e r_ratio Omega_pa Omega+ u_phi # (Phi_2 only printed when r_ratio != 1) _P2_COLS_ROT = [ "rho_c", "M", "M_0", "R", "Omega", "Omega_p", "T_W", "cJ_GMsun2", "I", "Phi_2", "h_plus", "h_minus", "Z_p", "Z_b", "Z_f", "omegac_Omega", "r_e_km", "r_ratio", "Omega_pa", "Omega_plus", "u_phi", ] _P2_COLS_STAT = [ "rho_c", "M", "M_0", "R", "Omega", "Omega_p", "T_W", "cJ_GMsun2", "I", "h_plus", "h_minus", "Z_p", "Z_b", "Z_f", "omegac_Omega", "r_e_km", "r_ratio", "Omega_pa", "Omega_plus", "u_phi", ] def _parse_p2(text: str) -> pd.DataFrame: """Parse horizontal (-p 2) sequence output into a DataFrame.""" rows = [] for line in text.splitlines(): s = line.strip() if not s: continue if s.startswith("-") or s.startswith("rho_c") or s.startswith("eos") \ or s.startswith("N=") or "MDIVxSDIV" in s: continue # Replace '---' with NaN sentinels, split, convert tokens = s.split() try: values = [np.nan if t == "---" else float(t) for t in tokens] except ValueError: continue rows.append(values) if not rows: return pd.DataFrame() width = max(len(r) for r in rows) cols = _P2_COLS_ROT if width == len(_P2_COLS_ROT) else _P2_COLS_STAT # pad short rows rows = [r + [np.nan]*(len(cols)-len(r)) for r in rows] df = pd.DataFrame(rows, columns=cols) for col in ("h_plus", "h_minus"): if col in df.columns: df[col] = df[col].clip(lower=0) return df # ------------------------------------------------------------------ # Core driver # ------------------------------------------------------------------ @dataclass class RNSResult: """Container for a single RNS run.""" task: str params: dict raw: str data: dict # parsed scalar results (always present) table: Optional[pd.DataFrame] = None # populated for sequences (-n>1) or -p 2 def __getitem__(self, key): return self.data[key] def get(self, key, default=None): return self.data.get(key, default) _FLAG_MAP = { "eos": "f", "task": "t", "e": "e", "e_max": "l", "n": "n", "r": "r", "M": "m", "M0": "z", "Omega": "o", "J": "j", "p": "p", "d": "d", "accuracy": "a", "fix_error": "b", "cf": "c", "N": "N", # polytropic index } def run_rns( task: str = "static", eos: Optional[str] = None, *, polytrope_N: Optional[float] = None, e: Optional[float] = None, e_max: Optional[float] = None, n: int = 1, r: Optional[float] = None, M: Optional[float] = None, M0: Optional[float] = None, Omega: Optional[float] = None, J: Optional[float] = None, p_format: int = 1, accuracy: Optional[float] = None, fix_error: Optional[float] = None, cf: Optional[float] = None, grid: str = "std", timeout: float = 300.0, quiet: bool = True, ) -> RNSResult: """ Drive one RNS computation. Returns an RNSResult. Parameters ---------- task : {'static','model','gmass','rmass','omega','jmoment','kepler','test'} The RNS task to run. eos : str Tabulated EOS name (looked up in BUILD_DIR/eos) or path. Mutually exclusive with `polytrope_N`. polytrope_N : float Polytropic index for polytropic EOS (sets -q poly -N ). e : float Central energy density in g/cm^3 (tabulated) or dimensionless (polytropic). The flag value in rns.c. e_max, n : float, int For sequences: end-point central density (`-l`) and number of models (`-n`). With n > 1, results are parsed via -p 2. r, M, M0, Omega, J : float The varying parameters: axes ratio, gravitational mass (M_sun), rest mass (M_sun), angular velocity (10^4 s^-1), angular momentum (cJ/GM_sun^2). Only set the one(s) the task requires. p_format : 1 or 2 Vertical (1) or horizontal (2). For sequences (n>1), 2 is forced. grid : {'std','high','vh'} Which compiled grid resolution to use. """ if n > 1: p_format = 2 # Build command line cmd: list[str] = [str(_binary(grid))] if polytrope_N is not None: cmd += ["-q", "poly", "-N", f"{polytrope_N}"] else: if eos is None: raise ValueError("Provide either `eos` or `polytrope_N`.") cmd += ["-q", "tab", "-f", _eos_path(eos)] cmd += ["-t", task, "-p", str(p_format), "-d", "0"] if e is not None: cmd += ["-e", f"{e:.10e}"] if e_max is not None: cmd += ["-l", f"{e_max:.10e}"] if n > 1: cmd += ["-n", str(n)] if r is not None: cmd += ["-r", f"{r:.10f}"] if M is not None: cmd += ["-m", f"{M:.10f}"] if M0 is not None: cmd += ["-z", f"{M0:.10f}"] if Omega is not None: cmd += ["-o", f"{Omega:.10f}"] if J is not None: cmd += ["-j", f"{J:.10f}"] if accuracy is not None: cmd += ["-a", f"{accuracy:.3e}"] if fix_error is not None: cmd += ["-b", f"{fix_error:.3e}"] if cf is not None: cmd += ["-c", f"{cf:.3f}"] if not quiet: print(" ".join(cmd)) proc = subprocess.run( cmd, capture_output=True, text=True, timeout=timeout ) raw = proc.stdout # Parse table = None data: dict = {} if p_format == 2: table = _parse_p2(raw) # Pick the last row as the canonical 'data' for convenience if len(table): data = table.iloc[-1].to_dict() else: data = _parse_p1(raw) return RNSResult( task=task, params=dict(eos=eos, polytrope_N=polytrope_N, e=e, e_max=e_max, n=n, r=r, M=M, M0=M0, Omega=Omega, J=J, p_format=p_format, grid=grid), raw=raw, data=data, table=table, ) # ------------------------------------------------------------------ # Convenience wrappers for the common sequence types # ------------------------------------------------------------------ def static_sequence(eos: str, e_min: float, e_max: float, n: int = 30, grid: str = "std") -> pd.DataFrame: """TOV (nonrotating) sequence over a range of central energy densities.""" df = run_rns(task="static", eos=eos, e=e_min, e_max=e_max, n=n, grid=grid).table return df.dropna(subset=["M"]).reset_index(drop=True) def kepler_sequence(eos: str, e_min: float, e_max: float, n: int = 30, grid: str = "std", fix_error: float = 1e-4) -> pd.DataFrame: """Mass-shedding (Keplerian) sequence.""" df = run_rns(task="kepler", eos=eos, e=e_min, e_max=e_max, n=n, fix_error=fix_error, grid=grid).table return df.dropna(subset=["M"]).reset_index(drop=True) def constant_rest_mass_sequence( eos: str, M0: float, e_min: float, e_max: float, n: int = 20, grid: str = "std", fix_error: float = 1e-4, ) -> pd.DataFrame: """Constant rest-mass sequence (an isolated NS evolutionary track).""" df = run_rns(task="rmass", eos=eos, M0=M0, e=e_min, e_max=e_max, n=n, fix_error=fix_error, grid=grid).table return df.dropna(subset=["M"]).reset_index(drop=True) def constant_M0_track_safe( eos: str, M0: float, e_min: float, e_max: float, n: int = 12, grid: str = "std", fix_error: float = 1e-4, per_model_timeout: float = 8.0, include_endpoints: bool = False, endpoint_n: int = 80, ) -> pd.DataFrame: """ Build a constant-M_0 track by calling rmass for each e_c separately, with a strict per-model timeout. Points where RNS cannot converge (because M_0 does not exist at that density) silently drop out. Much more robust than constant_rest_mass_sequence(...) for tracks near the Kepler envelope, where the rmass bisection can fail to terminate. The returned DataFrame uses the same column names as the -p 2 sequence helpers (rho_c in g/cm^3, R in km, etc.). If include_endpoints=True, append one nonrotating (TOV) and one Kepler endpoint by selecting the closest-M_0 models from dense static/kepler sequences over [e_min, e_max]. """ def _p1_to_row(d: dict[str, float]) -> dict[str, float]: return { "rho_c": d["e_c"] * 1e15, # back to g/cm^3 "M": d["M"], "M_0": d["M_0"], "R": d["R_e"], "Omega": d["Omega"], "Omega_p": d["Omega_p"], "T_W": d["T_W"], "cJ_GMsun2": d["cJ_GMsun2"], "I": d.get("I", np.nan), "Phi_2": d.get("Phi_2", np.nan), "h_plus": max(0.0, d.get("h_plus", 0.0)) if not np.isnan(d.get("h_plus", np.nan)) else np.nan, "h_minus": max(0.0, d.get("h_minus", 0.0)) if not np.isnan(d.get("h_minus", np.nan)) else np.nan, "Z_p": d["Z_p"], "Z_b": d["Z_b"], "Z_f": d["Z_f"], "omegac_Omega": d.get("omegac_Omega", np.nan), "r_e_km": d.get("r_e_km", np.nan), "r_ratio": d["r_ratio"], } def _nearest_m0_row(df: pd.DataFrame, target_m0: float) -> Optional[dict[str, float]]: if df is None or len(df) == 0 or "M_0" not in df.columns: return None d = df.dropna(subset=["M_0"]) if len(d) == 0: return None idx = (d["M_0"] - target_m0).abs().idxmin() return d.loc[idx].to_dict() e_vals = np.geomspace(e_min, e_max, n) rows = [] for ec in e_vals: try: r = run_rns(task="rmass", eos=eos, M0=M0, e=ec, p_format=1, fix_error=fix_error, grid=grid, timeout=per_model_timeout) d = r.data if not np.isnan(d.get("M", np.nan)): rows.append(_p1_to_row(d)) except subprocess.TimeoutExpired: pass if include_endpoints: static_df = static_sequence(eos, e_min=e_min, e_max=e_max, n=endpoint_n, grid=grid) kepler_df = kepler_sequence(eos, e_min=e_min, e_max=e_max, n=endpoint_n, grid=grid, fix_error=fix_error) srow = _nearest_m0_row(static_df, M0) krow = _nearest_m0_row(kepler_df, M0) if srow is not None: srow["endpoint"] = "static" rows.append(srow) if krow is not None: krow["endpoint"] = "kepler" rows.append(krow) if not rows: return pd.DataFrame() out = pd.DataFrame(rows) if "rho_c" in out.columns: out = out.sort_values("rho_c") return out.reset_index(drop=True) def constant_J_sequence( eos: str, J: float, e_min: float, e_max: float, n: int = 20, grid: str = "std", fix_error: float = 1e-4, ) -> pd.DataFrame: """Constant angular-momentum sequence.""" df = run_rns(task="jmoment", eos=eos, J=J, e=e_min, e_max=e_max, n=n, fix_error=fix_error, grid=grid).table return df.dropna(subset=["M"]).reset_index(drop=True) # ------------------------------------------------------------------ # EOS file reader # ------------------------------------------------------------------ def load_eos(eos: str) -> pd.DataFrame: """ Read a tabulated EOS file in RNS format. Columns: epsilon [g/cm^3], P [dyn/cm^2], H [cm^2/s^2], n_B [1/cm^3]. """ path = _eos_path(eos) with open(path) as f: first = f.readline().strip() n_pts = int(first.split()[0]) arr = np.loadtxt(f, max_rows=n_pts) return pd.DataFrame( arr, columns=["epsilon", "P", "H", "n_B"] ) # ------------------------------------------------------------------ # Full metric grid (-p 3) # ------------------------------------------------------------------ def parse_metric_grid(text: str) -> pd.DataFrame: """ Parse the -p 3 grid output. RNS prints 11 columns per data line: r(M) s cos(theta) rho gamma alpha omega 2M/r 2M/r nu_K nu_K We keep: s (col 1), mu=cos(theta) (col 2), and the four metric potentials rho/gamma/alpha/omega (cols 3-6). """ rows = [] for line in text.splitlines(): toks = line.split() if len(toks) != 11: continue try: vals = [float(t) for t in toks] except ValueError: continue # cols: 0=r(M), 1=s, 2=mu, 3=rho, 4=gamma, 5=alpha, 6=omega rows.append([vals[1], vals[2], vals[3], vals[4], vals[5], vals[6]]) return pd.DataFrame( rows, columns=["s", "mu", "rho_pot", "gamma_pot", "alpha_pot", "omega_pot"] ) # ------------------------------------------------------------------ # Small style helpers for matplotlib so the notebook keeps a consistent look # ------------------------------------------------------------------ def apply_paper_style(): import matplotlib as mpl mpl.rcParams.update({ "figure.figsize": (5.5, 4.0), "figure.dpi": 110, "savefig.dpi": 140, "font.size": 11, "axes.titlesize": 12, "axes.labelsize": 11, "legend.fontsize": 9, "lines.linewidth": 1.4, "axes.grid": True, "grid.alpha": 0.3, "grid.linestyle":"--", })