{ "cells": [ { "cell_type": "markdown", "id": "e38aaf25", "metadata": {}, "source": [ "# Grid PSR J0740 over distance, companion mass, cos i with PBDOT constraint and derived parameters" ] }, { "cell_type": "markdown", "id": "93129206", "metadata": {}, "source": [ "Try to reproduce Figure 6 from [Fonseca et al. 2021, ApJ, 915, L12](https://ui.adsabs.harvard.edu/abs/2021ApJ...915L..12F/abstract)\n", "\n", "We have downloaded wideband files from [Zenodo](https://zenodo.org/record/4773599#.YWdRDC-cbUI):\n", "* [J0740+6620.FCP+21.wb.tim](https://zenodo.org/record/4773599/files/J0740%2B6620.FCP%2B21.wb.tim?download=1)\n", "* [J0740+6620.FCP+21.wb.DMX3.0.par](https://zenodo.org/record/4773599/files/J0740%2B6620.FCP%2B21.wb.DMX3.0.par?download=1)" ] }, { "cell_type": "code", "execution_count": null, "id": "ec5d44fc", "metadata": {}, "outputs": [], "source": [ "# imports\n", "from astropy import units as u, constants as c\n", "import numpy as np\n", "\n", "import pint.config\n", "import pint.gridutils\n", "import pint.models.parameter as param\n", "import pint.residuals\n", "import pint.fitter\n", "import pint.derived_quantities\n", "from pint.models.model_builder import get_model, get_model_and_toas\n", "from pint.toa import get_TOAs\n", "\n", "import concurrent.futures\n", "import datetime\n", "import multiprocessing" ] }, { "cell_type": "markdown", "id": "e5e226a8", "metadata": {}, "source": [ "We want to set up a grid to cover:\n", "* distance $d$\n", "* $\\cos i$\n", "* companion mass $M_c$\n", "\n", "However, the intrinsic parameters are:\n", "* parallax $\\varpi$\n", "* $\\sin i$\n", "* companion mass $M2$\n", "* orbital decay rate $\\dot P_B$\n", "\n", "These depend on the grid parameters in various ways, some easy, some complicated. So:\n", "* $d=1/\\varpi$\n", "* $M_c = M2$\n", "* $\\sin i = \\sqrt{1-\\cos^2 i}$\n", "\n", "What about $\\dot P_B$? That is harder. We need to compute the expected $\\dot P_B$ from GR, and compare it with the measured value corrected for kinematic effects (Shklovskii effect, Galactic acceleration). And we also know that $\\dot P_B(GR)$ depends on both $M_c$ and pulsar mass $M_p$. But we can connect $M_c$ and $\\cos i$ to get $M_p$, and then $M_p$ and $M_c$ to get $\\dot P_B(GR)$. We compute the various kinematic effects using the same analytic versions in Fonseca et al.\n", "\n", "We will make an object to bundle all of these computations together to run in each iteration." ] }, { "cell_type": "code", "execution_count": null, "id": "8c7356fb", "metadata": {}, "outputs": [], "source": [ "class constrained_gridder:\n", " def __init__(\n", " self,\n", " fitter,\n", " R0=8.178 * u.kpc,\n", " Theta0=236.9 * u.km / u.s,\n", " rho0=1e-2 * u.Msun / u.pc ** 3,\n", " Sigma0=46 * u.Msun / u.pc ** 2,\n", " z0=0.18 * u.kpc,\n", " ):\n", " \"\"\"Compute the chi^2 given various input parameters that are then compared with the model\n", " \n", " Parameters\n", " ----------\n", " fitter : pint.fitter.Fitter\n", " R0 : astropy.units.Quantity, optional\n", " Galactic radius of the Sun\n", " Theta0 : astropy.units.Quantity, optional\n", " Galactic rotation speed for the Sun \n", " rho0 : astropy.units.Quantity, optional\n", " local mass density of the Galactic disk\n", " Sigma0 : astropy.units.Quantity, optional\n", " local mass surface density of the Galactic disk\n", " z0 : astropy.units.Quantity, optional\n", " local scale-height of the Galactic disk\n", " \"\"\"\n", " self.fitter = fitter\n", " self.e = np.sqrt(\n", " fitter.model.EPS1.quantity ** 2 + fitter.model.EPS2.quantity ** 2\n", " )\n", " self.coords_galactic = fitter.model.coords_as_GAL()\n", " self.mu = np.sqrt(\n", " self.coords_galactic.pm_l_cosb ** 2 + self.coords_galactic.pm_b ** 2\n", " )\n", " self.R0 = R0\n", " self.Theta0 = Theta0\n", " self.rho0 = rho0\n", " self.Sigma0 = Sigma0\n", " self.z0 = z0\n", "\n", " def Mp(self, Mc, cosi, d):\n", " \"\"\"Pulsar mass as a function of companion mass and inclination\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance (not used, just for consistency among function calls) \n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " i = np.arccos(cosi) * u.rad\n", " return pint.derived_quantities.pulsar_mass(\n", " self.fitter.model.PB.quantity, self.fitter.model.A1.quantity, Mc, i\n", " )\n", "\n", " def PBDOT_GR(self, Mc, cosi, d):\n", " \"\"\"Orbital decay due to GR, as a function of companion mass and inclination\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance (not used, just for consistency among function calls)\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " Mp = self.Mp(Mc, cosi, d)\n", " return pint.derived_quantities.pbdot(\n", " Mp, Mc, self.fitter.model.PB.quantity, self.e\n", " ).to(u.d / u.s)\n", "\n", " def PBDOT_mu(self, Mc, cosi, d):\n", " \"\"\"Apparent orbital decay due to proper motion (Shklovskii)\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \n", " Notes\n", " -----\n", " Proper motion is stored in the ``fitter`` object\n", " \"\"\"\n", " return ((self.mu ** 2 * d / c.c) * self.fitter.model.PB.quantity).to(\n", " u.d / u.s, equivalencies=u.dimensionless_angles()\n", " )\n", "\n", " def PBDOT_DR(self, Mc, cosi, d):\n", " \"\"\"Apparent orbital decay due to Galactic acceleration in the plane\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " # this is beta in Nice & Taylor (1995), after Eqn. 5\n", " kappa = (d / self.R0) * np.cos(self.coords_galactic.b) - np.cos(\n", " self.coords_galactic.l\n", " )\n", " # and Nice & Talor (1995), Eqn. 5\n", " return (\n", " -self.fitter.model.PB.quantity\n", " * (self.Theta0 ** 2 / c.c / self.R0)\n", " * (\n", " np.cos(self.coords_galactic.l)\n", " + kappa ** 2 / (np.sin(self.coords_galactic.l) ** 2 + kappa ** 2)\n", " )\n", " ).to(u.d / u.s)\n", "\n", " def PBDOT_z(self, Mc, cosi, d):\n", " \"\"\"Apparent orbital decay due to Galactic acceleration in the vertical direction\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " z = d * np.sin(self.coords_galactic.b)\n", " # This was:\n", " # Nice & Taylor (1995, ApJ, 441, 429), Eqn. 4\n", " # which used values from:\n", " # Kuijken and Gilmore (1989, MNRAS, 239, 605)\n", " # we will recast it in terms of the original measured quantities\n", " # so that they can be re-evaluated if needed\n", " # (they had be written in funny units)\n", " # it differs from the Nice & Taylor value by a few percent: maybe some\n", " # difference in constant conversions\n", " az = (\n", " -2\n", " * np.pi\n", " * c.G\n", " * (self.Sigma0 * z / np.sqrt(z ** 2 + self.z0 ** 2) + 2 * self.rho0 * z)\n", " )\n", " return (\n", " (az * self.fitter.model.PB.quantity / c.c) * np.sin(self.coords_galactic.b)\n", " ).to(u.s / u.s)\n", "\n", " def PBDOT(self, Mc, cosi, d):\n", " \"\"\"Total orbital decay due to GR, Shklovskii, Galactic acceleration\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass\n", " cosi : float\n", " cosine of inclination\n", " d : astropy.units.Quantity\n", " distance\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " return (\n", " self.PBDOT_GR(Mc, cosi, d)\n", " + self.PBDOT_DR(Mc, cosi, d)\n", " + self.PBDOT_mu(Mc, cosi, d)\n", " + self.PBDOT_z(Mc, cosi, d)\n", " )\n", "\n", " def parallax(self, Mc, cosi, d):\n", " \"\"\"parallax from distance\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass (not used, just for consistency among function calls)\n", " cosi : float\n", " cosine of inclination (not used, just for consistency among function calls)\n", " d : astropy.units.Quantity\n", " distance\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " return d.to(u.mas, equivalencies=u.parallax())\n", "\n", " def sini(self, Mc, cosi, d):\n", " \"\"\"sine of inclination from cos of inclination\n", " \n", " Parameters\n", " ----------\n", " Mc : astropy.units.Quantity\n", " Companion mass (not used, just for consistency among function calls)\n", " cosi : float\n", " cosine of inclination \n", " d : astropy.units.Quantity\n", " distance (not used, just for consistency among function calls)\n", " \n", " Returns\n", " -------\n", " astropy.units.Quantity\n", " \"\"\"\n", " return np.sqrt(1 - cosi ** 2)" ] }, { "cell_type": "code", "execution_count": null, "id": "c138a116", "metadata": {}, "outputs": [], "source": [ "print(f\"Starting at {datetime.datetime.now().isoformat()}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "941f12ee", "metadata": {}, "outputs": [], "source": [ "# Load in a basic dataset, construct a fitter, do an initial fit\n", "parfile = pint.config.examplefile(\"J0740+6620.FCP+21.wb.DMX3.0.par\")\n", "timfile = pint.config.examplefile(\"J0740+6620.FCP+21.wb.tim\")\n", "m, t = get_model_and_toas(parfile, timfile)\n", "\n", "f = pint.fitter.WidebandDownhillFitter(t, m)\n", "\n", "f.fit_toas()\n", "print(\n", " f\"Computed best-fit chi2={f.resids.chi2} at {datetime.datetime.now().isoformat()}\"\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "09aebb61", "metadata": {}, "outputs": [], "source": [ "# make the constrained gridder object out of the fitter\n", "cg = constrained_gridder(f)" ] }, { "cell_type": "code", "execution_count": null, "id": "4d79077e", "metadata": {}, "outputs": [], "source": [ "# set up the grid. The number of grid points will depend on your available CPU/time\n", "distance = np.linspace(0.5, 2, 1) * u.kpc\n", "cosi = np.linspace(0.03, 0.06, 1)\n", "Mc = np.linspace(0.23, 0.28, 1) * u.Msun\n", "outfile = \"J0740_PINT_grid.npz\"" ] }, { "cell_type": "code", "execution_count": null, "id": "e19828c5", "metadata": {}, "outputs": [], "source": [ "# may be needed for some executables\n", "# multiprocessing.freeze_support()" ] }, { "cell_type": "code", "execution_count": null, "id": "2ef0a863", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "# use this for SLURM\n", "# with cfut.SlurmExecutor(debug=True,keep_logs=False,additional_setup_lines=slurm_options) as executor:\n", "\n", "# use this for MPI\n", "# with MPIPoolExecutor() as executor:\n", "\n", "# default (could also just not supply any executor)\n", "with concurrent.futures.ProcessPoolExecutor() as executor:\n", "\n", " # the input grid is (Mc, cosi, distance)\n", " # that gets transformed to the intrinsic measured quantities (M2, SINI, PX, PBDOT)\n", " # via\n", " # M2 = Mc\n", " # sini = cg.sini(Mc, cosi, d)\n", " # parallax = cg.parallax(Mc, cosi, d)\n", " # PBDOT = cg.PBDOT(Mc, cosi, d)\n", " # note that all of those have the same input parameters,\n", " # even if they aren't needed for a specific calculation\n", " chi2, params = pint.gridutils.grid_chisq_derived(\n", " f,\n", " (\"M2\", \"SINI\", \"PX\", \"PBDOT\"),\n", " (lambda Mc, cosi, d: Mc, cg.sini, cg.parallax, cg.PBDOT),\n", " (Mc, cosi, distance),\n", " executor=executor,\n", " printprogress=False,\n", " )\n", " # save this if desired\n", " # np.savez(\n", " # outfile,\n", " # distance=distance,\n", " # cosi=cosi,\n", " # Mc=Mc,\n", " # chi2=chi2,\n", " # params=params,\n", " # )\n", " # print(f\"Wrote {outfile}\")\n", " print(\n", " f\"Finished {str(chi2.shape)}={len(chi2.flatten())} points at {datetime.datetime.now().isoformat()}\"\n", " )" ] }, { "cell_type": "code", "execution_count": null, "id": "1e421a7c", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "formats": "ipynb,py:percent" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }