Observation Planning Part 2 Exercise - Sophie Clark.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sophie Clark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Observation Planning Part 2 - Sensitivity Calculation\n",
    "\n",
    "1. Where is my target?\n",
    "2. When can I observe my target?\n",
    "3. How do I know when I've found my target (make a finder chart)?\n",
    "4. __How long do I need to observe?__\n",
    "5. How will I calibrate my data?\n",
    "\n",
    "In Observation Planning part 1, we covered parts 1,2,3 above.  Now we'll cover the 4th part, calculating observation times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, a brief review:  For each of the following questions, give a brief explanation of what tools you use to answer the first three questions:\n",
    "\n",
    "(1) Where is my target?\n",
    "\n",
    "We can query a database using astropy to get its RA and Dec.\n",
    "\n",
    "(2) When can I observe my target?\n",
    "\n",
    "Use astroplan to determine the sunrise and sunset time at your location, and then create an airmass plot to determine when each source will be visible. You can also create star charts to see where the source will be at different times during the observing run.\n",
    "\n",
    "(3) How do I know when I've found my target?\n",
    "\n",
    "Compare the image received with a finder image, obtained by using astroquery to obtain data from SkyView.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How long do I need to observe my target?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know approximately how long we _can_ observe each target, but how long do we _need_ per source?  To answer this question, we need to do some sensitivity calculations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we determine how bright our target is.  You can only obtain this information for known sources, but generally, we can set a reasonable limit based off of physical considerations.  In other words, as long as we know what we're looking for, we can make a pretty good guess of how long we need to look."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our observations will be performed in the visual, meaning we are interested in the source brightness in the R (red) and V (visual) filters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use SIMBAD to look up the typical brightness of the target source in magnitudes.  You can go on the [SIMBAD website](http://simbad.u-strasbg.fr/simbad/) to see how this works, but astroquery can also do it, so we're going to use astroquery here.  The astroquery SIMBAD module tells you what \"fields\" (data column names) are available to query on for each source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.simbad import Simbad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall from the previous workbook that we can use SIMBAD to retrieve information about stars, but the information included does not, by default, include the flux (or magnitude) of the source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=1</i>\n",
       "<table id=\"table1568908050240\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>MAIN_ID</th><th>RA</th><th>DEC</th><th>RA_PREC</th><th>DEC_PREC</th><th>COO_ERR_MAJA</th><th>COO_ERR_MINA</th><th>COO_ERR_ANGLE</th><th>COO_QUAL</th><th>COO_WAVELENGTH</th><th>COO_BIBCODE</th><th>SCRIPT_NUMBER_ID</th></tr></thead>\n",
       "<thead><tr><th></th><th>&quot;h:m:s&quot;</th><th>&quot;d:m:s&quot;</th><th></th><th></th><th>mas</th><th>mas</th><th>deg</th><th></th><th></th><th></th><th></th></tr></thead>\n",
       "<thead><tr><th>object</th><th>str13</th><th>str13</th><th>int16</th><th>int16</th><th>float32</th><th>float32</th><th>int16</th><th>str1</th><th>str1</th><th>object</th><th>int32</th></tr></thead>\n",
       "<tr><td>* P Cyg</td><td>20 17 47.2020</td><td>+38 01 58.552</td><td>11</td><td>11</td><td>0.120</td><td>0.130</td><td>90</td><td>A</td><td>O</td><td>2007A&amp;A...474..653V</td><td>1</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=1>\n",
       "MAIN_ID       RA           DEC      ...     COO_BIBCODE     SCRIPT_NUMBER_ID\n",
       "           \"h:m:s\"       \"d:m:s\"    ...                                     \n",
       " object     str13         str13     ...        object            int32      \n",
       "------- ------------- ------------- ... ------------------- ----------------\n",
       "* P Cyg 20 17 47.2020 +38 01 58.552 ... 2007A&A...474..653V                1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcyg_simbad = Simbad.query_object('P Cygni')\n",
    "pcyg_simbad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--NOTES--\n",
      "\n",
      "1. The parameter filtername must correspond to an existing filter. Filters include: B,V,R,I,J,K.  They are checked by SIMBAD but not astroquery.simbad\n",
      "\n",
      "2. Fields beginning with rvz display the data as it is in the database. Fields beginning with rv force the display as a radial velocity. Fields beginning with z force the display as a redshift\n",
      "\n",
      "3. For each measurement catalog, the VOTable contains all fields of the first measurement. When applicable, the first measurement is the mean one. \n",
      "\n",
      "Available VOTABLE fields:\n",
      "\n",
      "bibcodelist(y1-y2)\n",
      "biblio\n",
      "cel\n",
      "cl.g\n",
      "coo(opt)\n",
      "coo_bibcode\n",
      "coo_err_angle\n",
      "coo_err_maja\n",
      "coo_err_mina\n",
      "coo_qual\n",
      "coo_wavelength\n",
      "coordinates\n",
      "dec(opt)\n",
      "dec_prec\n",
      "diameter\n",
      "dim\n",
      "dim_angle\n",
      "dim_bibcode\n",
      "dim_incl\n",
      "dim_majaxis\n",
      "dim_minaxis\n",
      "dim_qual\n",
      "dim_wavelength\n",
      "dimensions\n",
      "distance\n",
      "distance_result\n",
      "einstein\n",
      "fe_h\n",
      "flux(filtername)\n",
      "flux_bibcode(filtername)\n",
      "flux_error(filtername)\n",
      "flux_name(filtername)\n",
      "flux_qual(filtername)\n",
      "flux_system(filtername)\n",
      "flux_unit(filtername)\n",
      "fluxdata(filtername)\n",
      "gcrv\n",
      "gen\n",
      "gj\n",
      "hbet\n",
      "hbet1\n",
      "hgam\n",
      "id(opt)\n",
      "ids\n",
      "iras\n",
      "irc\n",
      "iso\n",
      "iue\n",
      "jp11\n",
      "link_bibcode\n",
      "main_id\n",
      "measurements\n",
      "membership\n",
      "mesplx\n",
      "mespm\n",
      "mk\n",
      "morphtype\n",
      "mt\n",
      "mt_bibcode\n",
      "mt_qual\n",
      "otype\n",
      "otype(opt)\n",
      "otypes\n",
      "parallax\n",
      "plx\n",
      "plx_bibcode\n",
      "plx_error\n",
      "plx_prec\n",
      "plx_qual\n",
      "pm\n",
      "pm_bibcode\n",
      "pm_err_angle\n",
      "pm_err_maja\n",
      "pm_err_mina\n",
      "pm_qual\n",
      "pmdec\n",
      "pmdec_prec\n",
      "pmra\n",
      "pmra_prec\n",
      "pos\n",
      "posa\n",
      "propermotions\n",
      "ra(opt)\n",
      "ra_prec\n",
      "rot\n",
      "rv_value\n",
      "rvz_bibcode\n",
      "rvz_error\n",
      "rvz_qual\n",
      "rvz_radvel\n",
      "rvz_type\n",
      "rvz_wavelength\n",
      "sao\n",
      "sp\n",
      "sp_bibcode\n",
      "sp_nature\n",
      "sp_qual\n",
      "sptype\n",
      "td1\n",
      "typed_id\n",
      "ubv\n",
      "uvby\n",
      "uvby1\n",
      "v*\n",
      "velocity\n",
      "xmm\n",
      "z_value\n",
      "For more information on a field:\n",
      "Simbad.get_field_description ('field_name') \n",
      "Currently active VOTABLE fields:\n",
      " ['main_id', 'coordinates']\n"
     ]
    }
   ],
   "source": [
    "Simbad.list_votable_fields()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add the flux fields for the R and V filters, then we can see the two new columns added on the right:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=1</i>\n",
       "<table id=\"table1568911098688\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>MAIN_ID</th><th>RA</th><th>DEC</th><th>RA_PREC</th><th>DEC_PREC</th><th>COO_ERR_MAJA</th><th>COO_ERR_MINA</th><th>COO_ERR_ANGLE</th><th>COO_QUAL</th><th>COO_WAVELENGTH</th><th>COO_BIBCODE</th><th>FLUX_R</th><th>FLUX_V</th><th>SCRIPT_NUMBER_ID</th></tr></thead>\n",
       "<thead><tr><th></th><th>&quot;h:m:s&quot;</th><th>&quot;d:m:s&quot;</th><th></th><th></th><th>mas</th><th>mas</th><th>deg</th><th></th><th></th><th></th><th>mag</th><th>mag</th><th></th></tr></thead>\n",
       "<thead><tr><th>object</th><th>str13</th><th>str13</th><th>int16</th><th>int16</th><th>float32</th><th>float32</th><th>int16</th><th>str1</th><th>str1</th><th>object</th><th>float32</th><th>float32</th><th>int32</th></tr></thead>\n",
       "<tr><td>* P Cyg</td><td>20 17 47.2020</td><td>+38 01 58.552</td><td>11</td><td>11</td><td>0.120</td><td>0.130</td><td>90</td><td>A</td><td>O</td><td>2007A&amp;A...474..653V</td><td>4.28</td><td>4.82</td><td>1</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=1>\n",
       "MAIN_ID       RA           DEC      RA_PREC ...  FLUX_R  FLUX_V SCRIPT_NUMBER_ID\n",
       "           \"h:m:s\"       \"d:m:s\"            ...   mag     mag                   \n",
       " object     str13         str13      int16  ... float32 float32      int32      \n",
       "------- ------------- ------------- ------- ... ------- ------- ----------------\n",
       "* P Cyg 20 17 47.2020 +38 01 58.552      11 ...    4.28    4.82                1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Simbad.reset_votable_fields()\n",
    "Simbad.add_votable_fields('flux(R)', 'flux(V)')\n",
    "pcyg_simbad = Simbad.query_object('P Cygni')\n",
    "pcyg_simbad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "P Cygni is about 4.3 magnitudes in R.  What does that mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Magnitudes\n",
    "\n",
    "\n",
    "Magnitudes are a weird system with relative definitions.  The definition of magnitude is:\n",
    "\n",
    "$$ m = -2.5 \\log(F / F_0) $$\n",
    "\n",
    "$m$ is the observed (apparent) magnitude, $F$ is the source flux (in units of Jansky, ergs/cm$^2$/angstrom, W/m$^2$/angstrom, or erg/cm$^{2}$/Hz, usually), and $F_0$ is the flux at magnitude zero. \n",
    "\n",
    "__$F_0$ is a definition!__\n",
    "\n",
    "It was originally referenced to Vega, but more accurate observations have revealed that Vega is itself variable, so now the absolute standard is set by physical units.\n",
    "\n",
    "So, we have to look up $F_0$ or determine it ourselves.\n",
    "\n",
    "Note that each filter will have its own zero-point (also, terminology aside: the \"zero point of a filter\" and the \"photometric zero point\" are different things).  For example, for the visible $V$ band, we say:\n",
    "\n",
    "$$ m_V = -2.5 \\log(F_V / F_{0,V}) $$\n",
    "or for red:\n",
    "$$ m_R = -2.5 \\log(F_R / F_{0,R}) $$\n",
    "\n",
    "In the red filter, the zero-point is approximately 2500 Jy (which we just look up).  Therefore, we can compute the flux from P Cygni:\n",
    "\n",
    "$$ m_{P Cyg} = 4.3 = -2.5 \\log(F_{P Cyg} / 2500 \\mathrm{Jy}) $$\n",
    "\n",
    "Solve for $F_{P Cyg}$:\n",
    "\n",
    "$$ F_{P Cyg} = 10^{(-4.3/2.5)} \\times 2500 \\mathrm{Jy} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy import constants\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$47.636518 \\; \\mathrm{Jy}$"
      ],
      "text/plain": [
       "<Quantity 47.63651795 Jy>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_PCyg = 10**(-4.3/2.5) * 2500*u.Jy\n",
    "F_PCyg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tells us the amount of energy received per area per unit frequency at the telescope.  Next, we need to compute the number of photons we'll collect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1 Jansky is $10^{-23}$ ergs s$^{-1}$ cm$^{-2}$ Hz$^{-1}$ (https://en.wikipedia.org/wiki/Jansky)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1 \\times 10^{-23} \\; \\mathrm{\\frac{erg}{Hz\\,s\\,cm^{2}}}$"
      ],
      "text/plain": [
       "<Quantity 1.e-23 erg / (cm2 Hz s)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1*u.Jy).to(u.erg/u.s/u.cm**2/u.Hz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modeling the telescope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to describe our instrument now.  Our telescope is a 14\" diameter telescope, so its area is (pi * (radius)$^2$).  An inch is 2.54 centimeters.  So, our collecting area is $A_{CTO 14\"}=$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$993.14666 \\; \\mathrm{cm^{2}}$"
      ],
      "text/plain": [
       "<Quantity 993.14665903 cm2>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer here\n",
    "A_CTO = (np.pi)*((14/2)*2.54*u.cm)**2\n",
    "A_CTO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This should be 993.14666 cm$^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modeling the Filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine how many photons we will receive, we need to know which wavelengths of photons are hitting the detector, and which ones are energetic enough to produce electrons through the photoelectric effect.\n",
    "\n",
    "We do this by putting a filter in front of the detector, which prevents longer- or shorter-wavelength (lower- or higher-energy) photons from getting through."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will approximate that our red filter is $\\sim$1000 angstroms wide (in reality, there is a [wavelength-dependent filter function we could use](http://voservices.net/filter/), but we will take the simple approach for now).  The width of a filter is relative to the wavelength, and the red filter is centered at about 6500 angstroms.\n",
    "\n",
    "To convert a central wavelength from angstroms to Hertz, we use the definition $\\nu  \\lambda = c$, i.e., the frequency times the wavelength is the speed of light. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, 100 microns ($\\mu$m) is at 3 THz.  Astropy's `constants.c` is the speed of light:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2997924.6 \\; \\mathrm{\\frac{m}{s\\,\\mu m}}$"
      ],
      "text/plain": [
       "<Quantity 2997924.58 m / (s um)>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wavelength = 100*u.um\n",
    "frequency = constants.c  / wavelength\n",
    "frequency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the units on that are something funky, $m \\mu^{-1} s^{-1}$.  You can use astropy's `.to` to convert to a specific unit, or `.decompose` to reduce to a minimal expressible unit.  For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Decomposed frequency:  2997924580000.0 1 / s\n",
      "Frequency in THz:  2.99792458 THz\n"
     ]
    }
   ],
   "source": [
    "print(\"Decomposed frequency: \", frequency.decompose())\n",
    "print(\"Frequency in THz: \", frequency.to(u.THz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, returning to our red filter, for a wavelength $\\lambda = 6500$ angstroms, solve the equation $\\lambda \\nu = c$ for $\\nu$ in Hz or THz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$4.6121917 \\times 10^{14} \\; \\mathrm{Hz}$"
      ],
      "text/plain": [
       "<Quantity 4.61219166e+14 Hz>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer here\n",
    "nu = (constants.c)/(6500e-10*u.m)\n",
    "nu.to(u.Hz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should be $4.6121917×10^{14}$ Hz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filter width is $\\Delta\\lambda / \\lambda$.  We have defined the width $\\Delta \\lambda = 1000$ angstroms, and the central wavelength is $\\lambda=6500$ angstroms.  The fractional width of the filter is the same in both units, i.e.:\n",
    "$$ \\frac{\\Delta\\lambda}{\\lambda} = \\frac{\\Delta\\nu}{\\nu} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You solved for $\\nu$ above, so now determine $\\Delta \\nu$ (again, in Hz or THz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$7.0956795 \\times 10^{13} \\; \\mathrm{Hz}$"
      ],
      "text/plain": [
       "<Quantity 7.09567948e+13 Hz>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer here\n",
    "dlam = 1000*u.AA\n",
    "lam = 6500*u.AA # note we can't use the variable \"lambda\" because it's a reserved word in python\n",
    "dnu = (dlam/lam)*((constants.c)/(6500e-10*u.m))\n",
    "dnu.to(u.Hz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your answer should be $7.0956795×10^{13}$Hz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can convert our flux in Janskys to $10^{-23}$ ergs s$^{-1}$ cm$^{-2}$ Hz$^{-1}$, then multiply by the bandwidth in Hz and the collecting area in cm$^{2}$ to obtain the energy we will receive on the detector per second:\n",
    "\n",
    "$$ E_{rec} = F_{P Cyg} * A_{CTO 14\"} * \\Delta \\nu $$\n",
    "\n",
    "Compute $E_{rec}$ in units of erg $s^{-1}$ (`u.erg/u.s`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Quantity 7.09567948e+13 1 / s>,\n",
       " <Quantity 47.63651795 Jy>,\n",
       " <Quantity 993.14665903 cm2>)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# you should have these three variables defined already\n",
    "dnu, F_PCyg, A_CTO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$3.3569694 \\times 10^{-5} \\; \\mathrm{\\frac{erg}{s}}$"
      ],
      "text/plain": [
       "<Quantity 3.35696941e-05 erg / s>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer here (in erg/s)\n",
    "E_rec = F_PCyg*A_CTO*dnu\n",
    "E_rec.to(u.erg/u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should be: $3.3569694×10^{−5}$ ergs / s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CCD detector responds to individual photons, so we need to count how many photons are coming in.  The average wavelegnth of a photon is 6500 angstroms (approximately).  We can compute the energy per photon using\n",
    "\n",
    "$$ E_{phot} = h \\nu$$\n",
    "\n",
    "where $h$ is [Planck's constant](https://en.wikipedia.org/wiki/Planck_constant).  Using the above information, compute the mean energy per photon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(recall that frequency and wavelength are related as $c = \\nu / \\lambda$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.9074492066646196 eV\n",
      "3.05607054945989e-12 erg\n"
     ]
    }
   ],
   "source": [
    "# student answer here (compute the anwer in both eV and erg)\n",
    "h = constants.h\n",
    "E_phot = h*nu\n",
    "print(E_phot.to(u.eV))\n",
    "print(E_phot.to(u.erg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should obtain a number near 2 eV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given our energy received from the star $F_{E} =$ `E_rec` [erg/s] and the energy per photon $E_{phot}$ = `E_phot`[erg], how many photons are we receiving per second $F_{\\gamma}$?\n",
    "\n",
    "$F_{\\gamma} [\\mathrm{phot~s}^{-1}] = F_{E} / E_{phot}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$10984594 \\; \\mathrm{\\frac{1}{s}}$"
      ],
      "text/plain": [
       "<Quantity 10984593.97222623 1 / s>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer here (units (1/u.s))\n",
    "phot_per_s = (E_rec.to(u.erg/u.s))/(E_phot.to(u.erg))\n",
    "phot_per_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer should be: 10984594 s$^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CCD detectors saturate at approximately $2^{16} = 65536$ counts (though ours saturate at closer to 30,000).  Is the photon count greater than or less than this number? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The photon count is less than this number, since the number of counts indicates both photons and electrons. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above is greater than the saturation count, how many pixels do we need to spread the photons over if we want the image to stay below the saturation limit in a 1 second exposure?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$167.6116 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 167.61160236>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer\n",
    "npix = phot_per_s/(65536*(1/u.s))\n",
    "npix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CCDs at CTO are ST-402 cameras.  They have 765x510 pixels over 6.9x4.3mm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How big are the individual pixels, in microns? (should be about 9x9 microns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.019607843137257 um\n",
      "8.431372549019608 um\n"
     ]
    }
   ],
   "source": [
    "# student answer: compute the pixel size in both dimensions\n",
    "xsize = (6.9*u.mm)/(765)\n",
    "ysize = (4.3*u.mm)/(510)\n",
    "print(xsize.to(u.um))\n",
    "print(ysize.to(u.um))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By observing a pair of stars with known separation on the sky, we can determine what the _pixel scale_ of our CCD is.  \n",
    "\n",
    "However, for the purpose of this exercise, we will simply assume that the pixel scale is $\\approx0.5\"$ per pixel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Over what radius circle does our star need to be projected to avoid saturating any pixels?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate this in two steps:  First, we know that the _area_ we need to project onto must be ~168 pixels.  What _radius_ circle has an area of 168 pixels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.312732791431452"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer\n",
    "radius_pixels = (168/np.pi)**(1/2)\n",
    "radius_pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should be 7.3 pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does that size in pixels correspond to in arcseconds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the \"pixel scale\" is the number of arcseconds across a pixel,\n",
    "# so the units are \"arcseconds per pixel\",\n",
    "# but we leave the \"per pixel\" part implicit (otherwise the code gets confusing)\n",
    "pixel_scale = 0.5*u.arcsec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$3.6563664 \\; \\mathrm{{}^{\\prime\\prime}}$"
      ],
      "text/plain": [
       "<Quantity 3.6563664 arcsec>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radius_pixels * pixel_scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the radius we need to smear the source over, i.e., the radius in arcseconds we just calculated above, bigger or smaller than the expected source size (the \"point spread function\") given the typical seeing conditions at CTO ($\\sim2\"$)?\n",
    "\n",
    "Answer: This radius is greater than the expected source size. \n",
    "\n",
    "Given that answer, can we observe this target without saturating the CCD?  Justify your answer.\n",
    "\n",
    "Answer: No; since the radius is greater than the expected source size, we would need to be able to view a greater radius in order to avoid oversaturating the pixels. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration Times\n",
    "\n",
    "In the above, we were simply determining whether the source would be above or below the saturation level of our CCDs.  Let's now determine how long we need to integrate to achieve a target signal to noise ratio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We decide our integration time based on a target signal-to-noise ratio (SNR).  If our target SNR is 100, how long do we need to integrate?\n",
    "\n",
    "To a first approximation, we can assume that the signal-to-noise is limited by counting statistics.\n",
    "\n",
    "The signal $S$ rises as the count rate $\\dot{p}$ times the time,  $S = \\dot{p} t$.\n",
    "\n",
    "The noise $N$ is, in the Poisson noise limit, the square root of the signal (in photons), i.e., $N = \\sqrt{S} = \\sqrt{\\dot{p}t}$.  The SNR is therefore \n",
    "$$S/N = \\frac{\\dot{p} t}{\\sqrt{\\dot{p} t}} = \\sqrt{\\dot{p} t}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Warning / Note\n",
    "\n",
    "When we say that the uncertainty is the square root of the signal when counting photons, we have a bit of a units problem!\n",
    "\n",
    "$$ \\sigma = \\sqrt{\\lambda}$$\n",
    "\n",
    "The number of photons counted has units: number of photons.  The _uncertainty_ on that number has the same units, despite the fact that, naively, you would look at the above equation and say that the uncertainty has units of $\\sqrt{counts}$.\n",
    "\n",
    "You can't rely on the astropy unit system (or any automated unit system) to handle this equivalency; you need to take care to establish the correct units when dealing with uncertainties on counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculated the count rate above ($\\dot{p}$ = `phot_per_s`), so what integration time $t$ do we need for $SNR=100$ ?\n",
    "\n",
    "(solve the above equation, $SNR = \\sqrt{\\dot{p} t}$, for $t$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.00091036592 \\; \\mathrm{s}$"
      ],
      "text/plain": [
       "<Quantity 0.00091037 s>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer (in seconds)\n",
    "integration_time = (100**2)/phot_per_s\n",
    "integration_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should note that that is a pretty short time! (about 0.9 milliseconds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we account for read noise?  Let's assume that we are averaging over a 2 arcsecond Gaussian PSF.  First, how many pixels is that?\n",
    "\n",
    "The area of a Gaussian isn't exactly that of a circle.  $A_{2D \\mathrm{Gaussian}} = 2 \\pi \\sigma^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$25.132741 \\; \\mathrm{arcsec^{2}}$"
      ],
      "text/plain": [
       "<Quantity 25.13274123 arcsec2>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psf_area = 2 * np.pi * (2*u.arcsec)**2\n",
    "psf_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100.53096 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 100.53096491>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psf_area_pixels = psf_area * pixel_scale**-2\n",
    "psf_area_pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that our readnoise is 10 counts per pixel, and we're interested in the _sum_ over these pixels, what is the expected contribution to the noise on the sum from these pixels?\n",
    "\n",
    "Remember the propagation of error formula for summation:\n",
    "$$sum = \\Sigma x_i$$\n",
    "$$\\sigma_{sum}^2 = \\Sigma_i \\sigma_i^2$$\n",
    "$$\\sigma_{sum} = \\sqrt{\\Sigma_i \\sigma_i^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100.26513 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 100.26513099 adu>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer\n",
    "readnoise_per_pix = 10*u.adu/u.pix\n",
    "readnoise_sum = (psf_area_pixels*readnoise_per_pix**2)**(1/2)\n",
    "readnoise_sum*u.pix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By coincidence, that's about 100.265 ADU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given our calculation of the SNR and count noise above, what's our expected _total_ noise in the measurement over the PSF? (recall: the signal will be $t \\times \\dot{p}$, and we have calculated both of those quantities above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 100. adu>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer (note: this value should be in the same units as the read noise above - you might have to add those units!)\n",
    "signal = (integration_time*u.s)*(phot_per_s*(1/u.s))\n",
    "photon_noise = (signal)**(1/2)\n",
    "photon_noise*u.adu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$141.42136 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 141.42135624 adu>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer\n",
    "total_noise = (signal+(photon_noise)**2)**(1/2)\n",
    "total_noise*u.adu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Should be about sqrt(2) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that our noise is higher than we originally thought, we will not achieve our target signal to noise ratio if we use the `integration_time` above!  What will we achieve instead?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 100.>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# student answer\n",
    "signal = (integration_time*u.s)*(phot_per_s*(1/u.s))\n",
    "snr = signal**(1/2)\n",
    "snr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write down a different equation for the signal-to-noise ratio now that we know how much readnoise to expect for a stellar observation.\n",
    "\n",
    "We have just computed that our total noise is $N = {\\sqrt{\\dot{p} t + \\sigma_{rn}^2}}$, so our SNR is:\n",
    "\n",
    "$$S/N = \\frac{\\dot{p} t}{\\sqrt{\\dot{p} t + \\sigma_{rn}^2}} \\approx \\frac{\\dot{p} t}{\\sqrt{\\dot{p} t + 100}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This equation is no longer trivially solvable, but it can be rearranged nicely enough into a quadratic equation:\n",
    "$$SNR^2 \\left( \\dot{p} t + 100\\right) = \\left(\\dot{p}t\\right)^2$$ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "if we substitute $S = \\dot{p} t$, it looks a little simpler\n",
    "$$ S^2 - S \\cdot SNR^2  - 100 SNR^2 = 0$$\n",
    "then, we can use the quadratic formula to solve for the needed signal:\n",
    "$$S = SNR^2 \\pm \\frac{\\sqrt{SNR^4 + 400 SNR^2}}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15099.279489285953"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNR = 100\n",
    "target_signal = SNR**2 + (SNR**4 + 4*readnoise_sum.value*SNR**2)**0.5 / 2\n",
    "target_signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.0013745869 \\; \\mathrm{s}$"
      ],
      "text/plain": [
       "<Quantity 0.00137459 s>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integration_time_including_readnoise = target_signal / phot_per_s\n",
    "integration_time_including_readnoise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.5099279 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 1.50992795>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# how much longer is that?\n",
    "integration_time_including_readnoise / integration_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, accounting for readnoise, to get to the same SNR, we need to integrate for about twice as long (2.12x)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What would we need to be even more accurate in our prediction?\n",
    "\n",
    "Note that all of these calculations make a few assumptions that are generally not correct:\n",
    "\n",
    "1. We assumed the V-filter has a uniform transmission of 1 across its band.  In reality, the V-filter transmits closer to ~75% of the total band, so the total photon rate (in any given filter) is less than what we calculated above.\n",
    "\n",
    "2. We assumed that the CCD has perfect _quantum efficiency_.  The best CCDs have close to unity (one) efficiency, but not exactly.  Ours are probably in the 60-80% range.\n",
    "\n",
    "3. We assumed no flux losses to the atmosphere (absorption or scattering) or the telescope optics.  We know from class that the atmosphere absorbs anywhere from 10-30% of the incident flux if we point straight up on a clear night.\n",
    "\n",
    "Together these effects amount to a substantial loss of light.  Let's estimate how much:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.41999999999999993"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filter_efficiency = 0.75\n",
    "quantum_efficiency = 0.7\n",
    "atmosphere_loss = 0.2\n",
    "received_fraction = filter_efficiency * quantum_efficiency * (1-atmosphere_loss)\n",
    "received_fraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're getting less than 50% of the light recorded on our CCD.  That means we're liable to underestimate the brightness of our targets by a lot!  Next time, we'll go over how we account for that loss by _calibrating_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do it yourself for a more reasonable target.\n",
    "\n",
    "We noted above that P Cygni is really hard to observe, since it will saturate our detector.\n",
    "\n",
    "Now, let's repeat the above exercises for a different target: `TYC 2788-108-1`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer the following questions (using the code you learned above):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.  How bright is TYC 2788-108-1 in the [Johnson V-band](http://voservices.net/filter/filterdtls.aspx?filterid=60)?  Give your answer in magnitudes and in Janskys, assuming the zero-magnitude flux is 3500 Jy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Table length=1</i>\n",
       "<table id=\"table1568601241440\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>MAIN_ID</th><th>RA</th><th>DEC</th><th>RA_PREC</th><th>DEC_PREC</th><th>COO_ERR_MAJA</th><th>COO_ERR_MINA</th><th>COO_ERR_ANGLE</th><th>COO_QUAL</th><th>COO_WAVELENGTH</th><th>COO_BIBCODE</th><th>FLUX_R</th><th>FLUX_V</th><th>SCRIPT_NUMBER_ID</th></tr></thead>\n",
       "<thead><tr><th></th><th>&quot;h:m:s&quot;</th><th>&quot;d:m:s&quot;</th><th></th><th></th><th>mas</th><th>mas</th><th>deg</th><th></th><th></th><th></th><th>mag</th><th>mag</th><th></th></tr></thead>\n",
       "<thead><tr><th>object</th><th>str13</th><th>str13</th><th>int16</th><th>int16</th><th>float32</th><th>float32</th><th>int16</th><th>str1</th><th>str1</th><th>object</th><th>float32</th><th>float32</th><th>int32</th></tr></thead>\n",
       "<tr><td>TYC 2788-108-1</td><td>00 39 48.0059</td><td>+40 54 40.678</td><td>14</td><td>14</td><td>0.011</td><td>0.010</td><td>90</td><td>A</td><td>O</td><td>2020yCat.1350....0G</td><td>--</td><td>10.65</td><td>1</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Table length=1>\n",
       "   MAIN_ID           RA           DEC      ...  FLUX_R  FLUX_V SCRIPT_NUMBER_ID\n",
       "                  \"h:m:s\"       \"d:m:s\"    ...   mag     mag                   \n",
       "    object         str13         str13     ... float32 float32      int32      \n",
       "-------------- ------------- ------------- ... ------- ------- ----------------\n",
       "TYC 2788-108-1 00 39 48.0059 +40 54 40.678 ...      --   10.65                1"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Simbad.reset_votable_fields()\n",
    "Simbad.add_votable_fields('flux(R)', 'flux(V)')\n",
    "tyc_simbad = Simbad.query_object('TYC 2788-108-1')\n",
    "tyc_simbad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TYC 2788-108-1 is 10.65 magnitudes in the V-band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.19233931 \\; \\mathrm{Jy}$"
      ],
      "text/plain": [
       "<Quantity 0.19233931 Jy>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_tyc = 10**(-10.65/2.5) * 3500*u.Jy\n",
    "F_tyc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TYC 2788-108-1 is 0.19233931 Janskys in the V-band."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Using the same telescope as above, what is our expected received flux in ergs/second and in photons per second?  Assume the V-filter band width is 870 Angstroms and the band center is 5504 Angstroms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$5.4468106 \\times 10^{14} \\; \\mathrm{Hz}$"
      ],
      "text/plain": [
       "<Quantity 5.44681065e+14 Hz>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_CTO = (np.pi)*((14/2)*2.54*u.cm)**2\n",
    "wavelength = 5504*u.angstrom\n",
    "frequency = (constants.c  / wavelength).to(u.Hz)\n",
    "frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$8.6096026 \\times 10^{13} \\; \\mathrm{Hz}$"
      ],
      "text/plain": [
       "<Quantity 8.60960258e+13 Hz>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dlam = 870*u.angstrom\n",
    "lam = wavelength\n",
    "dnu = (dlam/lam)*((constants.c)/(5504e-10*u.m))\n",
    "dnu.to(u.Hz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.6446161 \\times 10^{-7} \\; \\mathrm{\\frac{erg}{s}}$"
      ],
      "text/plain": [
       "<Quantity 1.64461609e-07 erg / s>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_rec = F_tyc*A_CTO*dnu\n",
    "E_rec.to(u.erg/u.s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$53814.729 \\; \\mathrm{\\frac{1}{s}}$"
      ],
      "text/plain": [
       "<Quantity 53814.72925875 1 / s>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phot_per_s = (E_rec.to(u.erg/u.s))/(E_phot.to(u.erg))\n",
    "phot_per_s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. How many photons per second per pixel will the star produce if the point spread function, the area over which the star's light is spread, is the same as above ($\\sigma_{Gaussian} = 2\"$)?  How long can we integrate without saturating (assuming we saturate at 30,000 counts)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$4 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 4.>"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2*u.arcsec)/pixel_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50.26548245743669"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pixel_area = np.pi*4**2\n",
    "pixel_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1070.61 \\; \\mathrm{\\frac{1}{s}}$"
      ],
      "text/plain": [
       "<Quantity 1070.61002159 1 / s>"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phot_per_s/pixel_area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Over this area, there will be 1070.61 photons per second per pixel produced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.5574682 \\; \\mathrm{s}$"
      ],
      "text/plain": [
       "<Quantity 0.5574682 s>"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#30000 counts = photons per second * time\n",
    "t = (30000)/(phot_per_s)\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can integrate for approximately 0.56 seconds without saturating. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Again following the same assumptions (the same spread among pixels, same readnoise), how long do we need to integrate to get to a signal-to-noise ratio of 100?  Account for both photon statistics and readnoise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.18582273 \\; \\mathrm{s}$"
      ],
      "text/plain": [
       "<Quantity 0.18582273 s>"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integration_time = (100**2)/phot_per_s\n",
    "integration_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100.53096 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 100.53096491>"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psf_area = 2 * np.pi * (2*u.arcsec)**2\n",
    "psf_area_pixels = psf_area * pixel_scale**-2\n",
    "psf_area_pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100.26513 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 100.26513099 adu>"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "readnoise_per_pix = 10*u.adu/u.pix\n",
    "readnoise_sum = (psf_area_pixels*readnoise_per_pix**2)**(1/2)\n",
    "readnoise_sum*u.pix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 100. adu>"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "signal = (integration_time*u.s)*(phot_per_s*(1/u.s))\n",
    "photon_noise = (signal)**(1/2)\n",
    "photon_noise*u.adu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$141.42136 \\; \\mathrm{adu}$"
      ],
      "text/plain": [
       "<Quantity 141.42135624 adu>"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_noise = (signal+(photon_noise)**2)**(1/2)\n",
    "total_noise*u.adu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$100 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 100.>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "signal = (integration_time*u.s)*(phot_per_s*(1/u.s))\n",
    "snr = signal**(1/2)\n",
    "snr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15099.279489285953"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SNR = 100\n",
    "target_signal = SNR**2 + (SNR**4 + 4*readnoise_sum.value*SNR**2)**0.5 / 2\n",
    "target_signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.28057894 \\; \\mathrm{s}$"
      ],
      "text/plain": [
       "<Quantity 0.28057894 s>"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integration_time_including_readnoise = target_signal / phot_per_s\n",
    "integration_time_including_readnoise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to integrate for 0.28 seconds to get a SNR of 100. "
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}