{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example case of estimating observing time for assumed A10 clusters #\n",
    "\n",
    "This thread shows how to use _M2_ProposalTools_ to estimate the required time to observe clusters based on a single pointing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": "markdown",
   "metadata": {},
   "source": [
    "Let's just take the case of a few masses at high redshift ($z=1.0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "M500s = np.array([3.0,4.0,5.0,6.0])*1e14*u.M_sun\n",
    "lgM200   = 15.1                    # Log_10 M200\n",
    "M200     = np.power(10,15.1)       # Now in normal units\n",
    "###How to go to M500?\n",
    "M500s     = np.array([0.7]) * M200*u.M_sun \n",
    "zs    = np.array([0.411]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "h70            = 1                          # Normalization of Hubble parameter\n",
    "fwhm           = 10.0                       # Roughly the resolution of MUSTANG-2\n",
    "pixsize        = 2.0                        # Assume this pixel size\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",
    "scansize       = 3.0                        # Assume scan sizes of 3 arcminutes\n",
    "y2k            = -3.3                       # Approximate conversion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll just dive into some of the functions and get what we want: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M500 (1e14):  8.812477882559163\n",
      "Redshift:  0.411\n",
      "Y500 (1e12 sr):  111.41646803767371\n",
      "Unfiltered y peak:  0.0001712578226715335\n",
      "Beam-convolved y peak:  0.00016522446964474616\n",
      "Filtered y peak:  5.078609894123501e-05\n",
      "Filtered peak, uK:  -167.59412650607553\n",
      "Filtered peak, uJy/beam:  -128.56535731972917\n",
      "Time to 10 sigma peak detection (hrs):  19.108416934318672\n",
      "================================================\n",
      "Total Telescope Time:  38.25\n"
     ]
    }
   ],
   "source": [
    "TotalTime=0\n",
    "for M500,z in zip(M500s,zs):\n",
    "\n",
    "    y500       = WH.y_delta_from_mdelta(M500.value,z,delta=500,YMrel=\"A10\",ycyl=False,h70=h70)\n",
    "    \n",
    "    ymap       = WH.make_A10Map(M500,z,pixsize=pixsize,Dist=True)\n",
    "    mymap      = WH.smooth_by_M2_beam(ymap,pixsize=pixsize)\n",
    "    \n",
    "    tab        = WH.get_xfertab(scansize)\n",
    "    yxfer_init = FI.apply_xfer(mymap,tab,pixsize)\n",
    "    yxfer      = scipy.ndimage.filters.gaussian_filter(yxfer_init, pix_sigma*0.9)\n",
    "    uKmap      = yxfer*y2k*1e6\n",
    "    uJymap     = uKmap * (56/73.0)\n",
    "    MS         = 56.2\n",
    "    \n",
    "    #print(\"pixel size: \",pixsize)\n",
    "    SZpeak_uK  = np.min(uKmap)\n",
    "    SZpeak_uJy = np.min(uJymap)\n",
    "    SZpeak_yr  = np.max(ymap)\n",
    "    SZpeak_yc  = np.max(mymap)\n",
    "    SZpeak_yx  = np.max(yxfer)\n",
    "    ObsTime    = (MS*10/SZpeak_uJy)**2\n",
    "    TelTime    = np.ceil(ObsTime*8)/4.0\n",
    "    print(\"M500 (1e14): \",(M500/1e14).value)\n",
    "    print(\"Redshift: \",z)\n",
    "    #print(\"R500 (arcmin): \",Theta500*60*180/np.pi)\n",
    "    print(\"Y500 (1e-12 sr): \",y500*1e12)\n",
    "    print(\"Unfiltered y peak: \", SZpeak_yr)\n",
    "    print(\"Beam-convolved y peak: \", SZpeak_yc)\n",
    "    print(\"Filtered y peak: \", SZpeak_yx)\n",
    "    print(\"Filtered peak, uK: \", SZpeak_uK)\n",
    "    print(\"Filtered peak, uJy/beam: \", SZpeak_uJy)\n",
    "    print(\"Time to 10 sigma peak detection (hrs): \",ObsTime)\n",
    "    print(\"================================================\")\n",
    "    TotalTime += TelTime\n",
    "    \n",
    "print(\"Total Telescope Time: \", TotalTime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}