8GHz_j1427-4206_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radio flux density estimation of J1427-4206 @ 8GHz\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Please note:\n",
    "\n",
    "The data you will be processing in this and other tutorials is located in the \"data\" folder. In there you will\n",
    "find fits files for your calibrator [Hydra A](http://simbad.u-strasbg.fr/simbad/sim-id?Ident=NAME%20HYDRA%20A)\n",
    "and your target source [J1427-4206](http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=PKS+1424-418&submit=SIMBAD+search) at 2 GHz, 5 GHz, 8 GHz and 12 GHz frequencies. **For this exercize we will be working on\n",
    "8GHz data for our target source**.\n",
    "\n",
    "All the files you will need for the tutorial can be found at this [Github repository](https://github.com/Pfesi/AVN2019). **You need to download the file and work from that directory**. i.e. \"cd into that directory\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Listed below are the steps we are going to follow to calibrate our source and estimate the flux density\n",
    "\n",
    "# <a name=\"top\"></a>Task list:\n",
    "\n",
    "1. [Locate and open the observation file](#section1)\n",
    "2. [Examine the fits file](#section2)\n",
    "3. [Extract the drift scan data](#section3)\n",
    "4. [Measure the antenna temperature](#section4)\n",
    " * 4.1 [Convert counts to units of Kelvin](#section4.1)\n",
    " * 4.2 [Fit a baseline to remove systemic contributions](#section4.2)\n",
    " * 4.3 [Fit beam to get antenna temperature](#section4.3)\n",
    " * 4.4 [Amplitude correction](#section4.4)\n",
    "5. [Calculate the flux density](#section5)\n",
    "6. [Got to home page](AVN_2018_Drift_scan_continuum_tutorial.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 1. <a name=\"section1\">Locate and open the observation file</a> \n",
    "[back to top](#top)\n",
    "\n",
    "Select an observation to work on and open the fits file\n",
    "\n",
    "* The Python standard library has a [pyfits package](http://docs.astropy.org/en/stable/io/fits/index.html) we use for reading and writing FITS files and manipulating their contents. Fits files store information about a source observation in what it calls header data units (HDUs). We are going to load the package in order to begin processing our preferred file. You can read up on all the other interesting file manipulation methods by clicking on the link above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.io.fits as pyfits\n",
    "\n",
    "# Get the path name of the file you want to process\n",
    "fitsfile = 'data/J1427-4206/8GHz/2013d125_21h41m04s_Cont_mike_J1427-4206.fits'\n",
    "hdulist  = pyfits.open(fitsfile) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "source": [
    "# 2. <a name=\"section2\">Examine the fits file</a>\n",
    "[back to top](#top)\n",
    "\n",
    "It is always a good idea to familiarize yourself with the file contents by viewing the different types of information inside your file. To get an overview of the contents we use the info() method\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The info() method gives us a list of all the HDUs in this file. \n",
    "\n",
    "* The first header data unit (HDU) contains information on the  observation.  \n",
    "* The second one has information on the 2.5 cm (i.e. 12 GHz) feed system. This data is not necessarily up to date, so we generally disregard it. \n",
    "* The third unit has the noise diode firing, used to convert from raw counts to Kelvins. \n",
    "* The next three binary tables are the drift scans themselves, starting at the north offset position.\n",
    "* The last HDU contains more information about the observation.\n",
    "\n",
    "To view the content inside an HDU we use indexes. For example, to view the PRIMARY HDU, we use index 0.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist[0].header\n",
    "\n",
    "# You can play around with changing the indeces to see what \n",
    "# other information is stored in the other HDUs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. <a name=\"section3\">Extract the drift scan data</a>\n",
    "\n",
    "[back to top](#top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drift1 = hdulist[3] #HPN\n",
    "drift2 = hdulist[4] #ON\n",
    "drift3 = hdulist[5] #HPS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets view the drift1 HDU to see the file contents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drift1.header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the fields are actually just the position in different formats. We are interested in the output signal '**Count1**' and '**Count2**'.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. <a name=\"section4\">Measure the antenna temperature</a>\n",
    "[back to top](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 <a name=\"section4.1\">Convert counts to units of Kelvin </a>\n",
    "[back to top](#top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import plotting libraries\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#display figures inline \n",
    "%matplotlib inline  \n",
    "\n",
    "# Get the noise diode calibration HDU\n",
    "noise_cal = hdulist[2]\n",
    "\n",
    "# Construct an array for the x-axis in terms of right ascension\n",
    "ra = drift1.data['RA_J2000']\n",
    "\n",
    "count_K1 = noise_cal.header['HZPERK1']\n",
    "count_K2 = noise_cal.header['HZPERK2']\n",
    "\n",
    "# Get source counts and convert counts to antenna temperature\n",
    "north_scan_LCP = (drift1.data['Count1']/count_K1) - ((drift1.data['Count1'])[0]/count_K1) \n",
    "north_scan_RCP = (drift1.data['Count2']/count_K2) - ((drift1.data['Count2'])[0]/count_K2)\n",
    "on_scan_LCP    = (drift2.data['Count1']/count_K1) - ((drift2.data['Count1'])[0]/count_K1)\n",
    "on_scan_RCP    = (drift2.data['Count2']/count_K2) - ((drift2.data['Count2'])[0]/count_K2)\n",
    "south_scan_LCP = (drift3.data['Count1']/count_K1) - ((drift3.data['Count1'])[0]/count_K1)\n",
    "south_scan_RCP = (drift3.data['Count2']/count_K2) - ((drift3.data['Count2'])[0]/count_K2)\n",
    "            \n",
    "# Plot the data\n",
    "\n",
    "# HPN \n",
    "plt.figure(figsize=[15,5])\n",
    "plt.suptitle('North scan')\n",
    "\n",
    "ax = plt.subplot(121)\n",
    "plt.grid()\n",
    "plt.plot(ra, north_scan_LCP, label = 'LCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend()\n",
    "lim = plt.axis('tight')\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "plt.grid()\n",
    "plt.plot(ra, north_scan_RCP, label = 'RCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend()\n",
    "lim = plt.axis('tight')\n",
    "\n",
    "# ON\n",
    "plt.figure(figsize=[15,5])\n",
    "plt.suptitle('Central (on-source) scan')\n",
    "\n",
    "ax = plt.subplot(121)\n",
    "plt.grid()\n",
    "plt.plot(ra, on_scan_LCP, label = 'LCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend()\n",
    "lim = plt.axis('tight')\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "plt.grid()\n",
    "plt.plot(ra, on_scan_RCP, label = 'RCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend()\n",
    "lim = plt.axis('tight')\n",
    "\n",
    "# HPS\n",
    "plt.figure(figsize=[15,5])\n",
    "plt.suptitle('South scan')\n",
    "\n",
    "ax = plt.subplot(121)\n",
    "plt.grid()\n",
    "plt.plot(ra, south_scan_LCP, label = 'LCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend()\n",
    "lim = plt.axis('tight')\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "plt.plot(ra, south_scan_RCP, label = 'RCP')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Ta (K)')\n",
    "plt.legend(frameon=False)\n",
    "plt.grid()\n",
    "lim = plt.axis('tight')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 <a name=\"section4.2\">Fit a baseline to remove systemic contributions</a>\n",
    "[back to top](#top)\n",
    "\n",
    "We are only interested in the temperature contributed by the source, so we want to subract off the system temperature contribution.  We see that the baseline level changes through during the scan.\n",
    "\n",
    "What do we fit then?  This is a fairly strong source, we can clealy see the antenna beam pattern.  Looking carefully we can see a deflection point at the base of the beam.  This is the first null of the beam. At these points we will not be receiving any power from the source.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Locate the section you want to use for your fit. We are going to write a few methods to help us process the location selection for our fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getBasePts(x, y,len_scan,x1,x2,x3,x4): \n",
    "    '''\n",
    "        Get baseline points. Select points along the driftscan where you will fit\n",
    "        your baseline.\n",
    "    '''\n",
    "\n",
    "    xbleft = x[x1:x2]\n",
    "    ybleft = y[x1:x2]\n",
    "    xbright= x[len_scan-x3:len_scan-x4]\n",
    "    ybright= y[len_scan-x3:len_scan-x4]\n",
    "\n",
    "    left_base_pts = list(xbleft) + list(xbright)\n",
    "    right_base_pts = list(ybleft) + list(ybright)\n",
    "\n",
    "    return left_base_pts, right_base_pts \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the baseline and plot your data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def fitBasePts(x_base_pts,y_base_pts, x, y):\n",
    "    \"\"\"Fit the data to remove systematic contributions.\n",
    "    \"\"\"\n",
    "        \n",
    "    base_fit_coeffs = np.polyfit(x_base_pts, y_base_pts, 1)\n",
    "    base_fit_line = np.polyval(base_fit_coeffs,x_base_pts)   \n",
    "    data_fit = np.polyval(base_fit_coeffs, x)\n",
    "    data_fit = y - data_fit\n",
    "    \n",
    "    res,rms = residual(y_base_pts,base_fit_line)\n",
    "   \n",
    "    print (\"Fit = %.2fx + %d, rms error = %.3f\" %(base_fit_coeffs[0], base_fit_coeffs[1],rms))\n",
    "    \n",
    "    return data_fit, base_fit_line\n",
    "    \n",
    "# Get the residual and rms to estimate the errors in the fit\n",
    "def residual(model, data): \n",
    "    \n",
    "    from sklearn.metrics import mean_squared_error\n",
    "    \"\"\"\n",
    "        Calculate the residual between the model and the data.\n",
    "    \"\"\"\n",
    "\n",
    "    res = np.array(model - data)\n",
    "    rms = np.sqrt(mean_squared_error(data,model))\n",
    "\n",
    "    return res, rms\n",
    "\n",
    "def fitBaselineAndPlot(x,y,title,x1,x2,x3,x4,len_scan):\n",
    "    \"\"\"\n",
    "        Fit the baseline and plot your data\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    xb,yb = getBasePts(x,y,len_scan,x1,x2,x3,x4)\n",
    "    \n",
    "    # Fit the baseline points\n",
    "    fit, base_fit_line = fitBasePts(xb,yb, x, y)\n",
    "\n",
    "    # plot the data\n",
    "    plt.figure(figsize=[15,5])\n",
    "    plt.suptitle(title)\n",
    "\n",
    "    ax = plt.subplot(121)\n",
    "    plt.axhline(y=0, color='k', alpha= 0.5)\n",
    "    l1,= plt.plot(x,y, label = 'original data')\n",
    "    l2,= plt.plot(xb,yb,'r.', label = 'noise')\n",
    "    plt.plot(xb,base_fit_line, 'g', label = \"baseline fit\")\n",
    "    plt.xlabel('RA (deg)')\n",
    "    plt.ylabel('Ta (K)')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    lim = plt.axis('tight')\n",
    "\n",
    "    ax = plt.subplot(122)\n",
    "    l1,= plt.plot(x,fit, label = 'final data')\n",
    "    plt.xlabel('RA (deg)')\n",
    "    plt.ylabel('Ta (K)')\n",
    "    plt.axhline(y=0, color='k', alpha= 0.5)\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    lim = plt.axis('tight')\n",
    "    \n",
    "    return fit\n",
    "\n",
    "# Get the length of the scans\n",
    "len_scan = len(ra)\n",
    "\n",
    "# Set the baseline points, default is 100 points on each side\n",
    "# You can change the x1-x4 values to adjust the location of\n",
    "# your baseline fit.\n",
    "\n",
    " # YOU CAN MODIFY THESE VALUES\n",
    "x1 = \n",
    "x2 = \n",
    "x3 = \n",
    "x4 = \n",
    "\n",
    "print (\"Y = mx + c\")\n",
    "# Fit and plot \n",
    "\n",
    "# *NB* If you want to change individual baseline points, do it here.\n",
    "slcp  = fitBaselineAndPlot(ra, south_scan_LCP, 'South scan [LCP]',x1,x2,x3,x4,len_scan)\n",
    "srcp  = fitBaselineAndPlot(ra, south_scan_RCP, 'South scan [RCP]',x1,x2,x3,x4,len_scan)\n",
    "nlcp  = fitBaselineAndPlot(ra, north_scan_LCP, 'North scan [LCP]',x1,x2,x3,x4,len_scan)\n",
    "nrcp  = fitBaselineAndPlot(ra, north_scan_RCP, 'Northscan [RCP]',x1,x2,x3,x4,len_scan)\n",
    "onlcp = fitBaselineAndPlot(ra, on_scan_LCP, 'On scan [LCP]',x1,x2,x3,x4,len_scan)\n",
    "onrcp = fitBaselineAndPlot(ra, on_scan_RCP, 'On scan [RCP]',x1,x2,x3,x4,len_scan)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 <a name=\"section4.3\">Fit beam to get antenna temperature</a>\n",
    "[back to top](#top)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now going to fit the top of the beam to fit the \n",
    "antenna temperature. We will use a 2nd order polynomial fit\n",
    "at the top of the peak. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get location of peak\n",
    "def getAPeakPts(drift, percentage):\n",
    "    \"\"\"Get locations of scan peak where we are going to fit \n",
    "       We will be fitting the top 40% of our scan.\n",
    "    \"\"\"\n",
    "    peak_frac = float((100.0 - percentage)/100.0)\n",
    "    peak_max = max(drift)\n",
    "    peak_pts = np.where(drift > peak_frac * peak_max)[0]\n",
    "    return peak_pts\n",
    "\n",
    "def getBPeakPts(drift, percentage):\n",
    "    \"\"\"Get locations of scan peak where we are going to fit \n",
    "       We will be fitting the top 40% of our scan.\n",
    "    \"\"\"\n",
    "    peak_frac = float((100.0 - percentage)/100.0)\n",
    "    peak_min = min(drift)\n",
    "    peak_pts = np.where(drift < peak_frac * peak_min)[0]\n",
    "    return peak_pts\n",
    "\n",
    "# Fit the peak\n",
    "def fitPeak(peak_pts, x, y):\n",
    "    \"\"\"Fit the peak and estimate errors.\n",
    "    \"\"\"\n",
    "\n",
    "    peakfitcoeffs = np.polyfit(x[peak_pts], y[peak_pts], 2)\n",
    "    peakfitline = np.polyval(peakfitcoeffs,x[peak_pts])\n",
    "    res, rms = residual(peakfitline, y[peak_pts])\n",
    "\n",
    "    return peakfitline,rms, x[peak_pts], y[peak_pts]\n",
    "\n",
    "# Get peak points\n",
    "\n",
    "# YOU CAN MODIFY THIS VALUE\n",
    "peak_percentage = 40 # Change the percentage of the peak data you want to fit.\n",
    "\n",
    "# Get peak points, Abeam\n",
    "Apeak_pts_slcp  = getAPeakPts(slcp, peak_percentage)\n",
    "Apeak_pts_srcp  = getAPeakPts(srcp, peak_percentage)\n",
    "Apeak_pts_nlcp  = getAPeakPts(nlcp, peak_percentage)\n",
    "Apeak_pts_nrcp  = getAPeakPts(nrcp, peak_percentage)\n",
    "Apeak_pts_onlcp = getAPeakPts(onlcp, peak_percentage)\n",
    "Apeak_pts_onrcp = getAPeakPts(onrcp, peak_percentage)\n",
    "\n",
    "# Get peak points, Bbeam\n",
    "Bpeak_pts_slcp  = getBPeakPts(slcp, peak_percentage)\n",
    "Bpeak_pts_srcp  = getBPeakPts(srcp, peak_percentage)\n",
    "Bpeak_pts_nlcp  = getBPeakPts(nlcp, peak_percentage)\n",
    "Bpeak_pts_nrcp  = getBPeakPts(nrcp, peak_percentage)\n",
    "Bpeak_pts_onlcp = getBPeakPts(onlcp, peak_percentage)\n",
    "Bpeak_pts_onrcp = getBPeakPts(onrcp, peak_percentage)\n",
    "\n",
    "# Fit the peaks\n",
    "Afitpsl, Aslrms, Axpslcp, Aypslcp   = fitPeak(Apeak_pts_slcp, ra, slcp)\n",
    "Afitpsr, Asrrms, Axpsrcp, Aypsrcp   = fitPeak(Apeak_pts_srcp, ra, srcp)\n",
    "Afitpnl, Anlrms, Axpnlcp, Aypnlcp   = fitPeak(Apeak_pts_nlcp, ra, nlcp)\n",
    "Afitpnr, Anrrms, Axpnrcp, Aypnrcp   = fitPeak(Apeak_pts_nrcp, ra, nrcp)\n",
    "Afitponl, Aonlrms, Axpolcp, Aypolcp = fitPeak(Apeak_pts_onlcp, ra, onlcp)\n",
    "Afitponr, Aonrrms, Axporcp, Ayporcp = fitPeak(Apeak_pts_onrcp, ra, onrcp)\n",
    "\n",
    "Bfitpsl, Bslrms, Bxpslcp, Bypslcp   = fitPeak(Bpeak_pts_slcp, ra, slcp)\n",
    "Bfitpsr, Bsrrms, Bxpsrcp, Bypsrcp   = fitPeak(Bpeak_pts_srcp, ra, srcp)\n",
    "Bfitpnl, Bnlrms, Bxpnlcp, Bypnlcp   = fitPeak(Bpeak_pts_nlcp, ra, nlcp)\n",
    "Bfitpnr, Bnrrms, Bxpnrcp, Bypnrcp   = fitPeak(Bpeak_pts_nrcp, ra, nrcp)\n",
    "Bfitponl, Bonlrms, Bxpolcp, Bypolcp = fitPeak(Bpeak_pts_onlcp, ra, onlcp)\n",
    "Bfitponr, Bonrrms, Bxporcp, Byporcp = fitPeak(Bpeak_pts_onrcp, ra, onrcp)\n",
    "\n",
    "# Print out the ANTENNA TEMPERATURES\n",
    "print ('Ta_ASLCP = %.3f +- %.3f [K], Ta_BSLCP = %.3f +- %.3f [K]' %(max(Afitpsl), Aslrms, min(Bfitpsl), Bslrms))\n",
    "print ('Ta_ASRCP = %.3f +- %.3f [K], Ta_BSRCP = %.3f +- %.3f [K]' %(max(Afitpsr), Asrrms, min(Bfitpsr), Bsrrms))\n",
    "print ('Ta_ANLCP = %.3f +- %.3f [K], Ta_BNLCP = %.3f +- %.3f [K]' %(max(Afitpnl), Anlrms, min(Bfitpnl), Bnlrms))\n",
    "print ('Ta_ANRCP = %.3f +- %.3f K], Ta_BNRCP = %.3f +- %.3f K]' %(max(Afitpnr), Anrrms, min(Bfitpnr), Bnrrms))\n",
    "print ('Ta_AONLCP = %.3f +- %.3f [K], Ta_BONLCP = %.3f +- %.3f [K]' %(max(Afitponl), Aonlrms, min(Bfitponl), Bonlrms))\n",
    "print ('Ta_AONRCP = %.3f +- %.3f [K], Ta_BONRCP = %.3f +- %.3f [K]' %(max(Afitponr), Aonrrms, min(Bfitponr), Bonrrms))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotFit(ra, drift1, title, peak_pts1,fit1,drift2, peak_pts2,fit2, peak_pts3,fit3, peak_pts4,fit4, lab1,lab2,lab3,lab4,Axp1, Ayp1,Bxp1, Byp1,Axp2, Ayp2,Bxp2, Byp2):\n",
    "    \n",
    "    # plot the data\n",
    "    plt.figure(figsize=[15,5])\n",
    "    plt.suptitle(title)\n",
    "    \n",
    "    ax = plt.subplot(121)\n",
    "    l1,= plt.plot(ra,drift1, label = 'data')\n",
    "    plt.plot(Axp1,Ayp1,'y')\n",
    "    plt.plot(Bxp1,Byp1,'k')\n",
    "    l3,= plt.plot(ra[peak_pts1],fit1,'r', label=lab1)\n",
    "    l4,= plt.plot(ra[peak_pts2],fit2,'m', label=lab2)    \n",
    "    plt.xlabel('RA (deg)')\n",
    "    plt.ylabel('Ta (K)')\n",
    "    plt.legend()\n",
    "    plt.axhline(y=0, color='k', alpha= 0.5)\n",
    "    plt.grid()\n",
    "    lim = plt.axis('tight')\n",
    "\n",
    "    ax = plt.subplot(122)\n",
    "    l1,= plt.plot(ra,drift2, label = 'data')\n",
    "    plt.plot(Axp2,Ayp2,'y')\n",
    "    plt.plot(Bxp2,Byp2,'k')\n",
    "    l3,= plt.plot(ra[peak_pts3],fit3, 'r',label = lab3)\n",
    "    l4,= plt.plot(ra[peak_pts4],fit4, 'm', label = lab4)    \n",
    "    plt.xlabel('RA (deg)')\n",
    "    plt.ylabel('Ta (K)')\n",
    "    plt.axhline(y=0, color='k', alpha= 0.5)\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    lim = plt.axis('tight')\n",
    "\n",
    "plotFit(ra,slcp,'HPS PEAK FIT',Apeak_pts_slcp, Afitpsl,  srcp, Bpeak_pts_slcp, Bfitpsl, Apeak_pts_srcp, Afitpsr, Bpeak_pts_srcp, Bfitpsr, 'Ta_ASLCP  = %.3f +- %.3f [K]' %(max(Afitpsl), Aslrms),  'Ta_BSLCP = %.3f +- %.3f [K]' %(min(Bfitpsl), Bslrms),  'Ta_ASRCP = %.3f +- %.3f [K]' %(max(Afitpsr), Asrrms), 'Ta_BSRCP = %.3f +- %.3f [K]' %(min(Bfitpsr), Bsrrms),  Axpslcp, Aypslcp,Bxpslcp, Bypslcp,Axpsrcp, Aypsrcp,Bxpsrcp, Bypsrcp)\n",
    "plotFit(ra,nlcp,'HPN PEAK FIT',Apeak_pts_nlcp, Afitpnl,  nrcp, Bpeak_pts_nlcp, Bfitpnl, Apeak_pts_nrcp, Afitpnr, Bpeak_pts_nrcp, Bfitpnr, 'Ta_ANLCP  = %.3f +- %.3f [K]' %(max(Afitpnl), Anlrms),  'Ta_BNLCP = %.3f +- %.3f [K]' %(min(Bfitpnl), Bnlrms),  'Ta_ANRCP = %.3f [+- %.3f K]' %(max(Afitpnr), Anrrms), 'Ta_BNRCP = %.3f [+- %.3f K]' %(min(Bfitpnr), Bnrrms),  Axpnlcp, Aypnlcp,Bxpnlcp, Bypnlcp,Axpnrcp, Aypnrcp,Bxpnrcp, Bypnrcp)\n",
    "plotFit(ra,onlcp,'ON PEAK FIT',Apeak_pts_onlcp,Afitponl, onrcp,Bpeak_pts_onlcp,Bfitponl,Apeak_pts_onrcp,Afitponr,Bpeak_pts_onrcp,Bfitponr,'Ta_AONLCP = %.3f +- %.3f [K]' %(max(Afitponl),Aonlrms), 'Ta_BONLCP = %.3f +- %.3f [K]'%(min(Bfitponl),Bonlrms),'Ta_AONRCP = %.3f +- %.3f [K]' %(max(Afitponr),Aonrrms),'Ta_BONRCP = %.3f +- %.3f [K]'%(min(Bfitponr),Bonrrms), Axpolcp,Aypolcp,Bxpolcp, Bypolcp,Axporcp, Ayporcp,Bxporcp, Byporcp)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 <a name=\"section4.4\">Amplitude correction</a>\n",
    "[back to top](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see our pointing is off and we need to correct for that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correctAmp(hpn,errhpn,hps,errhps,on,erron):\n",
    "\n",
    "    ''' Calculate the pointing correction\n",
    "\n",
    "        Tcorr = Ton x exp[((ln(Thps) - ln(Thpn))**2)/(16*ln(2))]\n",
    "\n",
    "        Tcorr = Ton x pc\n",
    "\n",
    "        term1 = ((ln(Thps) - ln(Thpn))**2)\n",
    "        term2 = 16*ln(2)\n",
    "        term3 = term1/term2\n",
    "        pc = exp[term3]\n",
    "    '''\n",
    "\n",
    "    # Calculate the pointing correction\n",
    "    term1 = (np.log(abs(hps)) - np.log(abs(hpn)))**2\n",
    "    term2 = 16*np.log(2)\n",
    "    term3 = term1/term2\n",
    "    pc = np.exp(term3)\n",
    "            \n",
    "    # Calculate the corrected antenna temperature\n",
    "    Tcorr = abs(on * pc)\n",
    "\n",
    "    # Calculate the error in the pointing correction\n",
    "    if hps < 0.0:\n",
    "        pf1 = pc * 2.0 * term3 * (1.0/hps)\n",
    "    else:\n",
    "        pf1 = pc * 2.0 * term3 * (1.0/abs(hps))\n",
    "                                \n",
    "    if hpn < 0.0:\n",
    "        pf2 = pc * 2.0 * term3 * (-1.0/hpn)\n",
    "    else:\n",
    "        pf2 = pc * 2.0 * term3 * (1.0/abs(hpn))\n",
    "            \n",
    "    pcerr = np.sqrt((errhpn**2) + (errhps**2)) # err in pc\n",
    "\n",
    "    # Calculate the error in the corrected antenna temperature\n",
    "    Tcorrerr = np.sqrt((erron/on)**2 + (pcerr/pc)**2)*Tcorr\n",
    "            \n",
    "    return pc, pcerr, Tcorr, Tcorrerr\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the corrected amplitude\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Pointing correction\n",
    "Apcl, Apcerrl, ATcorrl, ATcorerrl = correctAmp(max(Afitpnl), Anlrms, max(Afitpsl), Aslrms, max(Afitponl), Aonlrms)\n",
    "Bpcl, Bpcerrl, BTcorrl, BTcorerrl = correctAmp(min(Bfitpnl), Bnlrms, min(Bfitpsl), Bslrms, min(Bfitponl), Bonlrms)\n",
    "\n",
    "Apcr, Apcerrr, ATcorrr, ATcorerrr = correctAmp(max(Afitpnr), Anrrms, max(Afitpsr), Asrrms, max(Afitponr), Aonrrms)\n",
    "Bpcr, Bpcerrr, BTcorrr, BTcorerrr = correctAmp(min(Bfitpnr), Bnrrms, min(Bfitpsr), Bsrrms, min(Bfitponr), Bonrrms)\n",
    "\n",
    "print ('Ta_ALCP corr = %.4f +- %.4f [K]' %(max(Afitponl), Aonlrms))\n",
    "print ('Ta_BLCP corr = %.4f +- %.4f [K]' %(min(Bfitponl), Bonlrms))\n",
    "print ('Ta_ARCP corr = %.4f +- %.4f [K]' %(max(Afitponr), Aonrrms)) \n",
    "print ('Ta_BRCP corr = %.4f +- %.4f [K]\\n' %(min(Bfitponr), Bonrrms))\n",
    "\n",
    "print ('Amplitude correction factor due to pointing error in A BEAM LCP = %.4f +- %f [K] ' %(Apcl, Apcerrl))\n",
    "print ('Amplitude correction factor due to pointing error in B BEAM LCP = %.4f +- %f [K] ' %(Bpcl, Bpcerrl))\n",
    "print ('Amplitude correction factor due to pointing error in A BEAM RCP = %.4f +- %f [K] ' %(Apcr, Apcerrr))\n",
    "print ('Amplitude correction factor due to pointing error in B BEAM RCP = %.4f +- %f [K]\\n' %(Bpcr, Bpcerrr))\n",
    "\n",
    "print ('Ta_ALCP corr = %.4f +- %.4f [K]' %(ATcorrl, ATcorerrl))\n",
    "print ('Ta_BLCP corr = %.4f +- %.4f [K]' %(BTcorrl, BTcorerrl))\n",
    "print ('Ta_ARCP corr = %.4f +- %.4f [K]' %(ATcorrr, ATcorerrr))\n",
    "print ('Ta_BRCP corr = %.4f +- %.4f [K]' %(BTcorrr, BTcorerrr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. <a name=\"section5\">Calculate the Flux density</a>\n",
    "[back to top](#top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fluxDensity(pss,errpss,ta,errta):\n",
    "\n",
    "    '''Calculate the source flux density'''\n",
    "    \n",
    "    flux = pss*ta\n",
    "    errflux = np.sqrt((errpss/pss)**2 + (errta/ta)**2) * flux\n",
    "\n",
    "    return flux, errflux\n",
    "\n",
    "def totfluxDensity(Afluxlcp,Aerrfluxlcp,Bfluxlcp,Berrfluxlcp,Afluxrcp,Aerrfluxrcp,Bfluxrcp,Berrfluxrcp):\n",
    "\n",
    "    '''Calculate the source flux density'''\n",
    "    \n",
    "    flux = ((Afluxlcp + Bfluxlcp + Afluxrcp + Bfluxrcp)/4.0)*2.0\n",
    "    errflux = np.sqrt((Aerrfluxlcp)**2 + (Berrfluxlcp)**2 + (Aerrfluxrcp)**2 + (Berrfluxrcp)**2 )\n",
    "\n",
    "    return flux, errflux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate flux density\n",
    "\n",
    "#PSS @ 5 GHz\n",
    "\n",
    "Apsslcp = \n",
    "Aerrpsslcp = \n",
    "Bpsslcp = \n",
    "Berrpsslcp = \n",
    "Apssrcp = \n",
    "Aerrpssrcp =  \n",
    "Bpssrcp = \n",
    "Berrpssrcp = \n",
    "\n",
    "\n",
    "#Flux\n",
    "Afluxlcp, Aerrfluxlcp = fluxDensity(Apsslcp,Aerrpsslcp,ATcorrl,ATcorerrl)\n",
    "Bfluxlcp, Berrfluxlcp = fluxDensity(Bpsslcp,Berrpsslcp,BTcorrl,BTcorerrl)\n",
    "Afluxrcp, Aerrfluxrcp = fluxDensity(Apssrcp,Aerrpssrcp,ATcorrr,ATcorerrr)\n",
    "Bfluxrcp, Berrfluxrcp = fluxDensity(Bpssrcp,Berrpssrcp,BTcorrr,BTcorerrr)\n",
    "\n",
    "\n",
    "totflux, errtotflux = totfluxDensity(Afluxlcp, Aerrfluxlcp,Bfluxlcp, Berrfluxlcp,Afluxrcp, Aerrfluxrcp,Bfluxrcp, Berrfluxrcp)\n",
    "\n",
    "print ('AFlux_LCP corr = %.3f +- %.3f [K]' %(Afluxlcp, Aerrfluxlcp))\n",
    "print ('BFlux_LCP corr = %.3f +- %.3f [K]' %(Bfluxlcp, Berrfluxlcp))\n",
    "print ('AFlux_RCP corr = %.3f +- %.3f [K]' %(Afluxrcp, Aerrfluxrcp))\n",
    "print ('BFlux_RCP corr = %.3f +- %.3f [K]' %(Bfluxrcp, Berrfluxrcp))\n",
    "print ('Total_flux = %.3f +- %.3f [K]' %(totflux, errtotflux))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Link to [AVN](background_info.ipynb) notebook "
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}