ToOSimulations_with_redback_arxiv_plots_updated.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook used with redback to generate ToO style plots for the publication version of the ToO report\n",
    "\n",
    "### You will need to install Redback. Instructions available at https://redback.readthedocs.io/en/latest/. I suggest installing from source via GitHub."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook Template Author: Meet Vyas <mjvyas2002@gmail.com>\n",
    "\n",
    "Code reviewed by Igor Andreoni and Ritwik Sharma "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Instructions to create the plots using Google Colab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### 1. Clone the repository in your colab environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!git clone https://github.com/nikhil-sarin/redback"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### 2. Make the cloned repository your current directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd redback/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### 3. Install the requirements for redback "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### 4. Also install the additional requirements "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install lalsuite nestle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### 5. Redback uses latex for character typesetting and therefore download latex for google colab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sudo apt-get install -y texlive texlive-latex-extra texlive-fonts-recommended dvipng cm-super"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This ends the google colab instructions, follow along on your local machine from here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Import all the required libraries and dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import redback\n",
    "import pandas as pd\n",
    "from redback.simulate_transients import SimulateOpticalTransient\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We use the Golden Strategy which comprises of a 3 filter deep strategy. This takes the form of a dataframe which specifies time from transient onset, right ascension, declination, filters and five sigma depth mags. The three cases for the Golden Strategy are given as follows \n",
    "\n",
    "> Night 0 timing considerations: we calculate that it will take 95 minutes (1.59 hrs) to scan a 100 sq deg skymap and the scheduling of the three epochs on Night 0 will depend on the observability window of the sky map (Twin ). Here we provide a timing strategy:\n",
    "> > ➢ Case 1: If Twin >=3*1.59 = 4.76 hrs; do all 3 scans, with the gap between them of ((Twin - 3*1.59)/2)\n",
    "> > \n",
    "> > ➢ Case 2: Else if 3.2 < Twin < 4.7 hrs; do 2 scans. The start time of the visits is separated by (Twin - 2*1.59), which is sufficient for asteroid rejection and short timescale evolution.\n",
    "> > \n",
    "> > ➢ Case 3: Else if Twin < 3.2 hrs then do 1x1.59 hrs scan + followed by 1 scan (0.53 hrs) in r-band only after the full filter complement is complete, to ensure asteroid rejection is reliable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the pointings file we have used the 5 sigma magnitudes calculated from etc.py based on the criteria for 120 seconds of exposure time on day 0 and 180 seconds on Day 1,2, and 3 as given in the Report, attached below for reference \n",
    "\n",
    " > ➢ Night 0: three scans across the skymap in 3 filters with 120-sec exposure for each filter\n",
    "gri8\n",
    "Each scan is estimated to take ~95 minutes.\n",
    " > \n",
    " > ➢ Night 1,2,3: one scan across the skymap in two filters (likely r+i) with 180 seconds per\n",
    "visit. Each night is estimated to take ~92 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### The manual sky pointings csv files need to be read using the pandas library "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the csv file for case 1 of the Golden Strategy\n",
    "pointings_case_1 = pd.read_csv(\"Golden_Strategy_Case_1.csv\")\n",
    "print(pointings_case_1)\n",
    "\n",
    "# Load the csv file for case 2 of the Golden Strategy\n",
    "pointings_case_2 = pd.read_csv(\"Golden_Strategy_Case_2.csv\")\n",
    "print(pointings_case_2)\n",
    "\n",
    "# Load the csv file for case 3 of the Golden Strategy\n",
    "pointings_case_3 = pd.read_csv(\"Golden_Strategy_Case_3.csv\")\n",
    "print(pointings_case_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### Below we will load the parameters for the kilonova model. The model takes redshift as a parameter therefore we will convert the luminosity distance to redshift which we do using the planck18 function using astropy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The redshift corresponding to a distance of 350.0 Mpc is: 0.07482 redshift\n"
     ]
    }
   ],
   "source": [
    "from astropy.cosmology import Planck18 as cosmo\n",
    "import astropy.units as u\n",
    "from astropy.cosmology import z_at_value\n",
    "\n",
    "# Ask the user for the distance in Mpc\n",
    "distance_mpc = float(input(\"Enter the distance in Mpc: \")) * u.Mpc\n",
    "\n",
    "# Calculate the redshift corresponding to the user-provided distance\n",
    "redshift = z_at_value(cosmo.luminosity_distance, distance_mpc)\n",
    "print(f\"The redshift corresponding to a distance of {distance_mpc} is: {redshift:.5f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### We are using the One Component kilonova model for the lightcurve simulations for a Faint kilonova from the redback library .\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> First we update the parameters for the model at a luminosity distance of 350 Mpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set model-specific keyword arguments.\n",
    "model_kwargs = {'output_format': 'magnitude'}\n",
    "\n",
    "# Specify a model. Here, any redback model can be referred to as a string. \n",
    "# If the user has their own model, they can pass a function instead. \n",
    "# There are over 100 models in redback, including options for kilonovae, GRB afterglows,\n",
    "# supernovae, TDEs, and other transient events.\n",
    "model = 'one_component_kilonova_model'\n",
    "\n",
    "# Load and sample from the default prior for the chosen model, generating parameters for a random kilonova.\n",
    "# This creates a base set of parameters that can be adjusted as needed.\n",
    "parameters = redback.priors.get_priors(model=model).sample()\n",
    "\n",
    "# We can also set the parameters manually. For the arxiv plots we use these parameters for the 350 Mpc kilonova event.\n",
    "parameters['mej'] = 0.05                 # Ejecta mass in solar masses\n",
    "parameters['t0_mjd_transient'] = 0       # Start time of the transient in MJD (Modified Julian Date)\n",
    "parameters['redshift'] = 0.074816848     # Redshift for the kilonova event at 350 Mpc \n",
    "parameters['t0'] = parameters['t0_mjd_transient']  # Set the transient's starting time\n",
    "parameters['temperature_floor'] = 3000   # Minimum temperature in Kelvin (floor value) to avoid extremely low temperatures in the model\n",
    "parameters['kappa'] = 1                  # Gray opacity, representing the opacity of the ejecta material\n",
    "parameters['vej'] = 0.2                  # Minimum initial velocity of the ejecta in units of c (speed of light)\n",
    "parameters['ra'] = 1.0                   # Right ascension in degrees for the kilonova's location\n",
    "parameters['dec'] = 1.5                  # Declination in degrees for the kilonova's location\n",
    "print(parameters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One can check the validity of the parameters used based on the priors defined in the library, the priors for the one component kilonova can be loaded as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "priors = redback.priors.get_priors(model=model)\n",
    "priors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Parameters for the model at a luminosity distance of 700 Mpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_kwargs = {'output_format': 'magnitude'}\n",
    "# Any redback model can be referred to as a string. \n",
    "# If the user has their own model, they can pass a function here instead. \n",
    "# There are over a 100 models implemented in redback, lots of models for kilonovae, GRB afterglows, \n",
    "# supernovae, TDEs and other things\n",
    "model = 'one_component_kilonova_model'\n",
    "# Load the default prior for this model in redback and sample from it to get 1 set of parameters. \n",
    "# We can sample from the default prior for this model for a random kilonova. \n",
    "parameters = redback.priors.get_priors(model=model).sample()\n",
    "\n",
    "# We fix a few parameters here to create a nice looking kilonova. \n",
    "# You can change any of the parameters here or add additional keyword arguments \n",
    "# to change some physical assumptions. Please refer to the documentation for this and units etc\n",
    "parameters['mej'] = 0.05 # Ejecta mass in solar masses\n",
    "parameters['t0_mjd_transient'] = 0 # Start time of the transient in MJD (Modified Julian Date)\n",
    "parameters['redshift'] = 0.14309633 # Redshift for 700 Mpc\n",
    "parameters['t0'] = parameters['t0_mjd_transient'] # Set the transient's starting time\n",
    "parameters['temperature_floor'] = 3000 # Minimum temperature in Kelvin (floor value) to avoid extremely low temperatures in the model\n",
    "parameters['kappa'] = 1 # Gray opacity, representing the opacity of the ejecta material\n",
    "parameters['vej'] = 0.2 # Minimum initial velocity of the ejecta in units of c (speed of light)\n",
    "parameters['ra'] = 1.0 # Right ascension in degrees for the kilonova's location\n",
    "parameters['dec'] = 1.5 # Declination in degrees for the kilonova's location\n",
    "print(parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Simulating Kilonovae for Different Cases\n",
    "\n",
    "We now simulate kilonovae for the 3 cases discussed above using the specified parameters for either the 350 Mpc or the 700 Mpc case, depending on which code block is used. This simulation utilizes the strategy designed in the corresponding CSV files.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate the transient for the Faint Kilonova case 1\n",
    "kn_sim_case_1 = SimulateOpticalTransient.simulate_transient(model='metzger_kilonova_model',\n",
    "                                       parameters=parameters, pointings_database=pointings_case_1,\n",
    "                                       survey=None, model_kwargs=model_kwargs,\n",
    "                                        end_transient_time=4, snr_threshold=5.0)\n",
    "\n",
    "# Simulate the transient for the Faint Kilonova case 2\n",
    "kn_sim_case_2 = SimulateOpticalTransient.simulate_transient(model='metzger_kilonova_model',\n",
    "                                       parameters=parameters, pointings_database=pointings_case_2,\n",
    "                                       survey=None, model_kwargs=model_kwargs,\n",
    "                                        end_transient_time=4, snr_threshold=5.0)\n",
    "\n",
    "# Simulate the transient for the Faint Kilonova case 3\n",
    "kn_sim_case_3 = SimulateOpticalTransient.simulate_transient(model='metzger_kilonova_model',\n",
    "                                       parameters=parameters, pointings_database=pointings_case_3,\n",
    "                                       survey=None, model_kwargs=model_kwargs,\n",
    "                                        end_transient_time=4, snr_threshold=5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can now print the observations that were simulated to look at the data that was generated. The data contains \n",
    "\n",
    "- `time`: \n",
    "  - The time of observation in mjd relative to the transient event start.\n",
    "\n",
    "- `magnitude`: \n",
    "  - The apparent AB magnitude of the observed object in the specified filter.\n",
    "\n",
    "- `e_magnitude`: \n",
    "  - The error or uncertainty in the measured magnitude.\n",
    "\n",
    "- `band`: \n",
    "  - The filter used for the observation. (e.g., lsstg, lsstr, lssti, etc.)\n",
    "\n",
    "- `system`: \n",
    "  - The photometric system used for magnitude calculation. In this case, it's the AB magnitude system.\n",
    "\n",
    "- `flux_density(mjy)`:\n",
    "  - The observed flux density in millijanskys (mJy).\n",
    "\n",
    "- `flux_density_error`: \n",
    "  - The uncertainty in the flux density measurement (mJy).\n",
    "\n",
    "- `flux(erg/cm2/s)`: \n",
    "  - The observed flux in units of erg/cm²/s.\n",
    "\n",
    "- `flux_error`: \n",
    "  - The uncertainty in the observed flux (erg/cm²/s).\n",
    "\n",
    "- `time (days)`: \n",
    "  - Same as the `time` column but shown with the label \"days\".\n",
    "\n",
    "- `detected`: \n",
    "  - A binary flag (1 or 0) indicating whether the transient was detected in the observation (1 for detected, 0 for non-detected).\n",
    "\n",
    "- `limiting_magnitude`: \n",
    "  - The limiting magnitude for the observation in the specified band. Observations fainter than this value may not be detected.\n",
    "\n",
    "### AB Magnitude System\n",
    "\n",
    "The AB magnitude system is a photometric system used to measure the brightness of astronomical objects. It is defined in terms of the flux density of the source in a specific bandpass. The AB magnitude $m_{\\text{AB}}$ is given by the following formula:\n",
    "\n",
    "$$\n",
    "m_{\\text{AB}} = -2.5 \\log_{10} \\left( \\frac{f_{\\nu}}{3631 \\ \\text{Jy}} \\right)\n",
    "$$\n",
    "\n",
    "\n",
    "Alternatively, with $f_{\\nu}$ still in janskys, the formula can also be written as:\n",
    "\n",
    "$$\n",
    "m_{\\text{AB}} = -2.5 \\log_{10} \\left( f_{\\nu} \\right) + 8.90\n",
    "$$\n",
    "\n",
    "Where:\n",
    "- $m_{\\text{AB}}$ is the AB magnitude.\n",
    "- $f_{\\nu}$ is the flux density of the source in a given band, measured in units of $\\text{Jy}$ (Jansky, $1 \\ \\text{Jy} = 10^{-26} \\ \\text{W} \\cdot \\text{m}^{-2} \\cdot \\text{Hz}^{-1}$).\n",
    "- $3631 \\ \\text{Jy}$ is the reference flux for an object with zero AB magnitude.\n",
    "\n",
    "### Explanation of Terms\n",
    "- **Flux density** ($f_{\\nu}$): The amount of energy per unit area per unit time per unit frequency being emitted by an astronomical source, typically in Janskys (Jy) for the AB system.\n",
    "- **Magnitude**: A logarithmic scale for measuring the brightness of an astronomical object. Lower magnitudes correspond to brighter objects.\n",
    "- **Reference Flux**: The AB magnitude system is designed so that an object with an equal flux density of $3631 \\ \\text{Jy}$ in any filter has an AB magnitude of 0.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the simulated observations for Golden Strategy Case-1\n",
    "print(kn_sim_case_1.observations)\n",
    "\n",
    "# Print the simulated observations for Golden Strategy Case-2\n",
    "print(kn_sim_case_2.observations)\n",
    "\n",
    "# Print the simulated observations for Golden Strategy Case-3\n",
    "print(kn_sim_case_3.observations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving Observations to a File\n",
    "\n",
    "We can also save the observations to a file using the `save_transient` method. This will save the observations to a CSV file in a `simulated` directory alongside the CSV file specifying the injection parameters.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the simulated observations for Golden Strategy Case-1\n",
    "kn_sim_case_1.save_transient(name='Golden_Strategy_Case_1')\n",
    "\n",
    "# Save the simulated observations for Golden Strategy Case-2\n",
    "kn_sim_case_2.save_transient(name='Golden_Strategy_Case_2')\n",
    "\n",
    "# Save the simulated observations for Golden Strategy Case-3\n",
    "kn_sim_case_3.save_transient(name='Golden_Strategy_Case_3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Calculating Filter Offsets and Plotting Simulated Light Curves\n",
    "\n",
    "Now that the data has been simulated and saved, we will define two custom functions. \n",
    "\n",
    "The first function will calculate the filter offsets for each scan based on the 5-sigma limiting magnitude of each band. Using the **etc.py** file, it will calculate the resulting exposure time in hours. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.ticker as ticker\n",
    "import etc\n",
    "\n",
    "def calculate_time_offsets_from_pointings(pointings, airmass=1.0):\n",
    "    \"\"\"\n",
    "    Calculate the filter time offsets from the manual pointings file using the fiveSigmaDepth.\n",
    "    Handles the filter mapping from 'lsstg', 'lsstr', etc., to 'g', 'r', etc.\n",
    "    \"\"\"\n",
    "    # Define the mapping from 'lsst' filters to 'u', 'g', 'r', 'i', 'z', 'y' as etc.py expects these\n",
    "    filter_mapping = {\n",
    "        'lsstu': 'u',\n",
    "        'lsstg': 'g',\n",
    "        'lsstr': 'r',\n",
    "        'lssti': 'i',\n",
    "        'lsstz': 'z',\n",
    "        'lssty': 'y'\n",
    "    }\n",
    "    \n",
    "    filter_time_offsets = {}\n",
    "    \n",
    "    for index, row in pointings.iterrows():\n",
    "        # Get the 'lsst' filter\n",
    "        lsst_filt = row['filter']\n",
    "        \n",
    "        # Map the 'lsst' filter to the corresponding 'u', 'g', 'r', etc.\n",
    "        filt = filter_mapping.get(lsst_filt, lsst_filt)  # Default to the original if no mapping found\n",
    "        \n",
    "        depth = row['fiveSigmaDepth']\n",
    "        \n",
    "        # Calculate the exposure time using the limiting magnitude (depth)\n",
    "        exptime = etc.get_exptime(depth, filt, X=airmass)\n",
    "        \n",
    "        # Convert exposure time to hours (assuming exptime is in seconds)\n",
    "        time_offset_hr = exptime / 3600\n",
    "        \n",
    "        # Store the offset for each filter, if multiple rows for the same filter, take the latest\n",
    "        filter_time_offsets[lsst_filt] = time_offset_hr  # Keep lsst format key\n",
    "    \n",
    "    return filter_time_offsets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second function will use these filter offsets, along with a visual hour delay from the transient onset and the start of the observation window, to plot the light curves simulated through the **redback** library in the Rubin ToO format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.ticker as ticker\n",
    "\n",
    "def redback_in_rubin_toO(kn_sim, visual_delay_hr, pointings, filters_color_dict, plot_title, xlim=[0, 4], ylim=[28, 18], n_interp=120):\n",
    "    \"\"\"\n",
    "    Function to plot light curves using simulated kilonova data and filter time offsets.\n",
    "    \n",
    "    Parameters:\n",
    "    - kn_sim: simulated kilonova object from redback\n",
    "    - visual_delay_hr: visual delay in hours from the transient onset\n",
    "    - pointings: pointings file with fiveSigmaDepth magnitudes\n",
    "    - filters_color_dict: dictionary mapping filters to their colors\n",
    "    - plot_title: title for the plot (should include distance and case information)\n",
    "    - xlim: x-axis limits\n",
    "    - ylim: y-axis limits\n",
    "    - n_interp: number of points for interpolation of the light curve\n",
    "    \"\"\"\n",
    "    \n",
    "    # Extract observations from kn_sim object for detected and upper limit data\n",
    "    detected_data = kn_sim.observations[kn_sim.observations['detected'] == 1.0]\n",
    "    upper_limits = kn_sim.observations[kn_sim.observations['detected'] != 1.0]\n",
    "    \n",
    "    # Calculate time offsets from the pointings file for each filter\n",
    "    filter_time_offsets = calculate_time_offsets_from_pointings(pointings)\n",
    "    print(filter_time_offsets)\n",
    "\n",
    "    # Initialize figure\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.set_ylim(ylim)\n",
    "    ax.set_xlim(xlim)\n",
    "\n",
    "    # Initialize flags for adding detection and upper limit labels only once\n",
    "    added_detection_label = False\n",
    "    added_upper_limit_label = False\n",
    "\n",
    "    # Plot the detected data points (detections) with visual delay and filter time offsets\n",
    "    for band in filters_color_dict.keys():\n",
    "        detected_band_data = detected_data[detected_data['band'] == band]\n",
    "        detected_times = detected_band_data['time (days)']\n",
    "        detected_magnitudes = detected_band_data['magnitude']\n",
    "        \n",
    "        if len(detected_times) > 0:  # Only plot detections if there are any\n",
    "            # Apply both visual delay and filter-specific time offset\n",
    "            shifted_times = detected_times + (visual_delay_hr / 24) + (filter_time_offsets.get(band, 0) / 24)\n",
    "            \n",
    "            if band == \"lsstz\":\n",
    "                if not added_detection_label:\n",
    "                    plt.scatter(shifted_times, detected_magnitudes, color=filters_color_dict[band],\n",
    "                                marker='o', s=100, edgecolor='black', linewidth=1.5, label=\"Det.\")\n",
    "                    added_detection_label = True\n",
    "                else:\n",
    "                    plt.scatter(shifted_times, detected_magnitudes, color=filters_color_dict[band],\n",
    "                                marker='o', s=100, edgecolor='black', linewidth=1.5, label=\"_nolegend_\")\n",
    "            else:\n",
    "                plt.scatter(shifted_times, detected_magnitudes, color=filters_color_dict[band],\n",
    "                            marker='o', s=100, edgecolor='black', linewidth=1.5, label=\"_nolegend_\")\n",
    "\n",
    "    # Check if the black band has no detections and add a placeholder for the legend\n",
    "    if \"lsstz\" not in detected_data['band'].values:\n",
    "        plt.scatter([], [], color='k', marker='o', s=100, edgecolor='black', linewidth=1.5, label=\"Det.\")\n",
    "\n",
    "    # Plot the non-detections (upper limits) with visual delay and filter time offsets\n",
    "    for band in filters_color_dict.keys():\n",
    "        up = upper_limits[upper_limits['band'] == band]\n",
    "        up_times = up['time (days)']\n",
    "        if len(up_times) > 0:\n",
    "            shifted_times = up_times + (visual_delay_hr / 24) + (filter_time_offsets.get(band, 0) / 24)\n",
    "            valid_times = shifted_times[shifted_times >= 0]  # Ensure no negative time values\n",
    "            valid_limits = up['limiting_magnitude'][shifted_times >= 0]\n",
    "\n",
    "            if band == \"lsstz\":\n",
    "                if not added_upper_limit_label:\n",
    "                    plt.scatter(valid_times, valid_limits, s=100, marker='v', edgecolor='blue', linewidth=1.5,\n",
    "                                color=filters_color_dict[band], label=\"UL\")\n",
    "                    added_upper_limit_label = True\n",
    "                else:\n",
    "                    plt.scatter(valid_times, valid_limits, s=100, marker='v', edgecolor='blue', linewidth=1.5,\n",
    "                                color=filters_color_dict[band], label=\"_nolegend_\")\n",
    "            else:\n",
    "                plt.scatter(valid_times, valid_limits, s=100, marker='v', edgecolor='blue', linewidth=1.5,\n",
    "                            color=filters_color_dict[band], label=\"_nolegend_\")\n",
    "\n",
    "    # Check if the black band has no upper limit and add a placeholder for the legend\n",
    "    if \"lsstz\" not in upper_limits['band'].values:\n",
    "        plt.scatter([], [], color='k', marker='v', s=100, edgecolor='blue', linewidth=1.5, label=\"UL\")\n",
    "\n",
    "    # Now plot the light curve model\n",
    "    tt = np.linspace(0, 4, n_interp)  # Adjust the time range for the model\n",
    "    parameters['output_format'] = 'magnitude' # Ensure the model output is in magnitude\n",
    "\n",
    "    for band in filters_color_dict.keys():\n",
    "        parameters['bands'] = band\n",
    "        out = redback.transient_models.kilonova_models.one_component_kilonova_model(tt, **parameters)\n",
    "        plt.plot(tt + (visual_delay_hr / 24) + (filter_time_offsets.get(band, 0) / 24), out,\n",
    "                 color=filters_color_dict[band], alpha=0.8, linewidth=2, label=f'{band[-1]}')\n",
    "\n",
    "    # Add horizontal dashed lines for each filter's fiveSigmaDepth magnitude as a reference\n",
    "\n",
    "    plt.axhline(y=24.5, color='r', linestyle='dotted', linewidth=2)\n",
    "    plt.axhline(y=23, color='yellow', linestyle='dotted', linewidth=2)\n",
    "    plt.axhline(y=25, color='g', linestyle='dotted', linewidth=2)\n",
    "\n",
    "    # Get current handles and labels\n",
    "    handles, labels = ax.get_legend_handles_labels()\n",
    "\n",
    "    # Filter out unwanted labels (like 'lssti', 'lsstg', 'lsstr'), these were appearing in the final legend when calling the plt.legend function\n",
    "    filtered_handles = []\n",
    "    filtered_labels = []\n",
    "    for handle, label in zip(handles, labels):\n",
    "        if label not in ['lssti', 'lsstg', 'lsstr']:  # Exclude these labels\n",
    "            filtered_handles.append(handle)\n",
    "            filtered_labels.append(label)\n",
    "\n",
    "    # Rebuild the legend with only the desired items and set the legend text color to black\n",
    "    plt.legend(filtered_handles, filtered_labels, loc='center left', bbox_to_anchor=(1.05, 0.5), fontsize=10, frameon=True)\n",
    "\n",
    "    # Set x-axis tick locator to intervals of 1 day\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "\n",
    "    # Customize the plot\n",
    "    plt.rcParams[\"font.family\"] = \"sans-serif\"\n",
    "    plt.rcParams[\"font.weight\"] = \"bold\"\n",
    "    plt.xlabel(\"Time from start of Rubin observing window (days)\", fontsize=13)\n",
    "    plt.ylabel(\"AB Magnitude\", fontsize=13)\n",
    "    \n",
    "    # Use the user-provided plot title\n",
    "    plt.title(plot_title, fontsize=13, weight='bold')\n",
    "\n",
    "    # Add annotation\n",
    "    plt.text(3.5 + (visual_delay_hr / 24), 19, \"Faint KN\", fontsize=20, weight='bold', color='black',\n",
    "             horizontalalignment='right', verticalalignment='top')\n",
    "\n",
    "    # Show or save the plot\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we can call the function to plot the lightcurves for different cases of the Golden Strategy \n",
    "\n",
    "> Run the below code cell if you are plotting the lightcurves for a luminosity distance of 350 Mpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 1\n",
    "\n",
    "plot_title1 = \"Gold Events: 3 Filter Deep Strategy (350 Mpc) Case 1\"\n",
    "redback_in_rubin_toO(kn_sim_case_1, visual_delay_hr=6, pointings=pointings_case_1, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title1)\n",
    "\n",
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 2\n",
    "\n",
    "plot_title2 = \"Gold Events: 3 Filter Deep Strategy (350 Mpc) Case 2\"\n",
    "redback_in_rubin_toO(kn_sim_case_2, visual_delay_hr=6, pointings=pointings_case_2, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title2)\n",
    "\n",
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 3\n",
    "\n",
    "plot_title3 = \"Gold Events: 3 Filter Deep Strategy (350 Mpc) Case 3\"\n",
    "redback_in_rubin_toO(kn_sim_case_3, visual_delay_hr=6, pointings=pointings_case_3, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Run the below code cell if you are plotting the lightcurves for a luminosity distance of 700 Mpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 1\n",
    "\n",
    "plot_title4 = \"Gold Events: 3 Filter Deep Strategy (700 Mpc) Case 1\"\n",
    "redback_in_rubin_toO(kn_sim_case_1, visual_delay_hr=6, pointings=pointings_case_1, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title4)\n",
    "\n",
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 2\n",
    "\n",
    "plot_title5 = \"Gold Events: 3 Filter Deep Strategy (700 Mpc) Case 2\"\n",
    "redback_in_rubin_toO(kn_sim_case_2, visual_delay_hr=6, pointings=pointings_case_2, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title5)\n",
    "\n",
    "# Plot the light curve for the Faint Kilonova with the Golden Strategy Case 3\n",
    "\n",
    "plot_title6 = \"Gold Events: 3 Filter Deep Strategy (700 Mpc) Case 3\"\n",
    "redback_in_rubin_toO(kn_sim_case_3, visual_delay_hr=6, pointings=pointings_case_3, filters_color_dict={'lsstg':'g', 'lsstu':'b', 'lsstr':'r', 'lsstz':'k', 'lssti':'yellow', 'lssty':'orange'}, plot_title=plot_title6)"
   ]
  }
 ],
 "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}