{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to get RMs from FITS images" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given appropriate metadata, it is possible to obtain ionospheric RMs from FITS images. Further, [Van Eck (2021)](https://github.com/CIRADA-Tools/FRion/blob/main/docs/source/Ionospheric_Correction.pdf) derived a method to obtain a time-integrated correction suitable for applying to image data. There is also Python implementation [`FRion`](https://github.com/CIRADA-Tools/FRion), which we have partly re-implemented here in `spinifex` for tight integration with the underlying ionospheric modelling routines.\n", "\n", "\n", "Van Eck (2021) Equation 3 provides the time-integrated effect on linear polarisation as:\n", "$$\n", "\\tilde\\Theta(\\lambda^2) = \\frac{1}{t_1 - t_0} \\int_{t_0}^{t_1} e^{2i\\lambda^2\\phi(t)}dt.\n", "$$\n", "Where $t_0$ and $t_1$ are the start and end times, respecitively, $\\lambda^2$ are the observed wavelengths-squared, and $\\phi$ is the time-dependent Faraday rotation. $\\tilde\\Theta$ is a complex quantity whose amplitude can be interpreted as the depolarisation caused by the ionosphere, and the phase as the change of polarisation angle." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this example, we'll generate a mock FITS file using `astropy`. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from __future__ import annotations\n", "\n", "from pathlib import Path\n", "from pprint import pprint\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from astropy import units as u\n", "from astropy.io import fits\n", "from astropy.visualization import quantity_support, time_support\n", "from spinifex.image_tools import (\n", " get_freq_from_fits,\n", " get_integrated_rm_from_fits,\n", " get_rm_from_fits,\n", ")\n", "\n", "_ = quantity_support()\n", "_ = time_support()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "example_fits_file = Path(\"example.fits\")\n", "\n", "header = fits.Header()\n", "header[\"NAXIS\"] = 4\n", "header[\"NAXIS1\"] = 10\n", "header[\"NAXIS2\"] = 10\n", "header[\"NAXIS3\"] = 6\n", "header[\"NAXIS4\"] = 1\n", "header[\"CTYPE1\"] = \"RA---SIN\"\n", "header[\"CRVAL1\"] = 0\n", "header[\"CRPIX1\"] = 5\n", "header[\"CDELT1\"] = -(1 / 3600)\n", "header[\"CUNIT1\"] = \"deg\"\n", "header[\"CTYPE2\"] = \"DEC--SIN\"\n", "header[\"CRVAL2\"] = 0\n", "header[\"CRPIX2\"] = 5\n", "header[\"CDELT2\"] = 1 / 3600\n", "header[\"CUNIT2\"] = \"deg\"\n", "header[\"CTYPE3\"] = \"FREQ\"\n", "header[\"CRVAL3\"] = 1.4e9\n", "header[\"CRPIX3\"] = 1\n", "header[\"CDELT3\"] = 1e6\n", "header[\"CUNIT3\"] = \"Hz\"\n", "header[\"CTYPE4\"] = \"STOKES\"\n", "header[\"CRVAL4\"] = 1\n", "header[\"CRPIX4\"] = 1\n", "header[\"CDELT4\"] = 1\n", "header[\"CUNIT4\"] = \"\"\n", "header[\"BUNIT\"] = \"Jy/beam\"\n", "header[\"DATE-OBS\"] = (\n", " \"2019-04-25T12:45:52.893302\" # pulled from a random SPICE-RACS file\n", ")\n", "header[\"DURATION\"] = 900\n", "header[\"TELESCOP\"] = \"ASKAP\"\n", "data = np.ones((1, 6, 10, 10), dtype=np.float32)\n", "hdu = fits.PrimaryHDU(data, header)\n", "\n", "hdu.writeto(example_fits_file, overwrite=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get the time-dependent RM - this returns the same object as the other `get_rm` functions\n", "rm = get_rm_from_fits(\n", " fits_path=example_fits_file,\n", " timestep=1\n", " * u.min, # timestep is optional, defaults to 15 minutes - but that's our total duration here!\n", ")\n", "\n", "\n", "fig, ax = plt.subplots()\n", "ax.errorbar(\n", " rm.times.datetime,\n", " rm.rm * u.rad / u.m**2,\n", " yerr=rm.rm_error * u.rad / u.m**2,\n", " fmt=\".\",\n", ")\n", "ax.set(\n", " xlabel=\"Time\",\n", " ylabel=f\"RM / {u.rad / u.m**2:latex_inline}\",\n", " title=\"Time-dependent RM\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get the time-integrated product\n", "integrated_rm = get_integrated_rm_from_fits(\n", " fits_path=example_fits_file,\n", " timestep=1\n", " * u.min, # timestep is optional, defaults to 15 minutes - but that's our total duration here!\n", ")\n", "# We'll use our helper function to get the frequency axis\n", "frequencies = get_freq_from_fits(example_fits_file)\n", "fig, ax = plt.subplots()\n", "ax2 = ax.twinx()\n", "ax.plot(frequencies.to(u.GHz), np.abs(integrated_rm.theta), label=r\"$|\\Theta|$\")\n", "ax.set(\n", " ylim=(0, 1.1),\n", " ylabel=\"Fractional polarisation\",\n", " xlabel=f\"Frequency / {u.GHz:latex_inline}\",\n", " title=\"Time-integrated polarisation effect\",\n", ")\n", "ax2.plot(\n", " frequencies.to(u.GHz),\n", " (np.angle(integrated_rm.theta) * u.rad).to(u.deg),\n", " color=\"tab:orange\",\n", " label=r\"$\\arg(\\Theta)$\",\n", ")\n", "ax2.set(ylim=(-90, 90), ylabel=\"Polarisation angle\")\n", "fig.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# The integrated RM object is a namedtuple with the following fields\n", "# These are mostly the same as the RM object, but averaged over time\n", "pprint(integrated_rm._asdict())" ] } ], "metadata": { "kernelspec": { "display_name": "spinifex", "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.12.9" } }, "nbformat": 4, "nbformat_minor": 2 }