feat: line-profile check of the RM closed forms and its figure for the paper appendix (7de93d64) · Commits · Gaetano Scandariato / PyAETNA

paper_benchmarks/make_figures.py

+41 −0

Changes for paper_benchmarks/make_figures.py: 41 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		@@ -532,6 +532,45 @@ def eiv_validation() -> None:
		print(f"wrote {out.relative_to(PAPER)}")


		def rm_line_profile() -> None:
		"""The three RM observables against brute-force line-profile simulations (Appendix B)."""
		d = json.loads(
		(
		BENCH / "toi5398_rm_pyorbit" / "results" / "line_profile_check.json"
		).read_text()
		)
		x = np.array(d["x_planet"])
		r = d["beta_kms"]["1.81"]
		fig, ax = plt.subplots(figsize=(COLUMN_WIDTH, 2.3))
		for key, color, marker, label in (
		("ccf", BLUE, "o", "Gaussian fit (Bou\u00e9 et al. 2013)"),
		("template", ORANGE, "s", "template (Hirano et al. 2011)"),
		):
		ax.plot(x, r[key], color=color, lw=1.5, label=label)
		ax.plot(
		x,
		r[f"{key}_ref"],
		marker,
		color=color,
		ms=3,
		mec="white",
		mew=0.5,
		ls="none",
		)
		ax.plot(x, r["flux_weighted"], color=AQUA, lw=1.5, ls="--", label="flux-weighted")
		ax.axhline(0, color=GRID, lw=0.8, zorder=0)
		ax.set_xlabel("planet position along the chord [stellar radii]")
		ax.set_ylabel(r"RM anomaly [m s$^{-1}$]")
		ax.legend(frameon=False, loc="lower left", fontsize=6, handlelength=1.5)
		ax.set_ylim(-75, 60)
		fig.tight_layout()
		out = FIG / "appendix" / "rm_observables.pdf"
		out.parent.mkdir(parents=True, exist_ok=True)
		fig.savefig(out)
		plt.close(fig)
		print(f"wrote {out.relative_to(PAPER)}")


		def main() -> None:
		branches()
		ns_calibration()
		@@ -546,6 +585,8 @@ def main() -> None:
		eiv_validation()
		if (BENCH / "scaling" / "results" / "summary.json").exists():
		scaling()
		if (BENCH / "toi5398_rm_pyorbit" / "results" / "line_profile_check.json").exists():
		rm_line_profile()


		if __name__ == "__main__":

paper_benchmarks/toi5398_rm_pyorbit/line_profile_check.py

0 → 100644

+98 −0

Changes for paper_benchmarks/toi5398_rm_pyorbit/line_profile_check.py: 98 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		"""The closed-form RM observables against brute-force simulations of the line profile (Appendix B).

		The disc of a TOI-5398-like star (v sin i = 6.79 km/s, quadratic limb darkening u1 = 0.58, u2 = 0.17) is divided
		into 801 x 801 cells, each contributing a Gaussian line shifted by its rotational velocity; a planet of radius ratio
		0.0905 crosses it at b = 0.36 (aligned orbit). At each position the in-transit line gives the RV:
		- ``ccf``: the centre of a Gaussian fitted to the line (HARPS-type pipelines), against Boue et al. (2013, Eq. 14);
		- ``template``: the shift that maximises the cross-correlation with the out-of-transit line, against Hirano et al.
		(2011, Eq. 16);
		- ``flux_weighted``: the first moment of the line, -f v_p / (1 - f), exact for the disc integration.
		The simulations are those of ``tests/unit/core/test_rm_ccf.py``. Writes results/line_profile_check.json
		(the curves and the largest difference in units of the amplitude, for three line widths beta).

		usage: python line_profile_check.py
		"""

		from __future__ import annotations

		import json
		import math
		import sys
		from pathlib import Path

		import jax

		jax.config.update("jax_enable_x64", True)
		import jax.numpy as jnp
		import numpy as np

		HERE = Path(__file__).resolve().parent
		sys.path.insert(0, str(HERE.parent.parent / "tests" / "unit" / "core"))
		from test_rm_ccf import (
		U1,
		U2,
		VSINI,
		K,
		_brute_force,
		_brute_force_gaussian_fit,
		)

		from PyAETNA.body_models import _boue2013_anomaly, _hirano2011_anomaly

		B_IMPACT = 0.36
		XP = np.round(np.arange(-0.9, 0.901, 0.05), 3) # planet centre on the disc (\|x\| < 0.93)
		BETAS = (1.81, 3.0, 6.0) # km/s; 1.81 is the HARPS-N line width used for TOI-5398


		def one_beta(beta: float) -> dict:
		sigma = beta / math.sqrt(2)
		rows = {
		k: [] for k in ("ccf", "ccf_ref", "template", "template_ref", "flux_weighted")
		}
		for xp in XP:
		g = _brute_force_gaussian_fit(float(xp), B_IMPACT, sigma)
		c = _brute_force(float(xp), B_IMPACT, beta)
		f, vp = jnp.array([g["f"]]), jnp.array([g["v_p"]])
		rows["ccf"].append(float(_boue2013_anomaly(f, vp, VSINI, sigma, K, U1, U2)[0]))
		rows["ccf_ref"].append(float(g["rv"]))
		rows["template"].append(
		float(
		_hirano2011_anomaly(
		jnp.array([c["f"]]), jnp.array([c["v_p"]]), jnp.array([c["cos_theta"]]),
		VSINI, beta, 0.0, 0.0, U1, U2,
		)[0]
		)
		) # fmt: skip
		rows["template_ref"].append(float(c["rv"]))
		rows["flux_weighted"].append(-g["f"] * g["v_p"] / (1 - g["f"]))
		out = {k: [round(v * 1e3, 3) for v in vals] for k, vals in rows.items()} # m/s
		for obs in ("ccf", "template"):
		ref = np.array(rows[f"{obs}_ref"])
		diff = np.abs(np.array(rows[obs]) - ref)
		out[f"{obs}_max_abs_diff_over_amplitude"] = float(
		diff.max() / np.abs(ref).max()
		)
		return out


		def main() -> None:
		result = {
		"star": {"v_sin_i_kms": VSINI, "u1": U1, "u2": U2},
		"planet": {"rp_by_rs": K, "impact_parameter": B_IMPACT, "lambda_deg": 0.0},
		"grid_cells": 801,
		"x_planet": XP.tolist(),
		"units": "m/s",
		"beta_kms": {str(b): one_beta(b) for b in BETAS},
		}
		for b, r in result["beta_kms"].items():
		print(
		f"beta {b}: ccf {100 * r['ccf_max_abs_diff_over_amplitude']:.1f} %, "
		f"template {100 * r['template_max_abs_diff_over_amplitude']:.1f} % of the amplitude"
		)
		(HERE / "results" / "line_profile_check.json").write_text(
		json.dumps(result, indent=1)
		)


		if __name__ == "__main__":
		main()

paper_benchmarks/toi5398_rm_pyorbit/results/line_profile_check.json

0 → 100644

+652 −0

Changes for paper_benchmarks/toi5398_rm_pyorbit/results/line_profile_check.json: 652 added lines, 0 removed lines.
Original line number	Diff line number	Diff line
		{
		"star": {
		"v_sin_i_kms": 6.79,
		"u1": 0.58,
		"u2": 0.17
		},
		"planet": {
		"rp_by_rs": 0.0905,
		"impact_parameter": 0.36,
		"lambda_deg": 0.0
		},
		"grid_cells": 801,
		"x_planet": [
		-0.9,
		-0.85,
		-0.8,
		-0.75,
		-0.7,
		-0.65,
		-0.6,
		-0.55,
		-0.5,
		-0.45,
		-0.4,
		-0.35,
		-0.3,
		-0.25,
		-0.2,
		-0.15,
		-0.1,
		-0.05,
		0.0,
		0.05,
		0.1,
		0.15,
		0.2,
		0.25,
		0.3,
		0.35,
		0.4,
		0.45,
		0.5,
		0.55,
		0.6,
		0.65,
		0.7,
		0.75,
		0.8,
		0.85,
		0.9
		],
		"units": "m/s",
		"beta_kms": {
		"1.81": {
		"ccf": [
		17.901,
		28.326,
		35.068,
		39.949,
		43.93,
		47.048,
		49.26,
		50.501,
		50.71,
		49.841,
		47.872,
		44.812,
		40.697,
		35.6,
		29.625,
		22.907,
		15.606,
		7.905,
		-0.0,
		-7.905,
		-15.606,
		-22.907,
		-29.625,
		-35.6,
		-40.697,
		-44.812,
		-47.872,
		-49.841,
		-50.71,
		-50.501,
		-49.26,
		-47.048,
		-43.93,
		-39.949,
		-35.068,
		-28.326,
		-17.901
		],
		"ccf_ref": [
		17.819,
		28.128,
		34.841,
		39.765,
		43.829,
		47.057,
		49.391,
		50.751,
		51.064,
		50.273,
		48.35,
		45.299,
		41.163,
		36.017,
		29.975,
		23.176,
		15.788,
		7.997,
		-0.002,
		-7.997,
		-15.788,
		-23.176,
		-29.975,
		-36.017,
		-41.163,
		-45.299,
		-48.35,
		-50.273,
		-51.064,
		-50.751,
		-49.391,
		-47.057,
		-43.828,
		-39.765,
		-34.841,
		-28.128,
		-17.819
		],
		"template": [
		20.852,
		32.461,
		39.047,
		42.747,
		44.798,
		45.427,
		44.824,
		43.189,
		40.726,
		37.629,
		34.069,
		30.187,
		26.092,
		21.86,
		17.548,
		13.188,
		8.802,
		4.404,
		-0.0,
		-4.404,
		-8.802,
		-13.188,
		-17.548,
		-21.86,
		-26.092,
		-30.187,
		-34.069,
		-37.629,
		-40.726,
		-43.189,
		-44.824,
		-45.427,
		-44.798,
		-42.747,
		-39.047,
		-32.461,
		-20.852
		],
		"template_ref": [
		20.527,
		31.546,
		39.101,
		41.68,
		43.127,
		46.596,
		43.211,
		42.108,
		40.424,
		38.29,
		32.824,
		30.125,
		27.269,
		21.313,
		18.295,
		12.242,
		9.169,
		3.087,
		0.0,
		-3.087,
		-9.169,
		-12.242,
		-18.295,
		-21.313,
		-27.269,
		-30.125,
		-32.824,
		-38.29,
		-40.424,
		-42.108,
		-43.211,
		-46.596,
		-43.127,
		-41.68,
		-39.101,
		-31.546,
		-20.527
		],
		"flux_weighted": [
		22.403,
		32.843,
		36.736,
		37.715,
		37.578,
		36.693,
		35.257,
		33.393,
		31.189,
		28.709,
		26.002,
		23.11,
		20.065,
		16.896,
		13.629,
		10.286,
		6.888,
		3.453,
		-0.0,
		-3.453,
		-6.888,
		-10.286,
		-13.629,
		-16.896,
		-20.065,
		-23.11,
		-26.002,
		-28.709,
		-31.189,
		-33.393,
		-35.257,
		-36.693,
		-37.578,
		-37.715,
		-36.736,
		-32.843,
		-22.403
		],
		"ccf_max_abs_diff_over_amplitude": 0.009547715643887672,
		"template_max_abs_diff_over_amplitude": 0.03586314808606435
		},
		"3.0": {
		"ccf": [
		19.075,
		29.642,
		35.837,
		39.848,
		42.825,
		44.891,
		46.073,
		46.373,
		45.789,
		44.324,
		41.997,
		38.842,
		34.909,
		30.268,
		25.007,
		19.228,
		13.047,
		6.593,
		-0.0,
		-6.593,
		-13.047,
		-19.228,
		-25.007,
		-30.268,
		-34.909,
		-38.842,
		-41.997,
		-44.324,
		-45.789,
		-46.373,
		-46.073,
		-44.891,
		-42.825,
		-39.848,
		-35.837,
		-29.642,
		-19.075
		],
		"ccf_ref": [
		19.039,
		29.537,
		35.731,
		39.791,
		42.84,
		44.991,
		46.261,
		46.644,
		46.13,
		44.717,
		42.419,
		39.27,
		35.321,
		30.643,
		25.327,
		19.48,
		13.221,
		6.681,
		0.0,
		-6.681,
		-13.221,
		-19.48,
		-25.327,
		-30.643,
		-35.321,
		-39.27,
		-42.419,
		-44.717,
		-46.13,
		-46.644,
		-46.261,
		-44.991,
		-42.84,
		-39.791,
		-35.731,
		-29.537,
		-19.039
		],
		"template": [
		20.001,
		30.84,
		36.875,
		40.465,
		42.86,
		44.235,
		44.668,
		44.218,
		42.938,
		40.887,
		38.133,
		34.746,
		30.805,
		26.388,
		21.577,
		16.454,
		11.098,
		5.587,
		-0.0,
		-5.587,
		-11.098,
		-16.454,
		-21.577,
		-26.388,
		-30.805,
		-34.746,
		-38.133,
		-40.887,
		-42.938,
		-44.218,
		-44.668,
		-44.235,
		-42.86,
		-40.465,
		-36.875,
		-30.84,
		-20.001
		],
		"template_ref": [
		19.968,
		30.505,
		37.698,
		40.195,
		41.861,
		42.818,
		43.12,
		42.809,
		41.92,
		40.495,
		38.579,
		33.222,
		30.478,
		27.401,
		21.049,
		17.477,
		10.741,
		6.898,
		0.0,
		-6.898,
		-10.741,
		-17.477,
		-21.049,
		-27.401,
		-30.478,
		-33.222,
		-38.579,
		-40.495,
		-41.92,
		-42.809,
		-43.12,
		-42.818,
		-41.861,
		-40.195,
		-37.698,
		-30.505,
		-19.968
		],
		"flux_weighted": [
		22.403,
		32.843,
		36.736,
		37.715,
		37.578,
		36.693,
		35.257,
		33.393,
		31.189,
		28.709,
		26.002,
		23.11,
		20.065,
		16.896,
		13.629,
		10.286,
		6.888,
		3.453,
		-0.0,
		-3.453,
		-6.888,
		-10.286,
		-13.629,
		-16.896,
		-20.065,
		-23.11,
		-26.002,
		-28.709,
		-31.189,
		-33.393,
		-35.257,
		-36.693,
		-37.578,
		-37.715,
		-36.736,
		-32.843,
		-22.403
		],
		"ccf_max_abs_diff_over_amplitude": 0.00918179253629234,
		"template_max_abs_diff_over_amplitude": 0.0358943649463406
		},
		"6.0": {
		"ccf": [
		20.984,
		31.557,
		36.552,
		38.903,
		40.111,
		40.446,
		40.041,
		38.984,
		37.341,
		35.166,
		32.509,
		29.42,
		25.947,
		22.142,
		18.056,
		13.743,
		9.258,
		4.658,
		-0.0,
		-4.658,
		-9.258,
		-13.743,
		-18.056,
		-22.142,
		-25.947,
		-29.42,
		-32.509,
		-35.166,
		-37.341,
		-38.984,
		-40.041,
		-40.446,
		-40.111,
		-38.903,
		-36.552,
		-31.557,
		-20.984
		],
		"ccf_ref": [
		21.021,
		31.614,
		36.646,
		39.042,
		40.297,
		40.674,
		40.308,
		39.281,
		37.659,
		35.494,
		32.837,
		29.736,
		26.241,
		22.404,
		18.277,
		13.915,
		9.376,
		4.718,
		0.0,
		-4.718,
		-9.376,
		-13.915,
		-18.277,
		-22.404,
		-26.241,
		-29.736,
		-32.837,
		-35.494,
		-37.659,
		-39.281,
		-40.308,
		-40.674,
		-40.297,
		-39.042,
		-36.646,
		-31.614,
		-21.021
		],
		"template": [
		21.089,
		31.725,
		36.751,
		39.119,
		40.342,
		40.686,
		40.289,
		39.236,
		37.592,
		35.413,
		32.746,
		29.642,
		26.15,
		22.32,
		18.204,
		13.858,
		9.337,
		4.698,
		-0.0,
		-4.698,
		-9.337,
		-13.858,
		-18.204,
		-22.32,
		-26.15,
		-29.642,
		-32.746,
		-35.413,
		-37.592,
		-39.236,
		-40.289,
		-40.686,
		-40.342,
		-39.119,
		-36.751,
		-31.725,
		-21.089
		],
		"template_ref": [
		20.744,
		31.162,
		37.665,
		39.317,
		40.168,
		40.407,
		40.128,
		39.391,
		38.243,
		36.72,
		31.858,
		29.69,
		27.252,
		21.579,
		18.706,
		12.672,
		9.516,
		3.279,
		0.0,
		-3.279,
		-9.516,
		-12.672,
		-18.706,
		-21.579,
		-27.252,
		-29.69,
		-31.858,
		-36.72,
		-38.243,
		-39.391,
		-40.128,
		-40.407,
		-40.168,
		-39.317,
		-37.665,
		-31.162,
		-20.744
		],
		"flux_weighted": [
		22.403,
		32.843,
		36.736,
		37.715,
		37.578,
		36.693,
		35.257,
		33.393,
		31.189,
		28.709,
		26.002,
		23.11,
		20.065,
		16.896,
		13.629,
		10.286,
		6.888,
		3.453,
		-0.0,
		-3.453,
		-6.888,
		-10.286,
		-13.629,
		-16.896,
		-20.065,
		-23.11,
		-26.002,
		-28.709,
		-31.189,
		-33.393,
		-35.257,
		-36.693,
		-37.578,
		-37.715,
		-36.736,
		-32.843,
		-22.403
		],
		"ccf_max_abs_diff_over_amplitude": 0.008074869614606014,
		"template_max_abs_diff_over_amplitude": 0.03511893144059537
		}
		}
		}
		No newline at end of file