{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9a1933a5",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Calculating the GBT time required to detect a galaxy cluster with MUSTANG-2\n",
    "\n",
    "## Authors: Charles Romero and Emily Moravec\n",
    "\n",
    "Last updated: 2024-12-30.\n",
    "\n",
    "## Usage\n",
    "You can use this notebook as a template to calculate the required GBT time to detect a cluster of a given mass and at a given redshift with MUSTANG-2. You will need to install M2_ProposalTools (see https://m2-tj.readthedocs.io/en/latest/index.html for details). This notebook is based on the \"Example case of estimating observing time for assumed [A10](https://ui.adsabs.harvard.edu/abs/2010A%26A...517A..92A/abstract) clusters\" example on M2_ProposalTools (https://m2-tj.readthedocs.io/en/latest/SinglePointing_M2.html which can be downloaded at https://github.com/CharlesERomero/M2_TJ/blob/master/docs/source/SinglePointing_M2.ipynb). \n",
    "\n",
    "## About\n",
    "Hello potential MUSTANG-2 observer that wants to know how much telescope time is required to detect a given a cluster of a specific mass at a particular redshift at a particular signal-to-noise ratio. This notebook walks you through creating an expected compton-y signal map for your given parameters ($z$ and $M_{500}$), running your map through the MUSTANG-2 transfer functions to take into account the filtering of MUSTANG-2, and convolving this model map with the MUSTANG-2 beam. It then calculates the required observing time for your cluster based on your SNR requirement and the peak of the simulated SZ signal from the physical parameters you give it ($z$ and $M_{500}$).\n",
    "\n",
    "## Note\n",
    "The time estimates calculated in this notebook are made assuming using the MIDAS pipeline for data processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a571eab",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Imports and constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9c9e63f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import M2_ProposalTools.WorkHorse as WH\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "import M2_ProposalTools.FilterImages as FI\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cffd75f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constants needed for calculations - it is recommended that you keep these default parameters.\n",
    "h70            = 1                          # Normalization of Hubble parameter\n",
    "fwhm           = 10.0                       # Roughly the resolution of MUSTANG-2\n",
    "pixsize        = 2.0                        # Map pixel size (arcseconds) - making a map - package does beam convolution and filtering\n",
    "s2f            = np.sqrt(8.0*np.log(2.0))   # Conversion between FWHM and sigma\n",
    "pix_sigma      = fwhm/(pixsize*s2f)         # Gaussian sigma, in pixel size\n",
    "y2k            = -3.3                       # Approximate conversion of compton-y to Kelvin. -3.3 is conservative whereas -3.4 is perhaps more accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8308cf",
   "metadata": {},
   "source": [
    "## MUSTANG-2 Parameters\n",
    "Change these parameters based on your science. Consult https://gbtdocs.readthedocs.io/en/latest/references/receivers/mustang2/mustang2_mapping.html for the corresponding row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3963caca",
   "metadata": {},
   "outputs": [],
   "source": [
    "scansize  = 3.5 # Scan size: 3' or 3.5' is recommended \n",
    "mapping_speed_uJy = 57 # mapping speed in uJy/beam root hour\n",
    "mapping_speed_uK = 74 # mapping speed in uK root hour\n",
    "K_to_J = (mapping_speed_uJy/mapping_speed_uK) # conversion of Kelvin to Janskys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d861f603",
   "metadata": {},
   "source": [
    "## Function for calculating the observing time for a single cluster\n",
    "\n",
    "We calculate the required observing time using the radiometer equation:\n",
    "\n",
    "$t \\propto 1/\\sigma^2$ \n",
    "\n",
    "where $t$ is observing time in hours and $\\sigma$ is the sensitivity/RMS/noise of the observation.\n",
    "\n",
    "We can set up a proportional relationship with the radiometer equation such that \n",
    "\n",
    "$t_2/t_1 \\propto (\\sigma_1/\\sigma_2)^2$ \n",
    "\n",
    "where $t_2$ is the required integration time that you are solving for, $\\sigma_1$ is the RMS corresponding to your desired map size (see https://gbtdocs.readthedocs.io/en/latest/references/receivers/mustang2/mustang2_mapping.html#mustang-2-mapping-information), $t_1=1$ as the mapping speed are in within 1 hour, and $\\sigma_2$ is your desired sensitivity. \n",
    "\n",
    "Note that using the fact that mapping speed = $ms = \\sigma_1 * (t_1)^{1/2}$ and the radiometer equation in a proportion like the above you can derive that $t_2 = (ms / \\sigma_2)^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a2a9416a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_observing_time(M500,z,SNR_required):\n",
    "    print(\"-\"*5,\"Cluster Properties\",\"-\"*5)\n",
    "    print(\"M500 (1e14):\",(M500/1e14).value)\n",
    "    print(\"Redshift:\",z)\n",
    "\n",
    "    # Create expected signal map\n",
    "    ymap = WH.make_A10Map(M500,z,pixsize=pixsize,Dist=True) # create compton-y map\n",
    "    mymap = WH.smooth_by_M2_beam(ymap,pixsize=pixsize) # convolve signal map (y) with MUSTANG-2 beam\n",
    "\n",
    "    # Filter signal map\n",
    "    tab = WH.get_xfertab(scansize) # 2D array of values used to create filtered map\n",
    "    filtered_model_ymap = FI.apply_xfer(mymap,tab,pixsize) # apply transfer function to compton-y map and create a filtered map\n",
    "    snr_smoothing_kernel = pix_sigma*0.9 # when an SNR map is created there is smoothing applied - it is recommended to keep this smoothing at 90% of M2 beam (so keep 0.9)\n",
    "    yxfer = scipy.ndimage.gaussian_filter(filtered_model_ymap, snr_smoothing_kernel)\n",
    "    \n",
    "    # Convert model map to observable values (uJy/beam or uK)\n",
    "    uKmap      = yxfer*y2k*1e6 # convert compton-y map to uK (1e6 is to convert from K to uK)\n",
    "    uJymap     = uKmap * K_to_J\n",
    "    \n",
    "    # Calculate peak value wrt the scale of interest (some values are negative)\n",
    "    print(\"-\"*5,\"Signal Peak Stats\",\"-\"*5)\n",
    "    SZpeak_yraw  = np.max(ymap) # peak of y-map\n",
    "    print(\"Unfiltered y peak:\", SZpeak_yraw)\n",
    "    SZpeak_yconvolved  = np.max(mymap)\n",
    "    print(\"Beam-convolved y peak:\", SZpeak_yconvolved)\n",
    "    SZpeak_yxfer  = np.max(yxfer)\n",
    "    print(\"Filtered y peak:\", SZpeak_yxfer)\n",
    "    SZpeak_uK  = np.min(uKmap) # this is negative so calculate minimum\n",
    "    print(\"Filtered peak, uK:\", SZpeak_uK)\n",
    "    SZpeak_uJy = np.min(uJymap) # this is negative so calculate minimum\n",
    "    print(\"Filtered peak, uJy/beam:\", SZpeak_uJy)\n",
    "    \n",
    "    # Calculate the observing time based on the Compton-y peak in uJy/beam\n",
    "    print(\"-\"*5,\"Telescope Time Stats\",\"-\"*5)\n",
    "    target_sensitivity = SZpeak_uJy/SNR_required # calculate the detection threshold based on the required SNR\n",
    "    print(\"Your target sensitivity or 1-sigma level in uJy/beam is:\",np.abs(target_sensitivity))\n",
    "    # from the radiometer and a proportional relationship (see math in comments above) we can calculate the observing time required\n",
    "    ObsTime = np.round((mapping_speed_uJy/target_sensitivity)**2,5) # from unit analysis this variable is in units of hours\n",
    "    print(\"On-source time to\",str(SNR_required),\"sigma peak detection (hrs):\",ObsTime)\n",
    "    TotalTimeRequest = ObsTime * 2 # include the required GBT high frequency factor of 2 to account for overheads\n",
    "    # GBT time gets scheduled in units of 0.25 hr so convert and \n",
    "    TelTime = np.ceil(TotalTimeRequest*4)/4.0\n",
    "\n",
    "    print(\"Total telescope time for this cluster (hrs):\",TelTime)\n",
    "    \n",
    "    return TelTime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5e073f4",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Calculate the telescope time for one cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "267d3e7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----- Cluster Properties -----\n",
      "M500 (1e14): 3.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 9.619897496889924e-05\n",
      "Beam-convolved y peak: 8.860120011199452e-05\n",
      "Filtered y peak: 4.829389221185476e-05\n",
      "Filtered peak, uK: -159.3698442991207\n",
      "Filtered peak, uJy/beam: -122.7578530412146\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 17.536836148744943\n",
      "On-source time to 7 sigma peak detection (hrs): 10.56446\n",
      "Total telescope time for this cluster (hrs): 21.25\n"
     ]
    }
   ],
   "source": [
    "# Physical parameters of your cluster\n",
    "M500 = 3*1e14*u.M_sun # Mass of cluster\n",
    "z = 1.0 # The redshift of the cluster\n",
    "obstime_single_cluster = calculate_observing_time(M500=M500,z=z,SNR_required=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa1752e0",
   "metadata": {},
   "source": [
    "Why do we require the SNR=7? Based on experience the MUSTANG-2 instrument team highly recommends that you adopt a $7\\sigma$ threshold, but depending on your exact science objective, a different threshold may be adopted. Feel free to discuss this with the instrument team.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34680d9a",
   "metadata": {},
   "source": [
    "## Calculate the telescope time for multiple clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6ee0fbe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Physical parameters of your cluster\n",
    "M500s = np.array([1,2,3,4,5,6,7,8])*1e14*u.M_sun # Mass of clusters\n",
    "zs = np.array([1.0,1.0,1.0,1.0,1.0,1.0,1.0,1]) # Matching redshifts for the M500s list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f1dc82b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 1.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 2.7103546504401935e-05\n",
      "Beam-convolved y peak: 2.4135944165655142e-05\n",
      "Filtered y peak: 1.4965757817167562e-05\n",
      "Filtered peak, uK: -49.387000796652956\n",
      "Filtered peak, uJy/beam: -38.04133845147592\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 5.434476921639417\n",
      "On-source time to 7 sigma peak detection (hrs): 110.01052\n",
      "Total telescope time for this cluster (hrs): 220.25\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 2.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 6.0365334122178264e-05\n",
      "Beam-convolved y peak: 5.530305199658348e-05\n",
      "Filtered y peak: 3.1817993020714536e-05\n",
      "Filtered peak, uK: -104.99937696835796\n",
      "Filtered peak, uJy/beam: -80.87789847562708\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 11.553985496518154\n",
      "On-source time to 7 sigma peak detection (hrs): 24.33807\n",
      "Total telescope time for this cluster (hrs): 48.75\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 3.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 9.619897496889924e-05\n",
      "Beam-convolved y peak: 8.860120011199452e-05\n",
      "Filtered y peak: 4.829389221185476e-05\n",
      "Filtered peak, uK: -159.3698442991207\n",
      "Filtered peak, uJy/beam: -122.7578530412146\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 17.536836148744943\n",
      "On-source time to 7 sigma peak detection (hrs): 10.56446\n",
      "Total telescope time for this cluster (hrs): 21.25\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 4.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 0.00013376973214360876\n",
      "Beam-convolved y peak: 0.00012461706064739533\n",
      "Filtered y peak: 6.52603290828719e-05\n",
      "Filtered peak, uK: -215.35908597347722\n",
      "Filtered peak, uJy/beam: -165.88470135794867\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 23.697814479706953\n",
      "On-source time to 7 sigma peak detection (hrs): 5.7854\n",
      "Total telescope time for this cluster (hrs): 11.75\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 5.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 0.00017266723341949052\n",
      "Beam-convolved y peak: 0.00016159098229221895\n",
      "Filtered y peak: 8.1695551395762e-05\n",
      "Filtered peak, uK: -269.5953196060146\n",
      "Filtered peak, uJy/beam: -207.66125969652475\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 29.66589424236068\n",
      "On-source time to 7 sigma peak detection (hrs): 3.69177\n",
      "Total telescope time for this cluster (hrs): 7.5\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 6.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 0.00021264462426804625\n",
      "Beam-convolved y peak: 0.00019890428730061485\n",
      "Filtered y peak: 9.761530156346558e-05\n",
      "Filtered peak, uK: -322.1304951594364\n",
      "Filtered peak, uJy/beam: -248.12754356875507\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 35.44679193839358\n",
      "On-source time to 7 sigma peak detection (hrs): 2.58581\n",
      "Total telescope time for this cluster (hrs): 5.25\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 7.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 0.0002535355941774716\n",
      "Beam-convolved y peak: 0.00023785709856947226\n",
      "Filtered y peak: 0.0001135652820943107\n",
      "Filtered peak, uK: -374.7654309112253\n",
      "Filtered peak, uJy/beam: -288.6706697559438\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 41.238667107991965\n",
      "On-source time to 7 sigma peak detection (hrs): 1.91047\n",
      "Total telescope time for this cluster (hrs): 4.0\n",
      "==================================================\n",
      "----- Cluster Properties -----\n",
      "M500 (1e14): 8.0\n",
      "Redshift: 1.0\n",
      "----- Signal Peak Stats -----\n",
      "Unfiltered y peak: 0.00029521992026628467\n",
      "Beam-convolved y peak: 0.0002786569572687501\n",
      "Filtered y peak: 0.0001299492606572224\n",
      "Filtered peak, uK: -428.8325601688339\n",
      "Filtered peak, uJy/beam: -330.3169720219396\n",
      "----- Telescope Time Stats -----\n",
      "Your target sensitivity or 1-sigma level in uJy/beam is: 47.188138860277085\n",
      "On-source time to 7 sigma peak detection (hrs): 1.4591\n",
      "Total telescope time for this cluster (hrs): 3.0\n",
      "==================================================\n",
      "Total Telescope Time (hrs):  321.75\n"
     ]
    }
   ],
   "source": [
    "TotalTime = 0 # Initiate total time counter\n",
    "for M500,z in zip(M500s,zs):\n",
    "    print(\"=\"*50)\n",
    "    obstime_cluster = calculate_observing_time(M500=M500,z=z,SNR_required=7)\n",
    "    TotalTime += obstime_cluster    \n",
    "print(\"=\"*50)\n",
    "print(\"Total Telescope Time (hrs): \", TotalTime)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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": 5
}