{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1: Basic conversion formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## License\n",
    "\n",
    "```\n",
    "Basic conversion formulas for Radioastronomy.\n",
    "Copyright (C) 2012+  Axel Jessner (jessner@mpifr.de)\n",
    "              2015+  Benjamin Winkel (bwinkel@mpifr.de)\n",
    "\n",
    "This program is free software; you can redistribute it and/or\n",
    "modify it under the terms of the GNU General Public License\n",
    "as published by the Free Software Foundation; either version 2\n",
    "of the License, or (at your option) any later version.\n",
    "\n",
    "This program is distributed in the hope that it will be useful,\n",
    "but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n",
    "GNU General Public License for more details.\n",
    "\n",
    "You should have received a copy of the GNU General Public License\n",
    "along with this program; if not, write to the Free Software\n",
    "Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bwinkel/local/anaconda3/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n",
      "/home/bwinkel/local/anaconda3/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from astropy import constants\n",
    "from pycraf import conversions as cnv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, pycraf defines the often used log-scales and quantities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R0 = 376.730 Ohm\n",
      "Erx = 74.768 dB(uV2 / m2)\n"
     ]
    }
   ],
   "source": [
    "cnv.dimless  # == u.Unit(1)\n",
    "cnv.dB  # == dBi = u.dB(dimless)\n",
    "cnv.dB_W  # == u.dB(u.W)\n",
    "cnv.dB_W_Hz  # == u.dB(u.W / u.Hz)\n",
    "cnv.dB_W_m2  # == u.dB(u.W / u.m ** 2)\n",
    "cnv.dB_W_m2_Hz  # == u.dB(u.W / u.Hz / u.m ** 2)\n",
    "cnv.dB_Jy_Hz  # == u.dB(u.Jy * u.Hz)\n",
    "cnv.dBm  # == dB_mW = u.dB(u.mW)\n",
    "cnv.dBm_MHz  # == dB_mW_MHz = u.dB(u.mW / u.MHz)\n",
    "cnv.dB_uV_m  # == u.dB(u.uV ** 2 / u.m ** 2)\n",
    "cnv.dB_1_m  # == u.dB(1. / u.m)  # for antenna factor\n",
    "\n",
    "cnv.R0  # == 1. * (con.mu0 / con.eps0)  ** 0.5\n",
    "cnv.Erx_unit  #= (1 * u.W / 4. / np.pi * R0) ** 0.5 / (1 * u.km)\n",
    "\n",
    "print('R0 = {0.value:.3f} {0.unit}'.format(cnv.R0.to(u.ohm)))\n",
    "print('Erx = {0.value:.3f} {0.unit}'.format(cnv.Erx_unit.to(cnv.dB_uV_m)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, that mathematically, the $\\mathrm{dB}_{\\mathrm{\\mu V / m}}$ scale operates on amplitude-squares, so for astropy we need to take this into account."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Exercises\n",
    "\n",
    "### 1-Jy source\n",
    "\n",
    "- Compute the power (in Watt) that is received by a 100-m radio telescope from an object having a flux density of 1 Jy. Assume: (1) 6-cm receiver with 500 MHz of bandwidth (2) effective antenna area is 3927 m$^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f = 4.997 GHz\n",
      "S = 5.000e-18 W / m2 = -173.0 dB(W / m2)\n",
      "P_rx = 1.964e-14 W = -107.1 dB(mW)\n"
     ]
    }
   ],
   "source": [
    "wavlen = 6 * u.cm\n",
    "freq = constants.c / wavlen\n",
    "print('f = {:.3f}'.format(freq.to(u.GHz)))\n",
    "\n",
    "bandwidth = 500 * u.MHz\n",
    "S_nu = 1 * u.Jy\n",
    "Aeff = 3927 * u.m ** 2\n",
    "\n",
    "S = S_nu * bandwidth\n",
    "print('S = {:.3e} = {:.1f}'.format(\n",
    "    S.to(u.Watt / u.m ** 2), S.to(cnv.dB_W_m2)\n",
    "    ))\n",
    "\n",
    "Prx = S * Aeff\n",
    "print('P_rx = {:.3e} = {:.1f}'.format(Prx.to(u.Watt), Prx.to(cnv.dBm)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Compute the energy that is received by the telescop, if this source is observed for one minute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E = 3.273e-16 h W = 3.273e-04 h pW = 1.178e-12 J\n"
     ]
    }
   ],
   "source": [
    "int_time = 1 * u.min\n",
    "E = Prx * int_time\n",
    "\n",
    "print('E = {:.3e} = {:.3e} = {:.3e}'.format(\n",
    "    E.to(u.Watt * u.h), E.to(u.picoWatt * u.h), E.to(u.Joule)\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Compute the energy that is received by the telescope, if a 40-Jy source was observed for 40 years (the time since inauguration). By how much could one lift-up a pencil (assume a mass of 10 g) with this energy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E = 9.914e-04 J\n",
      "E = 2.754e-07 h W = 275.393 h nW\n",
      "h = 1.0 cm\n"
     ]
    }
   ],
   "source": [
    "int_time = 40 * u.year\n",
    "S_nu = 40 * u.Jy\n",
    "E = S_nu * bandwidth * Aeff * int_time\n",
    "print('E = {:.3e}'.format(E.to(u.Joule)))\n",
    "mass = 10 * u.g\n",
    "\n",
    "print('E = {:.3e} = {:.3f}'.format(E.to(u.Watt * u.h), E.to(u.nanoWatt * u.h)))\n",
    "\n",
    "\n",
    "def h_from_E(mass, E):\n",
    "    # E_pot = mass * height * g\n",
    "    return E / mass / constants.g0\n",
    "\n",
    "print('h = {:.1f}'.format(h_from_E(mass, E).to(u.cm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell phone on the moon\n",
    "\n",
    "- Compute the spectral flux density (Jy) that is emitted by a cell phone (2 W, eirp, e.g., GSM-900) on the moon. We note that such mobile devices, being time-duplexed, are only active for at most 1/8 of the time (so-called TDMA operation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S_nu = 67.3 Jy\n"
     ]
    }
   ],
   "source": [
    "dist = 384400 * u.km\n",
    "Peirp = 2 * u.Watt\n",
    "bandwidth = 200 * u.kHz\n",
    "tdma_factor = 1 / 8\n",
    "\n",
    "S = cnv.powerflux_from_ptx(tdma_factor * Peirp, dist, 0 * cnv.dBi)\n",
    "S_nu = S / bandwidth\n",
    "print('S_nu = {:.1f}'.format(S_nu.to(u.Jy)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- How much power would one receive with a parabolic dish having 30 dB gain?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wave length = 33.3 cm\n",
      "Prx = 1.189e-18 W = -149.249 dB(mW)\n"
     ]
    }
   ],
   "source": [
    "freq = 900 * u.MHz\n",
    "wavelen = constants.c / freq\n",
    "print('wave length = {:.1f}'.format(wavelen.to(u.cm)))\n",
    "\n",
    "gain = 30 * cnv.dBi\n",
    "\n",
    "Prx = cnv.prx_from_powerflux(S, freq, gain)\n",
    "print('Prx = {:.3e} = {:.3f}'.format(Prx.to(u.W), Prx.to(cnv.dBm)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- How much power would one receive with a parabolic dish of 1-m size?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pycraf import antenna"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gmax = 16.5 dB\n"
     ]
    }
   ],
   "source": [
    "gain = antenna.ras_pattern(0 * u.deg, 1 * u.m, wavelen, eta_a=50 * u.percent)\n",
    "print('Gmax = {:.1f}'.format(gain.to(cnv.dBi)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prx = 5.287e-20 W = -162.768 dB(mW)\n"
     ]
    }
   ],
   "source": [
    "Prx = cnv.prx_from_powerflux(S, freq, gain)\n",
    "print('Prx = {:.3e} = {:.3f}'.format(Prx.to(u.W), Prx.to(cnv.dBm)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Given an uncooled dual-polarization receiver (Tsys ~ 100 K), how long would one need to integrate to detect the signal? (We assume spectroscopic observations with a channel width of 5 kHz and want a 5-sigma detection in each channel, at least.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Tsys = 100 * u.K\n",
    "chan_width = 5 * u.kHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T_A[phone] = 9.6 mK\n"
     ]
    }
   ],
   "source": [
    "# calculate increase in antenna temperature (per spectral channel)\n",
    "# induced by cell phone\n",
    "Prx_nu = Prx / bandwidth\n",
    "Tphone = cnv.t_a_from_prx_nu(Prx_nu)\n",
    "print('T_A[phone] = {:.1f}'.format(Tphone.to(u.mK)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T_A[rms] = 1.9 mK\n"
     ]
    }
   ],
   "source": [
    "# want a 5-sigma detection, i.e., need Tphone / Trms > 5\n",
    "Trms = Tphone / 5\n",
    "print('T_A[rms] = {:.1f}'.format(Trms.to(u.mK)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min. integration time, tau = 75.8 h\n"
     ]
    }
   ],
   "source": [
    "# with the radiometer equation\n",
    "# Trms = Tsys / sqrt(2 * tau * bandwidth)  # factor 2 only if dual-polarization Rx\n",
    "# we can calculate the necessary integration time:\n",
    "tau = (Tsys / Trms) ** 2 / 2 / chan_width\n",
    "print('Min. integration time, tau = {:.1f}'.format(tau.to(u.hour)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "- How long does it take with a 100-m class antenna?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gmax = 56.5 dB\n"
     ]
    }
   ],
   "source": [
    "gain = antenna.ras_pattern(0 * u.deg, 100 * u.m, wavelen, eta_a=50 * u.percent)\n",
    "print('Gmax = {:.1f}'.format(gain.to(cnv.dBi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use an alternative approach, just for fun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_eff = 3927.0 m2\n"
     ]
    }
   ],
   "source": [
    "A_eff = cnv.eff_area_from_gain(gain, freq)\n",
    "print('A_eff = {:.1f}'.format(A_eff.to(u.m ** 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T_A[phone] = 95.737 K\n",
      "T_A[rms] = 19.147 K\n"
     ]
    }
   ],
   "source": [
    "Tphone = cnv.t_a_from_powerflux_nu(S_nu, A_eff)\n",
    "Trms = Tphone / 5\n",
    "print('T_A[phone] = {:.3f}'.format(Tphone.to(u.K)))\n",
    "print('T_A[rms] = {:.3f}'.format(Trms.to(u.K)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min. integration time, tau = 2.7 ms\n"
     ]
    }
   ],
   "source": [
    "tau = (Tsys / Trms) ** 2 / 2 / chan_width\n",
    "print('Min. integration time, tau = {:.1f}'.format(tau.to(u.ms)))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}