yarp_planner.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# yarp planner\n",
    "\n",
    "This notebook will run through the various initial setup required for running yarp. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clip level \n",
    "\n",
    "The following cell will help you in estimating the right clip levels for the various calibrators. Usually the default value should work for the standard flux calibrators like 3C147 and 3C286. This tool will be useful for choosing a right clip level for the phase calibrator using the following steps:\n",
    "\n",
    "- Choose the number of flagcal loops\n",
    "- Choose the flux level of the calibrator\n",
    "- Choose a starting upper limit in units of n-sigma\n",
    "- Choose the rate at which the cliplevel falls as a function of the \"loop-index\" using an appropriate beta.\n",
    "- If you are satisfied with the clip levels feed the above values in the `parameters.yaml` file, except for the flux of the phase calibrator, which is choosen automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clip level 1 : [-19.16, 20.84]\n",
      "Clip level 2 : [-11.471444133449165, 13.151444133449164]\n",
      "Clip level 3 : [-8.429261135439395, 10.109261135439395]\n",
      "Clip level 4 : [-6.738582832551991, 8.41858283255199]\n",
      "Clip level 5 : [-5.642626386771051, 7.32262638677105]\n"
     ]
    }
   ],
   "source": [
    "# Finds the clip levels for the different flagcal loops\n",
    "\n",
    "import yaml\n",
    "import parameters as prms\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def clip_est(fluxscale, uplim, loop, beta):\n",
    "            return [fluxscale - uplim/loop**beta, fluxscale + uplim/loop**beta]\n",
    "    \n",
    "flagcal_loops = 5 # Number of flagcal loops\n",
    "fluxlevel = 0.84\n",
    "uplim = 20 #sigma\n",
    "beta = 0.7 # powerlaw fall wrt loops\n",
    "\n",
    "for i in range(5):\n",
    "    print('Clip level',i+1,':',  clip_est(fluxlevel, uplim, i+1, beta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imaging and Self-cal setup\n",
    "\n",
    "This cell will help you in choosing a source free region in your target field using the FIRST image. This regions helps in estimating rms noise which is used for keeping a track of the improvements in the various selfcal loops. For this: \n",
    "\n",
    "- Download the FIRST and/or NVSS images of the target field.\n",
    "- Create a source-free region file using casaviewer.\n",
    "- Save the region in the regions folder declared in the `parameters.yaml` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading the FIRST image ... https://skyview.gsfc.nasa.gov/tempspace/fits/skv24643656640318_1.fits\n",
      "Downloading the NVSS image ... https://skyview.gsfc.nasa.gov/tempspace/fits/skv24643656640318_2.fits\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args='wget https://skyview.gsfc.nasa.gov/tempspace/fits/skv24643656640318_1.fits -O /mnt/work/work/HI_DATA/AGC203001_NVSS.fits', returncode=0)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from astroquery.skyview import SkyView\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "import subprocess\n",
    "import yaml\n",
    "\n",
    "target_ra = 155.51375000 # deg\n",
    "target_dec = 13.76944444  # deg\n",
    "FIRST_name = 'AGC203001_FIRST.fits'\n",
    "NVSS_name = 'AGC203001_NVSS.fits'\n",
    "\n",
    "with open('parameters.yaml') as file:\n",
    "    params = yaml.safe_load(file)\n",
    "outdir = params['general']['outdir']\n",
    "\n",
    "outdir = '/mnt/work/work/HI_DATA/'\n",
    "c = SkyCoord(target_ra, target_dec, unit=\"deg\")\n",
    "paths = SkyView.get_image_list(position=c, survey=['VLA FIRST (1.4 GHz)', 'NVSS'], radius=26*u.arcmin, pixels='1024')\n",
    "\n",
    "\n",
    "print('Downloading the FIRST image ...', paths[0])\n",
    "subprocess.run('wget {} -O {}'.format(paths[0], outdir+FIRST_name), shell=True, check=True)\n",
    "\n",
    "print('Downloading the NVSS image ...', paths[1])\n",
    "subprocess.run('wget {} -O {}'.format(paths[0], outdir+NVSS_name), shell=True, check=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuum thermal rms calculator\n",
    "\n",
    "- Find the total number of visibilities in your data (either the avspc file or full resolution file).\n",
    "- Give the full bandwidth of your data.\n",
    "- Use the output thermal rms as an estimate for choosing the various clean thresholds in the imaging-selfcal loops.\n",
    "\n",
    "Note that this is merely an estimate, true values can be far from this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Continuum Thermal rms is 0.01449373138897218 mJy\n"
     ]
    }
   ],
   "source": [
    "# Continuum thermal rms calculator\n",
    "nvis = (1524844)\n",
    "bw = ((100*1000)/8192)*10**6 # Hz\n",
    "Tsys = 45\n",
    "G = 0.22\n",
    "interval = 10.7 # secs\n",
    "\n",
    "rms = (Tsys/(G))*1/(np.sqrt(nvis*bw*interval))*10**3 # \n",
    "\n",
    "print('Continuum Thermal rms is', rms, 'mJy')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selfcal planner\n",
    "- Calculates the SNR per antenna for a given solint.\n",
    "- Also calculates the minimum solint that can be used for p and aploops.\n",
    "\n",
    "Reference: https://arxiv.org/abs/1805.05266"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMS per antenna: 0.1631686244349691 mJy  \n",
      "\n",
      "Selfcal RMS for a solint of 30 mins: 0.7297122720634482 mJy \n",
      "\n",
      "Selfcal SNR for a solint of 30 mins: 13.704031551672333 \n",
      "\n",
      "Minimum solint for phase-only selfcal: 1.1774730239999998 mins \n",
      "\n",
      "Minimum solint for amp-phase selfcal: 13.083033599999997 mins \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Calculates the SNR per antenna for a given solint.. All the formulas are from the above reference.\n",
    "\n",
    "#----------------------------------INPUTS----------------------------------------------------\n",
    "\n",
    "rms = 0.032 # in milli Jy. The final continuum rms of the image. Since in the intermediate selfcal loops \n",
    "           # this is usually not achieved, it is better to choose a more conservative value.\n",
    "    \n",
    "nants = 29 # total number of working antennas\n",
    "tonsource = 8.19 # hrs\n",
    "solint = 30 # mins\n",
    "peak_flux = 10 # mJy\n",
    "\n",
    "snr_sc_p = 3 # minimum snr required for p solutions \n",
    "snr_sc_ap = 10 # minimum snr required for ap solutions\n",
    "\n",
    "#-----------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "rms_ant = rms*np.sqrt(nants - 3) # rms per antenna\n",
    "\n",
    "\n",
    "print('RMS per antenna:', rms_ant, 'mJy  \\n')\n",
    "\n",
    "rms_sc = rms_ant*np.sqrt((10*3600)/(solint*60))\n",
    "snr_sc = peak_flux/rms_sc\n",
    "\n",
    "print('Selfcal RMS for a solint of',solint,'mins:', rms_sc,'mJy \\n')\n",
    "print('Selfcal SNR for a solint of',solint,'mins:', snr_sc,'\\n')\n",
    "\n",
    "\n",
    "# Calculates the minimum solint that can be used for a sensible p and ap solutions.\n",
    "\n",
    "\n",
    "\n",
    "max_rms_sc_p = peak_flux/snr_sc_p\n",
    "max_rms_sc_ap = peak_flux/snr_sc_ap\n",
    "\n",
    "\n",
    "min_solint_p = (rms_ant/max_rms_sc_p)**2*tonsource*60 # in mins\n",
    "min_solint_ap = (rms_ant/max_rms_sc_ap)**2*tonsource*60 # in mins\n",
    "\n",
    "print('Minimum solint for phase-only selfcal:', min_solint_p, 'mins \\n')\n",
    "print('Minimum solint for amp-phase selfcal:', min_solint_ap, 'mins \\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cube thermal rms calculator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Per channel thermal rms is 0.28343240948950205 mJy\n"
     ]
    }
   ],
   "source": [
    "# thermal rms calculator based on number of visibilities\n",
    "\n",
    "nvis = 1193288 + 274310 #301956 #+113282\n",
    "ch_width = 100/8192*10**6\n",
    "c = 299792.458 # km/s\n",
    "vel_orig = c*(ch_width/(1420.405752*10**6))\n",
    "vel_cube = 7 #km/s\n",
    "ratio = vel_cube/vel_orig\n",
    "#print(ratio)\n",
    "Tsys = 45\n",
    "G = 0.22\n",
    "interval =  10.7 #s\n",
    "\n",
    "rms = (Tsys/(G))*1/(np.sqrt(nvis*ratio*ch_width*interval))*10**3 # 10.7 is the integration time of each uvdatapoint\n",
    "\n",
    "print('Per channel thermal rms is', rms, 'mJy')"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}