{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1b2c3d4",
   "metadata": {},
   "source": [
    "# Working with ICCUB Catalog and Waveform Integration\n",
    "\n",
    "This tutorial demonstrates how to work with waveforms from the ICCUB (Institute of Cosmos Sciences) numerical relativity catalog and perform waveform integration from Ψ₄ (Weyl scalar) to strain h.\n",
    "\n",
    "ICCUB simulations often provide Ψ₄ data, which needs to be integrated twice to obtain the gravitational wave strain. PyART provides tools for this integration.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5f6g7h8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from PyART.catalogs.icc import Waveform_ICC\n",
    "from PyART.analysis.integrate_multipole import Multipole\n",
    "from PyART.analysis.integrate_wave import IntegrateMultipole"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "i9j0k1l2",
   "metadata": {},
   "source": [
    "## Load ICCUB Waveform\n",
    "\n",
    "First, we load a waveform from the ICCUB catalog. We'll load it in two ways:\n",
    "1. Without integration (raw Ψ₄ data)\n",
    "2. With integration (computed strain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "m3n4o5p6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to ICCUB catalog data\n",
    "# Note: Adjust this path to your local data directory\n",
    "path = './local_data/icc/catalog'\n",
    "ID = 1  # Simulation ID\n",
    "\n",
    "# Integration options\n",
    "integr_opts = {\n",
    "    'f0': 0.002,         # Low-frequency cutoff\n",
    "    'extrap_psi4': True, # Extrapolate Psi4 to infinity\n",
    "    'method': 'FFI',     # Fixed Frequency Integration method\n",
    "    'window': [20, -20]  # Time window for integration\n",
    "}\n",
    "\n",
    "# Load without integration (raw data)\n",
    "print(\"Loading raw Psi4 data...\")\n",
    "icc = Waveform_ICC(path=path, ID=ID)\n",
    "\n",
    "# Load with integration\n",
    "print(\"Loading and integrating to strain...\")\n",
    "iicc = Waveform_ICC(\n",
    "    path=path, \n",
    "    ID=ID, \n",
    "    integrate=True,\n",
    "    integr_opts=integr_opts,\n",
    "    load_puncts=True  # Also load puncture trajectory data\n",
    ")\n",
    "\n",
    "print(\"\\nWaveform loaded successfully!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "q7r8s9t0",
   "metadata": {},
   "source": [
    "## Display Metadata\n",
    "\n",
    "Let's examine the simulation parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "u1v2w3x4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Simulation Metadata:\")\n",
    "print('=' * 60)\n",
    "for k, v in icc.metadata.items():\n",
    "    print(f'{k:20s} : {v}')\n",
    "print('=' * 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "y5z6a7b8",
   "metadata": {},
   "source": [
    "## Compare Stored vs Integrated Waveforms\n",
    "\n",
    "If the simulation provides both Ψ₄ and pre-computed strain, we can compare with our integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9d0e1f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find merger times\n",
    "try:\n",
    "    iicc_tmrg, _, _, _ = iicc.find_max(height=0.3)\n",
    "    icc_tmrg, _, _, _ = icc.find_max(height=0.3)\n",
    "except:\n",
    "    iicc_tmrg = 0\n",
    "    icc_tmrg = 0\n",
    "\n",
    "# Define waveform types to plot\n",
    "waves = ['hlm', 'dothlm', 'psi4lm']\n",
    "ylabs = [r'$h_{22}$', r'$\\dot{h}_{22}$', r'$\\psi_4^{22}$']\n",
    "\n",
    "npanels = len(waves)\n",
    "plt.figure(figsize=(12, 3*npanels))\n",
    "plt.suptitle(f'ICCUB Simulation ID: {ID:04d}', fontsize=18, y=0.995)\n",
    "\n",
    "for i, wave_name in enumerate(waves):\n",
    "    plt.subplot(npanels, 1, i+1)\n",
    "    \n",
    "    # Get integrated waveform time\n",
    "    if wave_name == 'psi4lm':\n",
    "        iicc_t = iicc.t_psi4\n",
    "    else:\n",
    "        iicc_t = iicc.u - iicc_tmrg\n",
    "    \n",
    "    # Plot stored waveform if available\n",
    "    icc_wave = getattr(icc, wave_name)\n",
    "    if icc_wave:\n",
    "        if wave_name == 'psi4lm':\n",
    "            icc_t = icc.t_psi4\n",
    "        else:\n",
    "            icc_t = icc.u - icc_tmrg\n",
    "        \n",
    "        plt.plot(icc_t, icc_wave[(2,2)]['A'], \n",
    "                lw=1.5, c='red', alpha=0.7, label='Stored')\n",
    "        plt.plot(icc_t, icc_wave[(2,2)]['real'], \n",
    "                lw=1.0, c='red', alpha=0.5)\n",
    "        plt.plot(icc_t, icc_wave[(2,2)]['imag'], \n",
    "                lw=1.0, c='red', ls=':', alpha=0.5)\n",
    "    \n",
    "    # Plot integrated waveform\n",
    "    iicc_wave = getattr(iicc, wave_name)\n",
    "    plt.plot(iicc_t, iicc_wave[(2,2)]['A'], \n",
    "            lw=1.5, c='blue', alpha=0.7, label='Integrated')\n",
    "    plt.plot(iicc_t, iicc_wave[(2,2)]['real'], \n",
    "            lw=1.0, c='blue', alpha=0.5)\n",
    "    plt.plot(iicc_t, iicc_wave[(2,2)]['imag'], \n",
    "            lw=1.0, c='blue', ls=':', alpha=0.5)\n",
    "    \n",
    "    plt.ylabel(ylabs[i], fontsize=14)\n",
    "    plt.xlabel(r'$t-t_{\\rm mrg}$ (M)', fontsize=14) if i == len(waves)-1 else None\n",
    "    plt.legend(fontsize=11)\n",
    "    plt.grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "g3h4i5j6",
   "metadata": {},
   "source": [
    "## Visualize Puncture Trajectories\n",
    "\n",
    "For scattering simulations, we can plot the trajectories of the two black holes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "k7l8m9n0",
   "metadata": {},
   "outputs": [],
   "source": [
    "if iicc.metadata['scat_angle'] is not None:\n",
    "    plt.figure(figsize=(10, 8))\n",
    "    \n",
    "    # Plot both puncture trajectories\n",
    "    plt.plot(iicc.puncts['x0'], iicc.puncts['y0'], \n",
    "            linewidth=2, label='BH 1', marker='o', markersize=3)\n",
    "    plt.plot(iicc.puncts['x1'], iicc.puncts['y1'], \n",
    "            linewidth=2, label='BH 2', marker='s', markersize=3)\n",
    "    \n",
    "    plt.xlabel(r'$x$ (M)', fontsize=16)\n",
    "    plt.ylabel(r'$y$ (M)', fontsize=16)\n",
    "    plt.title(f'Black Hole Trajectories (Scattering Angle: {iicc.metadata[\"scat_angle\"]:.2f}°)', \n",
    "             fontsize=16)\n",
    "    plt.legend(fontsize=14)\n",
    "    plt.grid(True, alpha=0.3)\n",
    "    plt.axis('equal')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "else:\n",
    "    print(\"This is not a scattering simulation - no trajectory plot available.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "o1p2q3r4",
   "metadata": {},
   "source": [
    "## Testing the Integration Method\n",
    "\n",
    "We can explicitly test the integration by comparing old and new integration methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "s5t6u7v8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Psi4 data\n",
    "t = icc.t_psi4\n",
    "psi4 = icc.psi4lm[(2,2)]['z']\n",
    "l, m = 2, 2\n",
    "M = 1\n",
    "\n",
    "# New integration method\n",
    "mode_new = IntegrateMultipole(\n",
    "    l, m, t, psi4, \n",
    "    **integr_opts,\n",
    "    mass=M, \n",
    "    radius=icc.r_extr, \n",
    "    integrand='psi4'\n",
    ")\n",
    "\n",
    "# Old integration method\n",
    "mode_old = Multipole(l, m, t, psi4, mass=M, radius=icc.r_extr, integrand='psi4')\n",
    "mode_old.hlm, mode_old.dothlm = mode_old.integrate_wave(integr_opts=integr_opts)\n",
    "\n",
    "# Plot comparison\n",
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "# Strain h\n",
    "plt.subplot(3, 2, 1)\n",
    "plt.plot(t, mode_old.hlm.real, c='r', ls='-', label='Old method', linewidth=2)\n",
    "plt.plot(t, mode_new.h.real, c='b', ls='--', label='New method', linewidth=2, alpha=0.7)\n",
    "plt.ylabel(r'Re[$h_{22}$]', fontsize=14)\n",
    "plt.legend(fontsize=12)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Strain Comparison', fontsize=14)\n",
    "\n",
    "plt.subplot(3, 2, 2)\n",
    "plt.plot(t, np.abs(mode_old.hlm - mode_new.h), c='green', label='Difference')\n",
    "plt.ylabel(r'$|h_{\\rm old} - h_{\\rm new}|$', fontsize=14)\n",
    "plt.yscale('log')\n",
    "plt.legend(fontsize=12)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Absolute Difference', fontsize=14)\n",
    "\n",
    "# Strain derivative\n",
    "plt.subplot(3, 2, 3)\n",
    "plt.plot(t, mode_old.dothlm.real, c='r', ls='-', linewidth=2)\n",
    "plt.plot(t, mode_new.dh.real, c='b', ls='--', linewidth=2, alpha=0.7)\n",
    "plt.ylabel(r'Re[$\\dot{h}_{22}$]', fontsize=14)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Strain Derivative', fontsize=14)\n",
    "\n",
    "plt.subplot(3, 2, 4)\n",
    "plt.plot(t, np.abs(mode_old.dothlm - mode_new.dh), c='green')\n",
    "plt.ylabel(r'$|\\dot{h}_{\\rm old} - \\dot{h}_{\\rm new}|$', fontsize=14)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Derivative Difference', fontsize=14)\n",
    "\n",
    "# Psi4\n",
    "plt.subplot(3, 2, 5)\n",
    "plt.plot(t, mode_old.psi.real, c='r', ls='-', linewidth=2)\n",
    "plt.plot(t, mode_new.psi4.real, c='b', ls='--', linewidth=2, alpha=0.7)\n",
    "plt.xlabel('Time (M)', fontsize=14)\n",
    "plt.ylabel(r'Re[$\\psi_4^{22}$]', fontsize=14)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Weyl Scalar', fontsize=14)\n",
    "\n",
    "plt.subplot(3, 2, 6)\n",
    "plt.plot(t, np.abs(mode_old.psi - mode_new.psi4), c='green')\n",
    "plt.xlabel('Time (M)', fontsize=14)\n",
    "plt.ylabel(r'$|\\psi_4^{\\rm old} - \\psi_4^{\\rm new}|$', fontsize=14)\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.title('Psi4 Difference', fontsize=14)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"\\nMaximum relative difference in h: {np.max(np.abs(mode_old.hlm - mode_new.h)/np.max(np.abs(mode_old.hlm))):.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "w9x0y1z2",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This tutorial demonstrated:\n",
    "- Loading waveforms from the ICCUB catalog\n",
    "- Integrating Ψ₄ (Weyl scalar) to obtain gravitational wave strain\n",
    "- Comparing stored and integrated waveforms\n",
    "- Visualizing puncture trajectories for scattering simulations\n",
    "- Testing and validating integration methods\n",
    "\n",
    "### Key Points\n",
    "\n",
    "- ICCUB provides high-quality numerical relativity data including scattering scenarios\n",
    "- The Fixed Frequency Integration (FFI) method provides accurate strain from Ψ₄\n",
    "- Integration requires careful choices of f₀ (low-frequency cutoff) and time windows\n",
    "- PyART's integration tools handle extrapolation and windowing automatically\n",
    "\n",
    "### Next Steps\n",
    "\n",
    "- Explore different integration methods\n",
    "- Compute scattering angles from trajectory data\n",
    "- Compare ICCUB waveforms with other catalogs\n",
    "- Perform mismatch calculations with EOB models"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}