{ "cells": [ { "cell_type": "markdown", "id": "22eb2263", "metadata": {}, "source": [ "# Coordinates and TwoPunctures from EOB IDs\n", "\n", "This tutorial has two goals:\n", "1. illustrate the main utilities in `PyART.analytic.CoordsChange`;\n", "2. show how to start from EOB initial data, map them to ADM with `eob_ID_to_ADM`, and assemble a TwoPunctures-ready initial-data dictionary.\n", "\n", "The low-level coordinate transforms are useful for inspection and sanity checks. For the actual EOB-to-ADM handoff, the main function to use is `eob_ID_to_ADM`." ] }, { "cell_type": "code", "execution_count": null, "id": "590351ee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TEOB available = True\n" ] } ], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'retina'\n", "\n", "from pathlib import Path\n", "from tempfile import TemporaryDirectory\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pprint import pprint\n", "\n", "from PyART.analytic import CoordsChange\n", "import PyART.analytic.coordschange as cc\n", "from PyART.numerical.twopuncts_ID import TwoPunctID\n", "\n", "try:\n", " import PyART.models.teob as teob\n", " HAVE_TEOB = True\n", " TEOB_IMPORT_ERROR = None\n", "except Exception as exc:\n", " HAVE_TEOB = False\n", " TEOB_IMPORT_ERROR = exc\n", "\n", "np.set_printoptions(precision=10, suppress=True)\n", "\n", "print(f\"TEOB available = {HAVE_TEOB}\")\n", "if TEOB_IMPORT_ERROR is not None:\n", " print(f\"TEOB import error: {TEOB_IMPORT_ERROR}\")\n", " print(\"To run this tutorial, install TEOBResumS and ensure it is in your PYTHONPATH.\")" ] }, { "cell_type": "markdown", "id": "fe128183", "metadata": {}, "source": [ "## 1. Low-level coordinate changes\n", "\n", "`CoordsChange` exposes polar/cartesian transforms and the 2PN EOB/ADM maps. The small example below performs both round trips and reports the residuals." ] }, { "cell_type": "code", "execution_count": 2, "id": "c6b97087", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "polar -> cartesian\n", "[11.4640378695 3.5462424799 -0.1323894757 0.3603019167]\n", "\n", "cartesian -> polar residuals\n", "[ 0. -0. -0. 0.]\n", "\n", "EOB -> ADM\n", "q_ADM = [10.5861678003 3.2569829327]\n", "p_ADM = [-0.1413954526 0.3907243109]\n", "\n", "round-trip residuals\n", "max |qe_back - qe| = 0.00652476010515457\n", "max |pe_back - pe| = 0.000497622413047405\n" ] } ], "source": [ "r0, phi0, pr0, pphi0 = 12.0, 0.3, -0.02, 4.6\n", "x, y, px, py = CoordsChange.Polar2Cartesian(r0, phi0, pr0, pphi0)\n", "r1, phi1, pr1, pphi1 = CoordsChange.Cartesian2Polar(x, y, px, py)\n", "\n", "q = 1.2\n", "nu = q / (1.0 + q) ** 2\n", "qe = np.array([x, y])\n", "pe = np.array([px, py])\n", "\n", "qa, pa = CoordsChange.Eob2Adm(qe, pe, nu, PN_order=2)\n", "qe_back, pe_back = CoordsChange.Adm2Eob(qa, pa, nu, PN_order=2)\n", "\n", "print(\"polar -> cartesian\")\n", "print(np.array([x, y, px, py]))\n", "print()\n", "print(\"cartesian -> polar residuals\")\n", "print(np.array([r1 - r0, phi1 - phi0, pr1 - pr0, pphi1 - pphi0]))\n", "print()\n", "print(\"EOB -> ADM\")\n", "print(\"q_ADM =\", qa)\n", "print(\"p_ADM =\", pa)\n", "print()\n", "print(\"round-trip residuals\")\n", "print(\"max |qe_back - qe| =\", np.max(np.abs(qe_back - qe)))\n", "print(\"max |pe_back - pe| =\", np.max(np.abs(pe_back - pe)))" ] }, { "cell_type": "markdown", "id": "c45502ca", "metadata": {}, "source": [ "## 2. Reusable EOB and TwoPunctures helpers\n", "\n", "The functions do the following:\n", "\n", "- `ADM_from_EOB(...)` builds a TEOB waveform, plots the $(2,2)$ amplitude, computes $t_{\\rm mrg} = -u_0$, and converts the initial state to ADM with `eob_ID_to_ADM`.\n", "- `TP_from_ADM(...)` maps that ADM dictionary into the TwoPunctures-style initial-data dictionary.\n", "- `TP_diagnostics_from_ADM(...)` keeps a few geometric diagnostics explicit, including the center offset.\n", "- `TP_parfile_from_ADM(...)` writes a TwoPunctures parfile from the same `adm_id` data when the configuration is compatible with the current `TwoPunctID` helper." ] }, { "cell_type": "code", "execution_count": 3, "id": "fe4b4729", "metadata": {}, "outputs": [], "source": [ "def ADM_from_EOB(\n", " M=1.0,\n", " q=1.0,\n", " chi1=(0.0, 0.0, 0.0),\n", " chi2=(0.0, 0.0, 0.0),\n", " f0=0.001,\n", " e0=0.1,\n", " plot_waveform=True,\n", " verbose=True,\n", "):\n", " \"\"\"Initial data based on EOB prediction, mapped to ADM.\"\"\"\n", " if not HAVE_TEOB:\n", " raise RuntimeError(\n", " \"TEOBResumS is not available in this environment, so ADM_from_EOB cannot run.\"\n", " )\n", "\n", " pars = teob.CreateDict(\n", " M=M,\n", " q=q,\n", " chi1x=chi1[0],\n", " chi1y=chi1[1],\n", " chi1z=chi1[2],\n", " chi2x=chi2[0],\n", " chi2y=chi2[1],\n", " chi2z=chi2[2],\n", " f0=f0,\n", " ecc=e0,\n", " )\n", " eob_wave = teob.Waveform_EOB(pars=pars)\n", "\n", " if plot_waveform:\n", " fig, ax = plt.subplots(figsize=(8, 4))\n", " ax.plot(eob_wave.u, eob_wave.hlm[(2, 2)][\"A\"])\n", " ax.set_xlabel(\"t / M\")\n", " ax.set_ylabel(r\"$|h_{22}|$\")\n", " ax.set_title(\"TEOB waveform amplitude\")\n", " ax.grid(True, alpha=0.3)\n", " plt.show()\n", "\n", " t_mrg = -eob_wave.u[0]\n", "\n", " adm_redundant_id = cc.eob_ID_to_ADM(eob_wave, verbose=verbose)\n", " adm_id = {\n", " \"px\": adm_redundant_id[\"px\"],\n", " \"py\": adm_redundant_id[\"py\"],\n", " \"D\": adm_redundant_id[\"D\"],\n", " \"M\": M,\n", " \"q\": q,\n", " \"chi1x\": chi1[0],\n", " \"chi1y\": chi1[1],\n", " \"chi1z\": chi1[2],\n", " \"chi2x\": chi2[0],\n", " \"chi2y\": chi2[1],\n", " \"chi2z\": chi2[2],\n", " \"x1\": adm_redundant_id[\"x1\"],\n", " \"x2\": adm_redundant_id[\"x2\"],\n", " \"x_offset\": adm_redundant_id[\"x_offset\"],\n", " }\n", " return adm_id, t_mrg\n", "\n", "\n", "def TP_from_ADM(adm_id):\n", " \"\"\"Build a TwoPunctures-style initial-data dictionary from ADM IDs.\"\"\"\n", " D = adm_id[\"D\"]\n", " M = adm_id[\"M\"]\n", " q = adm_id[\"q\"]\n", "\n", " assert q >= 1.0\n", "\n", " mp = M * q / (1 + q)\n", " mm = M / (1 + q)\n", "\n", " xp = D * mm\n", " xm = -D * mp\n", "\n", " Spx = adm_id[\"chi1x\"] * mp**2\n", " Spy = adm_id[\"chi1y\"] * mp**2\n", " Spz = adm_id[\"chi1z\"] * mp**2\n", "\n", " Smx = adm_id[\"chi2x\"] * mm**2\n", " Smy = adm_id[\"chi2y\"] * mm**2\n", " Smz = adm_id[\"chi2z\"] * mm**2\n", "\n", " Ppx = adm_id[\"px\"]\n", " Ppy = adm_id[\"py\"]\n", " Ppz = 0.0\n", "\n", " Pmx = -Ppx\n", " Pmy = -Ppy\n", " Pmz = -Ppz\n", "\n", " initial_data = {\n", " \"bh_0_m\": mp,\n", " \"bh_0_madm\": mp,\n", " \"bh_0_x\": xp,\n", " \"bh_0_px\": Ppx,\n", " \"bh_0_py\": Ppy,\n", " \"bh_0_sx\": Spx,\n", " \"bh_0_sy\": Spy,\n", " \"bh_0_sz\": Spz,\n", " \"bh_1_m\": mm,\n", " \"bh_1_madm\": mm,\n", " \"bh_1_x\": xm,\n", " \"bh_1_px\": Pmx,\n", " \"bh_1_py\": Pmy,\n", " \"bh_1_sx\": Smx,\n", " \"bh_1_sy\": Smy,\n", " \"bh_1_sz\": Smz,\n", " \"give_bare_mass\": False,\n", " \"bitant\": False,\n", " }\n", " return initial_data\n", "\n", "\n", "def TP_diagnostics_from_ADM(adm_id):\n", " \"\"\"Expose a few geometric quantities that are useful for checks.\"\"\"\n", " D = adm_id[\"D\"]\n", " M = adm_id[\"M\"]\n", " q = adm_id[\"q\"]\n", "\n", " mp = M * q / (1 + q)\n", " mm = M / (1 + q)\n", " xp = D * mm\n", " xm = -D * mp\n", " half_D = D * M / 2.0\n", " center_offset = xp - half_D\n", "\n", " return {\n", " \"xp\": xp,\n", " \"xm\": xm,\n", " \"half_D\": half_D,\n", " \"center_offset\": center_offset,\n", " }" ] }, { "cell_type": "code", "execution_count": 4, "id": "2e5f1806", "metadata": {}, "outputs": [], "source": [ "def TP_parfile_from_ADM(\n", " adm_id,\n", " outdir,\n", " parfile_name=\"from_eob_ids.par\",\n", " npoints_A=8,\n", " npoints_B=8,\n", " npoints_phi=6,\n", " verbose=False,\n", "):\n", " \"\"\"Write a TwoPunctures parfile from ADM IDs using TwoPunctID.\"\"\"\n", " if abs(adm_id[\"chi1x\"]) > 0 or abs(adm_id[\"chi1y\"]) > 0:\n", " raise ValueError(\n", " \"TwoPunctID currently supports only aligned spins in the parfile helper.\"\n", " )\n", " if abs(adm_id[\"chi2x\"]) > 0 or abs(adm_id[\"chi2y\"]) > 0:\n", " raise ValueError(\n", " \"TwoPunctID currently supports only aligned spins in the parfile helper.\"\n", " )\n", "\n", " D = float(adm_id[\"D\"])\n", " px = float(adm_id[\"px\"])\n", " py = float(adm_id[\"py\"])\n", " P = float(np.hypot(px, py))\n", " L = float(D * py)\n", "\n", " tp = TwoPunctID(\n", " q=float(adm_id[\"q\"]),\n", " D=D,\n", " L=L,\n", " chi1z=float(adm_id[\"chi1z\"]),\n", " chi2z=float(adm_id[\"chi2z\"]),\n", " npoints_A=npoints_A,\n", " npoints_B=npoints_B,\n", " npoints_phi=npoints_phi,\n", " outdir=outdir,\n", " verbose=verbose,\n", " )\n", " parfile = tp.create_TP_parfile(P, parfile=parfile_name)\n", " return Path(outdir) / parfile" ] }, { "cell_type": "markdown", "id": "fac569ce", "metadata": {}, "source": [ "## 3. Run the workflow on a concrete configuration\n", "\n", "We now reproduce the same flow as in your driver script: choose intrinsic parameters, call `ADM_from_EOB`, convert them with `TP_from_ADM`, and inspect the final TwoPunctures dictionary.\n", "\n", "The cell also checks a few consistency relations, including the agreement between the geometric center offset and the `x_offset` returned by `eob_ID_to_ADM`." ] }, { "cell_type": "code", "execution_count": 5, "id": "b05d0bfb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ADM IDs\n", "D = 10.00583246818054\n", "px = 0.0\n", "py = 0.0828681049253926\n", "x1 = 3.3352774893935133\n", "x2 = -6.670554978787027\n", "x_offset = -1.6676387446967564\n", "t_mrg = 625.3477422287092\n", "\n", "TwoPunctures initial_data\n", "bh_0_m = 0.6666666666666666\n", "bh_0_x = 3.335277489393513\n", "bh_0_px = 0.0\n", "bh_0_py = 0.0828681049253926\n", "bh_0_sz = 0.0\n", "bh_1_m = 0.3333333333333333\n", "bh_1_x = -6.670554978787026\n", "bh_1_px = -0.0\n", "bh_1_py = -0.0828681049253926\n", "bh_1_sz = 0.0\n", "give_bare_mass = False\n", "bitant = False\n", "\n", "Derived diagnostics\n", "xp = 3.335277489393513\n", "xm = -6.670554978787026\n", "half_D = 5.00291623409027\n", "center_offset = -1.6676387446967569\n", "\n", "Full initial_data dictionary\n", "{'bh_0_m': 0.6666666666666666,\n", " 'bh_0_madm': 0.6666666666666666,\n", " 'bh_0_px': np.float64(0.0),\n", " 'bh_0_py': np.float64(0.0828681049253926),\n", " 'bh_0_sx': 0.0,\n", " 'bh_0_sy': 0.0,\n", " 'bh_0_sz': 0.0,\n", " 'bh_0_x': np.float64(3.335277489393513),\n", " 'bh_1_m': 0.3333333333333333,\n", " 'bh_1_madm': 0.3333333333333333,\n", " 'bh_1_px': np.float64(-0.0),\n", " 'bh_1_py': np.float64(-0.0828681049253926),\n", " 'bh_1_sx': 0.0,\n", " 'bh_1_sy': 0.0,\n", " 'bh_1_sz': 0.0,\n", " 'bh_1_x': np.float64(-6.670554978787026),\n", " 'bitant': False,\n", " 'give_bare_mass': False}\n" ] } ], "source": [ "intrinsic = {\n", " \"test\": {\n", " \"q\": 2.0,\n", " \"chi1x\": 0.0,\n", " \"chi1y\": 0.0,\n", " \"chi1z\": 0.0,\n", " \"chi2x\": 0.0,\n", " \"chi2y\": 0.0,\n", " \"chi2z\": 0.0,\n", " \"f0\": 0.01006,\n", " \"ecc\": 0.1,\n", " }\n", "}\n", "\n", "ID = \"test\"\n", "adm_id, tmax = ADM_from_EOB(\n", " M=1.0,\n", " q=intrinsic[ID][\"q\"],\n", " chi1=[\n", " intrinsic[ID][\"chi1x\"],\n", " intrinsic[ID][\"chi1y\"],\n", " intrinsic[ID][\"chi1z\"],\n", " ],\n", " chi2=[\n", " intrinsic[ID][\"chi2x\"],\n", " intrinsic[ID][\"chi2y\"],\n", " intrinsic[ID][\"chi2z\"],\n", " ],\n", " f0=intrinsic[ID][\"f0\"],\n", " e0=intrinsic[ID][\"ecc\"],\n", " plot_waveform=False,\n", " verbose=False,\n", " )\n", "\n", "initial_data = TP_from_ADM(adm_id)\n", "tp_diagnostics = TP_diagnostics_from_ADM(adm_id)\n", "\n", "assert np.isclose(initial_data[\"bh_0_x\"] - initial_data[\"bh_1_x\"], adm_id[\"D\"])\n", "assert np.isclose(initial_data[\"bh_0_px\"], -initial_data[\"bh_1_px\"])\n", "assert np.isclose(initial_data[\"bh_0_py\"], -initial_data[\"bh_1_py\"])\n", "assert np.isclose(tp_diagnostics[\"center_offset\"], adm_id[\"x_offset\"])\n", "\n", "print(\"ADM IDs\")\n", "for key in (\"D\", \"px\", \"py\", \"x1\", \"x2\", \"x_offset\"):\n", " print(f\"{key:10s} = {adm_id[key]}\")\n", "print(f\"{'t_mrg':10s} = {tmax}\")\n", "\n", "print(\"\\nTwoPunctures initial_data\")\n", "for key in (\n", " \"bh_0_m\",\n", " \"bh_0_x\",\n", " \"bh_0_px\",\n", " \"bh_0_py\",\n", " \"bh_0_sz\",\n", " \"bh_1_m\",\n", " \"bh_1_x\",\n", " \"bh_1_px\",\n", " \"bh_1_py\",\n", " \"bh_1_sz\",\n", " \"give_bare_mass\",\n", " \"bitant\",\n", "):\n", " print(f\"{key:14s} = {initial_data[key]}\")\n", "\n", "print(\"\\nDerived diagnostics\")\n", "for key, value in tp_diagnostics.items():\n", " print(f\"{key:14s} = {value}\")\n", "\n", "print(\"\\nFull initial_data dictionary\")\n", "pprint(initial_data)" ] }, { "cell_type": "markdown", "id": "3beb135c", "metadata": {}, "source": [ "## 4. Create a TwoPunctures parfile from the same ADM IDs\n", "\n", "The `initial_data` dictionary is useful when you want the puncture data in Python. In addition, the same `adm_id` can be turned into a TwoPunctures parfile through `TwoPunctID.create_TP_parfile`.\n", "\n", "For the current `TwoPunctID` helper, this parfile path is restricted to aligned-spin configurations, so the helper raises if the in-plane spin components are nonzero." ] }, { "cell_type": "code", "execution_count": 6, "id": "94b93e8a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "parfile written to: /tmp/tmpwgk547lu/from_eob_ids.par\n", "\n", "center_offset1 actual = -1.667638744697e+00 expected = -1.667638744697e+00\n", "par_P_minus1 actual = 0.000000000000e+00 expected = -0.000000000000e+00\n", "par_P_minus2 actual = -8.286810492539e-02 expected = -8.286810492539e-02\n", "par_P_plus1 actual = -0.000000000000e+00 expected = 0.000000000000e+00\n", "par_P_plus2 actual = 8.286810492539e-02 expected = 8.286810492539e-02\n", "par_b actual = 5.002916234090e+00 expected = 5.002916234090e+00\n" ] } ], "source": [ "with TemporaryDirectory() as tmpdir:\n", " parfile_path = TP_parfile_from_ADM(adm_id, outdir=tmpdir)\n", "\n", " selected_keys = {\n", " \"par_b\",\n", " \"par_P_plus1\",\n", " \"par_P_plus2\",\n", " \"par_P_minus1\",\n", " \"par_P_minus2\",\n", " \"center_offset1\",\n", " }\n", "\n", " parfile_values = {}\n", " for line in parfile_path.read_text().splitlines():\n", " if \"=\" not in line or line.lstrip().startswith(\"#\"):\n", " continue\n", " key, value = line.split(\"=\", 1)\n", " key = key.strip()\n", " if key in selected_keys:\n", " parfile_values[key] = float(value)\n", "\n", "expected_parfile_values = {\n", " \"par_b\": adm_id[\"D\"] / 2.0,\n", " \"par_P_plus1\": adm_id[\"px\"],\n", " \"par_P_plus2\": adm_id[\"py\"],\n", " \"par_P_minus1\": -adm_id[\"px\"],\n", " \"par_P_minus2\": -adm_id[\"py\"],\n", " \"center_offset1\": adm_id[\"x_offset\"],\n", "}\n", "\n", "print(f\"parfile written to: {parfile_path}\")\n", "print()\n", "for key in sorted(selected_keys):\n", " actual = parfile_values[key]\n", " expected = expected_parfile_values[key]\n", " print(f\"{key:14s} actual = {actual: .12e} expected = {expected: .12e}\")\n", " assert np.isclose(actual, expected)" ] }, { "cell_type": "markdown", "id": "fd27346c", "metadata": {}, "source": [ "## Summary\n", "\n", "In this notebook we have shown how to convert EOB initial data into ADM puncture data, and then into a TwoPunctures initial-data dictionary.\n", "We also showed how to write a TwoPunctures parfile from the same ADM data, and verified the key consistency relations between the geometric center offset, the puncture positions, and the parfile entries.\n", "\n", "Note: the parfile helper currently follows the capabilities of `TwoPunctID`, so it assumes that the\n", "punctures start on the x-axis and with no z-component of the momenta (i.e. `\"center_offset2\": 0.0`, `\"center_offset3\": 0.0` and `\"par_P_plus3\": 0.0`, `\"par_P_minus3\": 0.0`)\n", "in the generated parfile." ] }, { "cell_type": "markdown", "id": "8a73410f", "metadata": {}, "source": [] } ], "metadata": { "kernelspec": { "display_name": "pyart-dev-3.11", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.11" } }, "nbformat": 4, "nbformat_minor": 5 }