dva_noisesrc_gain_correction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e4659258",
   "metadata": {},
   "source": [
    "## Apply noise-diode-based gain corrections to azimuth and raster scans\n",
    "### A. Ordog, July 2022\n",
    "### NOTE: please do not modify\n",
    "#### July 12, 2022: updated to include RR gains, set outlier threshold based on stdev, write out new files with all 4 products corrected\n",
    "#### July 18, 2022: modified to use polynomial fit to gain vs time rather than interpolation\n",
    "#### Oct 5, 2022: modified to record corrected scans in noise source units\n",
    "#### Oct 26, 2022: changed noise source outlier thresholding, added possibility of multiple evening/morning raster scans\n",
    "#### Nov 22, 2022: corrected 'off' source calculation to not just use endpoints on either side of interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "07766661",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imp\n",
    "import os\n",
    "import subprocess\n",
    "import h5py\n",
    "import numpy as np\n",
    "from astropy.time import Time\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import griddata\n",
    "import datetime\n",
    "import matplotlib.dates as mdates\n",
    "from matplotlib.dates import HourLocator as HourLocator\n",
    "from matplotlib.dates import MinuteLocator as MinuteLocator\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from astropy import units as u\n",
    "from astropy.time import TimeDelta\n",
    "from scipy import interpolate\n",
    "from operator import itemgetter\n",
    "from itertools import groupby\n",
    "from astropy.convolution import convolve, Box1DKernel\n",
    "import gc\n",
    "\n",
    "phase = '1'\n",
    "day ='21'\n",
    "trim_rast_exist = False\n",
    "trim_az_exist = True\n",
    "Nm = 1 # number of morning raster scans\n",
    "Ne = 1 # number of evening raster scans\n",
    "\n",
    "#### Change the directory to where the files are located\" ####\n",
    "#dir_in = '/media/ordoga/DVA_data/survey_phase1_day'+day+'/'\n",
    "#dir_out = '/media/ordoga/DVA_data/survey_phase1_day'+day+'/'\n",
    "##############################################################\n",
    "\n",
    "### Use these directories on elephant: ###################\n",
    "dir_in_rast  = '/srv/data/dva/survey_raster/'\n",
    "dir_out_rast = '/srv/data/dva/survey_raster_noise_corr/'\n",
    "dir_in_az    = '/srv/data/dva/survey_azimuth_scans/'\n",
    "dir_out_az    = '/srv/data/dva/survey_azimuth_scans_noise_corr/'\n",
    "###########################################################"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae57e12a",
   "metadata": {},
   "source": [
    "## Read in files listing start and stop times for azimuth and raster scans:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c9132ff4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "raster 1: ['2022-06-30T02:40:02.316993'] ['2022-06-30T03:37:07.839786']\n",
      "\n",
      "2491 2022-06-30T03:54:28.153 2022-06-30T04:12:30.329\n",
      "2492 2022-06-30T04:14:12.932 2022-06-30T04:32:15.053\n",
      "1197 2022-06-30T04:37:33.211 2022-06-30T04:55:35.319\n",
      "1198 2022-06-30T04:57:18.001 2022-06-30T05:15:20.323\n",
      "1199 2022-06-30T05:17:38.664 2022-06-30T05:35:41.006\n",
      "2712 2022-06-30T05:40:22.475 2022-06-30T05:58:24.550\n",
      "1417 2022-06-30T06:03:06.860 2022-06-30T06:21:09.224\n",
      "1418 2022-06-30T06:23:27.746 2022-06-30T06:41:29.518\n",
      "1493 2022-06-30T06:46:11.995 2022-06-30T07:04:13.791\n",
      "1494 2022-06-30T07:05:56.687 2022-06-30T07:23:58.984\n",
      "1495 2022-06-30T07:25:41.527 2022-06-30T07:43:43.927\n",
      "1496 2022-06-30T07:46:02.123 2022-06-30T08:04:03.981\n",
      "1497 2022-06-30T08:05:47.068 2022-06-30T08:23:48.912\n",
      "1498 2022-06-30T08:25:31.642 2022-06-30T08:43:33.545\n",
      "0203 2022-06-30T08:48:51.641 2022-06-30T09:06:53.689\n",
      "1716 2022-06-30T09:11:36.149 2022-06-30T09:29:38.011\n",
      "1717 2022-06-30T09:31:20.914 2022-06-30T09:49:22.944\n",
      "1718 2022-06-30T09:51:41.198 2022-06-30T10:09:43.126\n",
      "1719 2022-06-30T10:11:26.268 2022-06-30T10:29:28.179\n",
      "1864 2022-06-30T10:35:22.433 2022-06-30T10:53:24.462\n",
      "1865 2022-06-30T10:55:07.155 2022-06-30T11:13:09.413\n",
      "1866 2022-06-30T11:15:27.592 2022-06-30T11:33:29.927\n",
      "0571 2022-06-30T11:38:12.079 2022-06-30T11:56:14.399\n",
      "\n",
      "raster 2: ['2022-06-30T11:59:40.290185'] ['2022-06-30T12:56:47.138683']\n"
     ]
    }
   ],
   "source": [
    "scan_id = []    # The scan id number\n",
    "scan_start = []  # Start time of the scan (UTC)\n",
    "scan_stop = []   # Stop time of the scan (UTC)\n",
    "\n",
    "raster1_start = []\n",
    "raster2_start = []\n",
    "raster1_stop = []\n",
    "raster2_stop = []\n",
    "\n",
    "rast_list = ['a','b','c']\n",
    "\n",
    "# Read in the azimuth scan data and store it in arrays:\n",
    "with open(dir_in_az+'DVAsurvey_phase'+phase+'_day0'+day+'.txt') as fp:\n",
    "    for line in fp:       \n",
    "        scan_id.append(int(line.split()[0]))\n",
    "        scan_start.append(line.split()[1]+'T'+line.split()[2][0:12])\n",
    "        scan_stop.append(line.split()[3]+'T'+line.split()[4][0:12])\n",
    "              \n",
    "for i in range(0,Ne):\n",
    "    \n",
    "    if Ne == 1:\n",
    "        rast1name = 'dva_survey_phase'+phase+'_day0'+day+'_raster1.txt'\n",
    "    else:\n",
    "        rast1name = 'dva_survey_phase'+phase+'_day0'+day+'_raster1'+rast_list[i]+'.txt'\n",
    "        \n",
    "    with open(dir_in_rast+rast1name) as fp:\n",
    "        for line in fp:  \n",
    "            raster1_start.append(line.split()[3])\n",
    "            raster1_stop.append(line.split()[4])\n",
    "            \n",
    "for i in range(0,Nm):\n",
    "    \n",
    "    if Nm == 1:\n",
    "        rast2name = 'dva_survey_phase'+phase+'_day0'+day+'_raster2.txt'\n",
    "    else:\n",
    "        rast2name = 'dva_survey_phase'+phase+'_day0'+day+'_raster2'+rast_list[i]+'.txt'\n",
    "        \n",
    "    with open(dir_in_rast+rast2name) as fp:\n",
    "        for line in fp:  \n",
    "            raster2_start.append(line.split()[3])\n",
    "            raster2_stop.append(line.split()[4])\n",
    "                \n",
    "print('raster 1:',raster1_start,raster1_stop)\n",
    "print('')\n",
    "for i in range(0,len(scan_id)):\n",
    "    print(f\"{scan_id[i]:04}\",scan_start[i],scan_stop[i])\n",
    "print('')\n",
    "print('raster 2:',raster2_start,raster2_stop)\n",
    "\n",
    "# Convert start and stop times to Modified Julian Day (MJD).\n",
    "scan_start_mjd = Time(scan_start, format='isot',scale='utc').mjd\n",
    "scan_stop_mjd  = Time(scan_stop,  format='isot',scale='utc').mjd\n",
    "raster1_start_mjd = Time(raster1_start, format='isot',scale='utc').mjd\n",
    "raster1_stop_mjd  = Time(raster1_stop,  format='isot',scale='utc').mjd\n",
    "raster2_start_mjd = Time(raster2_start, format='isot',scale='utc').mjd\n",
    "raster2_stop_mjd  = Time(raster2_stop,  format='isot',scale='utc').mjd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78e4c5eb",
   "metadata": {},
   "source": [
    "## Read in data and stitch together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "fb10fd8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def concatenate_data(file,RR,LL,reRL,imRL,dec,ra,el,az,t,noise):\n",
    "    \n",
    "    dataset = file['data']['beam_0']['band_SB0']['scan_0']\n",
    "    \n",
    "    idx1 = len(t)\n",
    "    \n",
    "    # Add the position and time data to the corresponding arrays:\n",
    "    dec = np.concatenate([dec,dataset['metadata']['declination']])\n",
    "    ra = np.concatenate([ra,dataset['metadata']['right_ascension']])\n",
    "    el = np.concatenate([el,dataset['metadata']['elevation']])\n",
    "    az = np.concatenate([az,dataset['metadata']['azimuth']])\n",
    "    t = np.concatenate([t,dataset['metadata']['utc']])\n",
    "    noise = np.concatenate([noise,dataset['metadata']['noise_state']])\n",
    "\n",
    "    idx2 = len(t)\n",
    "    \n",
    "    # Add the spectrometer power data to the corresponding arrays:\n",
    "    RR = np.concatenate([RR,dataset['data'][:,0,:]],axis=0)\n",
    "    LL = np.concatenate([LL,dataset['data'][:,1,:]],axis=0)\n",
    "    reRL = np.concatenate([reRL,dataset['data'][:,2,:]],axis=0)\n",
    "    imRL = np.concatenate([imRL,dataset['data'][:,3,:]],axis=0)\n",
    "    \n",
    "    return RR,LL,reRL,imRL,dec,ra,el,az,t,noise,idx1,idx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a4f4cab0",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-26-166fc0cc7ced>\u001b[0m in \u001b[0;36mconcatenate_data\u001b[0;34m(file, RR, LL, reRL, imRL, dec, ra, el, az, t, noise)\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mLL\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mLL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m     \u001b[0mreRL\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mreRL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m     \u001b[0mimRL\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mimRL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mRR\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mLL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mreRL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mimRL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdec\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mra\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maz\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnoise\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0midx1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0midx2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mh5py/_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mh5py/_objects.pyx\u001b[0m in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m/UBC-O/ashill01/anaconda3/envs/env_2021-05/lib/python3.9/site-packages/h5py/_hl/dataset.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, args, new_dtype)\u001b[0m\n\u001b[1;32m    788\u001b[0m         \u001b[0mmspace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh5s\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcreate_simple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mselection\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    789\u001b[0m         \u001b[0mfspace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mselection\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mid\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 790\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mid\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmspace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfspace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdxpl\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_dxpl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    791\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    792\u001b[0m         \u001b[0;31m# Patch up the output for NumPy\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "t = []\n",
    "az = []\n",
    "dec = []\n",
    "ra = []\n",
    "el = []\n",
    "noise = []\n",
    "idx1_rast1 = []\n",
    "idx2_rast1 = []\n",
    "idx1_rast2 = []\n",
    "idx2_rast2 = []\n",
    "\n",
    "# Use one of the scans to get the list of frequencies:\n",
    "scan0 = f\"{scan_id[0]:04}\"\n",
    "file = h5py.File(dir_in_az+'dva_survey_phase'+phase+'_raw_'+scan0+'.h5','r')\n",
    "freq = file['data']['beam_0']['band_SB0']['frequency'][:]/1e6\n",
    "\n",
    "# Create empty arrays for the power data:\n",
    "RR = np.empty([0,len(freq)])\n",
    "LL = np.empty([0,len(freq)])\n",
    "reRL = np.empty([0,len(freq)])\n",
    "imRL = np.empty([0,len(freq)])\n",
    "\n",
    "# Raster scan 1:\n",
    "for i in range(0,Ne):\n",
    "    if Ne == 1:\n",
    "        file = h5py.File(dir_in_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster1'+'.h5','r')\n",
    "    else:\n",
    "        file = h5py.File(dir_in_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster1'+rast_list[i]+'.h5','r')\n",
    "\n",
    "    rast1_file_t = [Time(file['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][0],\n",
    "                        format='isot',scale='utc').mjd,\n",
    "                    Time(file['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][-1],\n",
    "                        format='isot',scale='utc').mjd]\n",
    "    RR,LL,reRL,imRL,dec,ra,el,az,t,noise,idx1,idx2 = concatenate_data(file,RR,LL,reRL,imRL,\n",
    "                                                            dec,ra,el,az,t,noise)\n",
    "    print(idx1,idx2)\n",
    "    idx1_rast1.append(idx1)\n",
    "    idx2_rast1.append(idx2)\n",
    "\n",
    "# Loop through all the scans in the \"scan_num\" list:\n",
    "for i in scan_id:\n",
    "    file = h5py.File(dir_in_az+'dva_survey_phase'+phase+'_raw_'+f\"{i:04}\"+'.h5','r')\n",
    "    print(i,file)\n",
    "    RR,LL,reRL,imRL,dec,ra,el,az,t,noise,idx1,idx2 = concatenate_data(file,RR,LL,reRL,imRL,\n",
    "                                                            dec,ra,el,az,t,noise)\n",
    "    print(idx1,idx2)\n",
    "    \n",
    "# Raster scan 2:\n",
    "for i in range(0,Nm):\n",
    "    if Nm == 1:\n",
    "        file = h5py.File(dir_in_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster2'+'.h5','r')\n",
    "    else:\n",
    "        file = h5py.File(dir_in_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster2'+rast_list[i]+'.h5','r')\n",
    "\n",
    "    rast2_file_t = [Time(file['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][0],\n",
    "                        format='isot',scale='utc').mjd,\n",
    "                    Time(file['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][-1],\n",
    "                        format='isot',scale='utc').mjd]\n",
    "    RR,LL,reRL,imRL,dec,ra,el,az,t,noise,idx1,idx2 = concatenate_data(file,RR,LL,reRL,imRL,\n",
    "                                                            dec,ra,el,az,t,noise)\n",
    "    print(idx1,idx2)\n",
    "    idx1_rast2.append(idx1)\n",
    "    idx2_rast2.append(idx2)\n",
    "\n",
    "t_set_plt = Time(t, format='isot',scale='utc').mjd\n",
    "\n",
    "print(RR.shape)\n",
    "print(t_set_plt.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "e4a70378",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[59760.11062500001, 59760.15228472224]\n",
      "[]\n",
      "[]\n",
      "\n",
      "[59759.5022569445, 59759.5466944445]\n",
      "[]\n",
      "[]\n",
      "\n",
      "[59759.10434028 59759.10434722 59759.10435417 ... 59759.54668056\n",
      " 59759.5466875  59759.54669444]\n",
      "59759.32992716642\n"
     ]
    }
   ],
   "source": [
    "print(rast1_file_t)\n",
    "print(idx1_rast1)\n",
    "print(idx2_rast1)\n",
    "print('')\n",
    "print(rast2_file_t)\n",
    "print(idx1_rast2)\n",
    "print(idx2_rast2)\n",
    "print('')\n",
    "print(t_set_plt)\n",
    "t_mean = np.nanmean(t_set_plt)\n",
    "print(t_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c97cdc3",
   "metadata": {},
   "source": [
    "## Read in temperature data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "6d15cacf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def month_to_num(month_name):\n",
    "    if month_name == 'Jan': month_num = '01'\n",
    "    if month_name == 'Feb': month_num = '02'\n",
    "    if month_name == 'Mar': month_num = '03'\n",
    "    if month_name == 'Apr': month_num = '04'\n",
    "    if month_name == 'May': month_num = '05'\n",
    "    if month_name == 'Jun': month_num = '06'\n",
    "    if month_name == 'Jul': month_num = '07'\n",
    "    if month_name == 'Aug': month_num = '08'\n",
    "    if month_name == 'Sep': month_num = '09'\n",
    "    if month_name == 'Oct': month_num = '10'\n",
    "    if month_name == 'Nov': month_num = '11'\n",
    "    if month_name == 'Dec': month_num = '12'\n",
    "    return(month_num)\n",
    "\n",
    "i = 0\n",
    "t_weath = []\n",
    "temp_C = []\n",
    "\n",
    "with open('/srv/data/dva/weather_data/weather_survey_phase'+phase+'.txt') as fp:\n",
    "    for line in fp:\n",
    "        t_weath.append(str( line.split()[2]+'-'+month_to_num(line.split()[1])+'-'+line.split()[0]+\n",
    "                  'T'+line.split()[3]))\n",
    "        temp_C.append(line.split()[4])        \n",
    "\n",
    "temp_C = np.array(temp_C,dtype=float)\n",
    "t_weath_fix = Time(t_weath, format='isot',scale='utc')\n",
    "t_weath_plt = t_weath_fix.mjd\n",
    "\n",
    "#plt.plot(t_weath_plt,temp_C)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39248381",
   "metadata": {},
   "source": [
    "## New code to calculate the power deflection at each noise source instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "42a618a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "################################################\n",
    "# The number of integrations to include on either \n",
    "# side of each noise source instance:\n",
    "n_off = 5\n",
    "################################################\n",
    "\n",
    "# Make arrays for noise source deflection:\n",
    "LL_dnoise = []\n",
    "RR_dnoise = []\n",
    "reRL_dnoise = []\n",
    "imRL_dnoise = []\n",
    "t_noise = []\n",
    "\n",
    "wnoise = np.where(noise == 1)[0]\n",
    "print(wnoise)\n",
    "for k,g in groupby(enumerate(wnoise),lambda x:x[0]-x[1]):\n",
    "\n",
    "    try:\n",
    "\n",
    "        group = np.array(list(map(itemgetter(1),g)))\n",
    "        #print(group)\n",
    "        #print(group[1:-1])\n",
    "        #middle = [group[int(np.floor((len(group)-1)/2))] ,group[int(np.ceil((len(group)-1)/2))]]\n",
    "        #print(middle)\n",
    "        offleft = [group[0]-n_off,group[0]-2]\n",
    "        offright = [group[-1]+2,group[-1]+n_off]\n",
    "        #print(offleft)\n",
    "        #print(offright)\n",
    "        #print('')\n",
    "    \n",
    "        LL_noise = np.nanmedian(LL[group[1:-1],:],axis=0)\n",
    "        LL_off = (np.nanmedian(LL[offleft[0]:offleft[-1],:],axis=0) + np.nanmedian(LL[offright[0]:offright[-1],:],axis=0))/2.\n",
    "    \n",
    "        RR_noise = np.nanmedian(RR[group[1:-1],:],axis=0)\n",
    "        RR_off = (np.nanmedian(RR[offleft[0]:offleft[-1],:],axis=0) + np.nanmedian(RR[offright[0]:offright[-1],:],axis=0))/2.\n",
    "    \n",
    "        reRL_noise = np.nanmedian(reRL[group[1:-1],:],axis=0)\n",
    "        reRL_off = (np.nanmedian(reRL[offleft[0]:offleft[-1],:],axis=0) + np.nanmedian(reRL[offright[0]:offright[-1],:],axis=0))/2.\n",
    "    \n",
    "        imRL_noise = np.nanmedian(imRL[group[1:-1],:],axis=0)\n",
    "        imRL_off = (np.nanmedian(imRL[offleft[0]:offleft[-1],:],axis=0) + np.nanmedian(imRL[offright[0]:offright[-1],:],axis=0))/2.\n",
    "    \n",
    "        t_noise.append(np.nanmedian(t_set_plt[group[1:-1]]))\n",
    "\n",
    "        LL_dnoise.append(LL_noise - LL_off)\n",
    "        RR_dnoise.append(RR_noise - RR_off)\n",
    "        reRL_dnoise.append(reRL_noise - reRL_off)\n",
    "        imRL_dnoise.append(imRL_noise - imRL_off)\n",
    "        \n",
    "    except:\n",
    "        pass\n",
    "\n",
    "LL_dnoise_arr = np.empty([len(LL_dnoise),len(LL_dnoise[0])])\n",
    "RR_dnoise_arr = np.empty([len(RR_dnoise),len(RR_dnoise[0])])\n",
    "reRL_dnoise_arr = np.empty([len(reRL_dnoise),len(reRL_dnoise[0])])\n",
    "imRL_dnoise_arr = np.empty([len(imRL_dnoise),len(imRL_dnoise[0])])\n",
    "\n",
    "for i in range(0,len(LL_dnoise)):\n",
    "    LL_dnoise_arr[i,:] = LL_dnoise[i]\n",
    "    RR_dnoise_arr[i,:] = RR_dnoise[i]\n",
    "    reRL_dnoise_arr[i,:] = reRL_dnoise[i]\n",
    "    imRL_dnoise_arr[i,:] = imRL_dnoise[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d840da8",
   "metadata": {},
   "source": [
    "## Check that the noise source was being fired where expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "452d0629",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x and y must be the same size",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-31-c7be6a7ca384>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mii\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m7000\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt_set_plt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnoise\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_xlim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mraster1_start_mjd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mraster1_stop_mjd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0maxs2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtwinx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/UBC-O/ashill01/anaconda3/envs/env_2021-05/lib/python3.9/site-packages/matplotlib/__init__.py\u001b[0m in \u001b[0;36minner\u001b[0;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1410\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1411\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1412\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msanitize_sequence\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1413\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1414\u001b[0m         \u001b[0mbound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnew_sig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbind\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/UBC-O/ashill01/anaconda3/envs/env_2021-05/lib/python3.9/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mscatter\u001b[0;34m(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, edgecolors, plotnonfinite, **kwargs)\u001b[0m\n\u001b[1;32m   4367\u001b[0m         \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4368\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4369\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"x and y must be the same size\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4370\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4371\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0ms\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: x and y must be the same size"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6oAAAHWCAYAAAB+A3SNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwh0lEQVR4nO3df6hn910n/uerGZPCWkdwuqxkZkzAKTibCuleYqV/tJDuMgnMDHx1S/KlaCV0+IIRdxXZiFIlsrBV1EWIP+5iqRZsjP1DZnUkQo2UFRNyQ3dDk5JyGZdmopCkxvtPaWPW1/ePz6fL9eYm99yZz+dzz7n38YCBzznn/fmcJ4fXXO7rnvN5v6u7AwAAAGPxjoMOAAAAANtpVAEAABgVjSoAAACjolEFAABgVDSqAAAAjIpGFQAAgFHZs1Gtqk9V1ctV9aW3OF5V9RtVtVlVz1bV+xYfEwAAgKNiyB3VTyc59zbH70lyZv7vUpLfGnLiqro0ZByMnVrmMFDHHBZqmcNAHXNY3Egt79modvcXkvz92wy5mOT3e+bJJN9ZVd894NyT+w9YVecPOsN+TC1vMs3MmVgtT/Eay7wSk6rjZHrXeGp5k2lmzsRqeYrXWOaVmFQdJ5O8xpPLPLW8c8trVAe4NcmL27avzfcdRlMrjqnlTaaZeWqmeI1lZjdTu8ZTy5tMM/PUTPEay8xupniNp5Z5anlvSHX33oOqbkvyJ919xy7H/iTJf+nu/zHf/nyS/9TdG7uMvZR5V33LLbf8mzvueNPHjdrW1laOHz9+0DEGm1reZJqZX3nllbz73e8+6BiDTfEay7x8U6vjZHrXeGp5k2lmnlotT/Eay7x8U6vjZHrXOJle5qnlTZJnnnnmm0m2z3W03t3rQ957bAHnfynJqW3bJ+f73mQeaj1J1tbWemPjTb0sAAAAh0BVfam7167nvYt49Pdykh+Zz/77/iRb3f13C/hcAAAAjqA976hW1WeTfCjJiaq6luQXknxbknT3bye5kuTeJJtJvp7kx5YVFgAAgMNvz0a1u+/f43gn+fGFJQIAAOBIW8SjvwAAALAwK29Uq+p8Va1vbW2t+tQAAACszvGqWr+eNWAHLU+zDGb9BQAAOLyq6pmDnPUXAAAAFkajCgAAwKhoVAEAABgVjSoAAACjolEFAABgVCxPAwAAwDJYngYAAIDxsDwNAAAAh4ZGFQAAgFHRqAIAADAqGlUAAABGxay/AAAALINZfwEAABgPs/4CAABwaAxqVKvqXFW9UFWbVfXQLsdPV9UTVfXFqnq2qu5dfFQAAACOgj0b1aq6KckjSe5JcjbJ/VV1dsewn0/yWHffmeS+JL+56KAAAAAcDUPuqN6VZLO7r3b360keTXJxx5hO8h3z18eT/O3iIgIAAHCUHBsw5tYkL27bvpbkB3aM+cUkf15VP5HkXyT58ELSAQAAcOQsajKl+5N8urtPJrk3yWeq6k2fXVWXqmqjqjZeeeWVBZ0aAACAETrxrf5v/u/S0DcOuaP6UpJT27ZPzvdt90CSc0nS3X9dVe9MciLJy9sHdfd6kvVktjzN0JAAAABMzqvLXJ7m6SRnqur2qro5s8mSLu8Y89UkdydJVX1fkncmccsUAACAfduzUe3uN5I8mOTxJF/ObHbf56rq4aq6MB/200k+XlX/K8lnk3ysu90xBQAAYN+GPPqb7r6S5MqOfZ/Y9vr5JB9YbDQAAACOokVNpgQAAAALsfJGtarOV9X61tbWqk8NAADA6hyvqvWqOr/fN9ZBfZV0bW2tNzY2DuTcAAAALFdVPbPMWX8BAABgZTSqAAAAjIpGFQAAgFHRqAIAADAqZv0FAABgGcz6CwAAwHiY9RcAAIBDQ6MKAADAqGhUAQAAGBWNKgAAAKOiUQUAAGBULE8DAADAMlieBgAAgPFY+vI0VXWuql6oqs2qeugtxnykqp6vqueq6g+uJwwAAAAc22tAVd2U5JEk/zbJtSRPV9Xl7n5+25gzSX42yQe6+7Wq+pfLCgwAAMDhNuSO6l1JNrv7ane/nuTRJBd3jPl4kke6+7Uk6e6XFxsTAACAo2JIo3prkhe3bV+b79vuPUneU1V/VVVPVtW5RQUEAADgaNnz0d99fM6ZJB9KcjLJF6rqvd39D9sHVdWlJJeS5PTp0ws6NQAAACN0oqq2z6C73t3rQ944pFF9Kcmpbdsn5/u2u5bkqe7+xyR/U1VfyaxxfXr7oHmo9WQ26++QgAAAAEzSq8uc9ffpJGeq6vaqujnJfUku7xjzx5ndTU1VncjsUeCr1xMIAACAo23PRrW730jyYJLHk3w5yWPd/VxVPVxVF+bDHk/ytap6PskTSX6mu7+2rNAAAAAcXtV9ME/grq2t9cbGxt4DAQAAmJyqemaZj/4CAADAyqy8Ua2q81W1vrW1tepTAwAAsDrHq2q9qs7v940e/QUAAGDhPPoLAADAoaFRBQAAYFQ0qgAAAIyKRhUAAIBR0agCAAAwKpanAQAAYBksTwMAAMB4WJ4GAACAQ0OjCgAAwKhoVAEAABgVjSoAAACjYtZfAAAAlsGsvwAAAIyHWX8BAAA4NAY1qlV1rqpeqKrNqnrobcb9UFV1VV1X1wwAAAB7NqpVdVOSR5Lck+Rskvur6uwu496V5CeTPLXokAAAABwdQ+6o3pVks7uvdvfrSR5NcnGXcb+U5JNJvrHAfAAAABwxQxrVW5O8uG372nzf/1VV70tyqrv/dIHZAAAAOIKO3egHVNU7kvxako8NGHspyaUkOX369I2eGgAAgPE6UVXbl3pZ7+71IW8c0qi+lOTUtu2T833f8q4kdyT5y6pKkn+V5HJVXejuf7b+zDzUejJbnmZIQAAAACbp1WUuT/N0kjNVdXtV3ZzkviSXv3Wwu7e6+0R339bdtyV5MsmbmlQAAAAYYs9GtbvfSPJgkseTfDnJY939XFU9XFUXlh0QAACAo2XQd1S7+0qSKzv2feItxn7oxmMBAABwVA159BcAAABWZuWNalWdr6r1ra2tVZ8aAACA1TleVetVdX6/b6zug5l8d21trTc2zLcEAABwGFXVM8uc9RcAAABWRqMKAADAqGhUAQAAGBWNKgAAAKNi1l8AAACWway/AAAAjIdZfwEAADg0NKoAAACMikYVAACAUdGoAgAAMCoaVQAAAEbF8jQAAAAsg+VpAAAAGA/L0wAAAHBoDGpUq+pcVb1QVZtV9dAux3+qqp6vqmer6vNV9T2LjwoAAMBRsGejWlU3JXkkyT1Jzia5v6rO7hj2xSRr3f39ST6X5JcXHRQAAICjYcgd1buSbHb31e5+PcmjSS5uH9DdT3T31+ebTyY5udiYAAAAHBVDGtVbk7y4bfvafN9beSDJn+12oKouVdVGVW288sorw1MCAAAwNSe+1f/N/10a+sZji0xRVR9Nspbkg7sd7+71JOvJbNbfRZ4bAACAUXn1emf9HdKovpTk1Lbtk/N9/0xVfTjJzyX5YHd/83rCAAAAwJBHf59Ocqaqbq+qm5Pcl+Ty9gFVdWeS30lyobtfXnxMAAAAjoo9G9XufiPJg0keT/LlJI9193NV9XBVXZgP+5Uk357kj6rqf1bV5bf4OAAAAHhbg76j2t1XklzZse8T215/eMG5AAAAOKKGPPoLAAAAK7PyRrWqzlfV+tbW1qpPDQAAwOocr6r1qjq/3zdW98GsErO2ttYbGxsHcm4AAACWq6qeud7laTz6CwAAwKhoVAEAABgVjSoAAACjolEFAABgVMz6CwAAwDKY9RcAAIDxMOsvAAAAh4ZGFQAAgFHRqAIAADAqGlUAAABGxay/AAAALINZfwEAABgPs/4CAABwaAxqVKvqXFW9UFWbVfXQLsdvqao/nB9/qqpuW3hSAAAAjoQ9G9WquinJI0nuSXI2yf1VdXbHsAeSvNbd35vk15N8ctFBAQAAOBqG3FG9K8lmd1/t7teTPJrk4o4xF5P83vz155LcXVW1uJgAAAAcFUMa1VuTvLht+9p8365juvuNJFtJvmsRAQEAADhajq3yZFV1Kcml+eY3q+pLqzz/AhzPrAmfiqnlTaaZ+USSVw86xD5M8RrLvHxTq+Nketd4anmTaWaeWi1P8RrLvHxTq+Nketc4mV7mqeVNkjuqavtSL+vdvT7kjUMa1ZeSnNq2fXK+b7cx16rqWGYX8Ws7P2geaj1JqmrjeqcqPihVtd7dl/YeOQ5Ty5tMNvOkanmi11jmJZtaHSeTvMaTyptMNvOkanmi11jmJZtaHSfTu8bJ9DJPLW9yY7U85NHfp5Ocqarbq+rmJPclubxjzOUkPzp//cNJ/qIPaoHW5frvBx1gn6aWN5lm5qmZ4jWWmd1M7RpPLW8yzcxTM8VrLDO7meI1nlrmqeW9ITWkn6yqe5P81yQ3JflUd//nqno4yUZ3X66qdyb5TJI7k/x9kvu6++oenzm5vxTBbtQyh4E65rBQyxwG6pjD4kZqedB3VLv7SpIrO/Z9YtvrbyT59/s896Bnk2EC1DKHgTrmsFDLHAbqmMPiumt50B1VAAAAWJUh31EFAACAlVl6o1pV56rqhararKqHdjl+S1X94fz4U1V127IzwX4NqOOfqqrnq+rZqvp8VX3PQeSEvexVy9vG/VBVdVX5jhSjM6SOq+oj85/Lz1XVH6w6Iwwx4PeL01X1RFV9cf47xr0HkRPeTlV9qqpefqulR2vmN+Z1/mxVvW/I5y61Ua2qm5I8kuSeJGeT3F9VZ3cMeyDJa939vUl+Pcknl5kJ9mtgHX8xyVp3f3+SzyX55dWmhL0NrOVU1buS/GSSp1abEPY2pI6r6kySn03yge7+10n+w6pzwl4G/kz++SSPdfedma288ZurTQmDfDrJubc5fk+SM/N/l5L81pAPXfYd1buSbHb31e5+PcmjSS7uGHMxye/NX38uyd1VVUvOBfuxZx139xPd/fX55pOZrTcMYzPkZ3KS/FJmfzT8xirDwUBD6vjjSR7p7teSpLtfXnFGGGJILXeS75i/Pp7kb1eYDwbp7i9ktvLLW7mY5Pd75skk31lV373X5y67Ub01yYvbtq/N9+06prvfSLKV5LuWnAv2Y0gdb/dAkj9baiK4PnvW8vxxnFPd/aerDAb7MORn8nuSvKeq/qqqnqyqt/tLPxyUIbX8i0k+WlXXMluB4ydWEw0War+/SycZuDwNMExVfTTJWpIPHnQW2K+qekeSX0vysQOOAjfqWGaPmH0osydcvlBV7+3ufzjIUHAd7k/y6e7+1ar6wSSfqao7uvufDjoYLNuy76i+lOTUtu2T8327jqmqY5k91vC1JeeC/RhSx6mqDyf5uSQXuvubK8oG+7FXLb8ryR1J/rKq/neS9ye5bEIlRmbIz+RrSS539z92998k+UpmjSuMyZBafiDJY0nS3X+d5J1JTqwkHSzOoN+ld1p2o/p0kjNVdXtV3ZzZl8Av7xhzOcmPzl//cJK/aIu7Mi571nFV3ZnkdzJrUn0XirF621ru7q3uPtHdt3X3bZl93/pCd28cTFzY1ZDfLf44s7upqaoTmT0KfHWFGWGIIbX81SR3J0lVfV9mjeorK00JN+5ykh+Zz/77/iRb3f13e71pqY/+dvcbVfVgkseT3JTkU939XFU9nGSjuy8n+d3MHmPYzOxLuPctMxPs18A6/pUk357kj+ZzgX21uy8cWGjYxcBahlEbWMePJ/l3VfV8kv+T5Ge629NajMrAWv7pJP+tqv5jZhMrfcwNHcamqj6b2R8HT8y/T/0LSb4tSbr7tzP7fvW9STaTfD3Jjw36XLUOAADAmCz70V8AAADYF40qAAAAo6JRBQAAYFQ0qgAAAIyKRhUAAIBR0agCAAAwKns2qlX1qap6uaq+9BbHq6p+o6o2q+rZqnrfkBNX1aX9hoUxUsscBuqYw0ItcxioYw6LG6nlIXdUP53k3NscvyfJmfm/S0l+a+C5J/cfsKrOH3SG/Zha3mSamTOxWp7iNZZ5JSZVx8n0rvHU8ibTzJyJ1fIUr7HMKzGpOk4meY0nl3lqeeeW16h29xeS/P3bDLmY5Pd75skk31lV3329gUZuasUxtbzJNDNPzRSvsczsZmrXeGp5k2lmnpopXmOZ2c0Ur/HUMk8t7w2p7t57UNVtSf6ku+/Y5difJPkv3f0/5tufT/Kfuntjl7GXMu+qb7nlln9zxx1v+rhR29rayvHjxw86xmBTy5tMM/Mrr7ySd7/73QcdY7ApXmOZl29qdZxM7xpPLW8yzcxTq+UpXmOZl29qdZxM7xon08s8tbxJ8swzz3wzyfavkK539/qQ9x5bTqTdzUOtJ8na2lpvbLyplwUAAOAQqKovdffa9bx3EbP+vpTk1Lbtk/N9AAAAsG+LaFQvJ/mR+ey/70+y1d1/t4DPBQAA4Aja89Hfqvpskg8lOVFV15L8QpJvS5Lu/u0kV5Lcm2QzydeT/NiywgIAAHD47dmodvf9exzvJD++sEQAAAAcaYt49BcAAAAWZuWNalWdr6r1ra2tVZ8aAACA1TleVetVte81YAeto7oMlqcBAAA4vKrqmYNcngYAAAAWRqMKAADAqGhUAQAAGBWNKgAAAKNi1l8AAACWway/AAAAjIdZfwEAADg0NKoAAACMikYVAACAUdGoAgAAMCoaVQAAAEbF8jQAAAAsg+VpAAAAGI+lL09TVeeq6oWq2qyqh3Y5frqqnqiqL1bVs1V17/WEAQAAgD0b1aq6KckjSe5JcjbJ/VV1dsewn0/yWHffmeS+JL+56KAAAAAcDUPuqN6VZLO7r3b360keTXJxx5hO8h3z18eT/O3iIgIAAHCUHBsw5tYkL27bvpbkB3aM+cUkf15VP5HkXyT58ELSAQAAcOQsatbf+5N8urtPJrk3yWeq6k2fXVWXqmqjqjZeeeWVBZ0aAACAETrxrf5v/u/S0DcOuaP6UpJT27ZPzvdt90CSc0nS3X9dVe9MciLJy9sHdfd6kvVkNuvv0JAAAABMzqvLnPX36SRnqur2qro5s8mSLu8Y89UkdydJVX1fkncmccsUAACAfduzUe3uN5I8mOTxJF/ObHbf56rq4aq6MB/200k+XlX/K8lnk3ysD2qBVgAAACZtyKO/6e4rSa7s2PeJba+fT/KBxUYDAADgKFrUZEoAAACwECtvVKvqfFWtb21trfrUAAAArM7xqlqvqvP7fWMd1FdJ19bWemNj40DODQAAwHJV1TPLnPUXAAAAVkajCgAAwKhoVAEAABgVjSoAAACjolEFAABgVCxPAwAAwDJYngYAAIDxsDwNAAAAh4ZGFQAAgFHRqAIAADAqGlUAAABGxay/AAAALINZfwEAABgPs/4CAABwaAxqVKvqXFW9UFWbVfXQW4z5SFU9X1XPVdUfLDYmAAAAR8WxvQZU1U1JHknyb5NcS/J0VV3u7ue3jTmT5GeTfKC7X6uqf7mswAAAABxuQ+6o3pVks7uvdvfrSR5NcnHHmI8neaS7X0uS7n55sTEBAAA4KoY0qrcmeXHb9rX5vu3ek+Q9VfVXVfVkVZ1bVEAAAACOlj0f/d3H55xJ8qEkJ5N8oare293/sH1QVV1KcilJTp8+vaBTAwAAMEInqmr7Ui/r3b0+5I1DGtWXkpzatn1yvm+7a0me6u5/TPI3VfWVzBrXp7cPmodaT2bL0wwJCAAAwCS9uszlaZ5Ocqaqbq+qm5Pcl+TyjjF/nNnd1FTVicweBb56PYEAAAA42vZsVLv7jSQPJnk8yZeTPNbdz1XVw1V1YT7s8SRfq6rnkzyR5Ge6+2vLCg0AAMDhVd0H8wTu2tpab2xs7D0QAACAyamqZ5b56C8AAACszMob1ao6X1XrW1tbqz41AAAAq3O8qtar6vx+3+jRXwAAABbOo78AAAAcGhpVAAAARkWjCgAAwKhoVAEAABgVs/4CAACwDGb9BQAAYDzM+gsAAMChoVEFAABgVDSqAAAAjIpGFQAAgFHRqAIAADAqlqcBAABgGSxPAwAAwHgsfXmaqjpXVS9U1WZVPfQ2436oqrqqrisMAAAA7NmoVtVNSR5Jck+Ss0nur6qzu4x7V5KfTPLUokMCAABwdAy5o3pXks3uvtrdryd5NMnFXcb9UpJPJvnGAvMBAABwxAxpVG9N8uK27Wvzff9XVb0vyanu/tMFZgMAAOAIOnajH1BV70jya0k+NmDspSSXkuT06dM3emoAAADG60RVbZ9Bd72714e8cUij+lKSU9u2T873fcu7ktyR5C+rKkn+VZLLVXWhu//ZtL7zUOvJbNbfIQEBAACYpFeXOevv00nOVNXtVXVzkvuSXP7Wwe7e6u4T3X1bd9+W5Mkkb2pSAQAAYIg9G9XufiPJg0keT/LlJI9193NV9XBVXVh2QAAAAI6WQd9R7e4rSa7s2PeJtxj7oRuPBQAAwFE15NFfAAAAWBmNKgAAAKOy8ka1qs5X1frW1taqTw0AAMDqHK+q9ao6v983VvfBrBKztrbWGxsmBgYAADiMquqZZS5PAwAAACujUQUAAGBUNKoAAACMikYVAACAUTHrLwAAAMtg1l8AAADGw6y/AAAAHBoaVQAAAEZFowoAAMCoaFQBAAAYFbP+AgAAsAxm/QUAAGA8zPoLAADAoTGoUa2qc1X1QlVtVtVDuxz/qap6vqqerarPV9X3LD4qAAAAR8GejWpV3ZTkkST3JDmb5P6qOrtj2BeTrHX39yf5XJJfXnRQAAAAjoYhd1TvSrLZ3Ve7+/Ukjya5uH1Adz/R3V+fbz6Z5ORiYwIAAHBUDGlUb03y4rbta/N9b+WBJH92I6EAAAA4uo4t8sOq6qNJ1pJ88C2OX0pyKUlOnz69yFMDAAAwLieqavtSL+vdvT7kjUMa1ZeSnNq2fXK+75+pqg8n+bkkH+zub+72QfNQ68lseZohAQEAAJikV5e5PM3TSc5U1e1VdXOS+5Jc3j6gqu5M8jtJLnT3y9cTBAAAAJIBjWp3v5HkwSSPJ/lykse6+7mqeriqLsyH/UqSb0/yR1X1P6vq8lt8HAAAALytQd9R7e4rSa7s2PeJba8/vOBcAAAAHFFDHv0FAACAlVl5o1pV56tqfWtra9WnBgAAYHWOV9V6VZ3f7xur+2Am311bW+uNjY29BwIAADA5VfXMMmf9BQAAgJXRqAIAADAqGlUAAABGRaMKAADAqGhUAQAAGBXL0wAAALAMlqcBAABgPCxPAwAAwKGhUQUAAGBUNKoAAACMikYVAACAUTHrLwAAAMtg1l8AAADGw6y/AAAAHBqDGtWqOldVL1TVZlU9tMvxW6rqD+fHn6qq2xaeFAAAgCNhz0a1qm5K8kiSe5KcTXJ/VZ3dMeyBJK919/cm+fUkn1x0UAAAAI6GIXdU70qy2d1Xu/v1JI8mubhjzMUkvzd//bkkd1dVLS4mAAAAR8WxAWNuTfLitu1rSX7grcZ09xtVtZXku5K8un1QVV1Kcmm++c2q+tL1hD5Ax5NMabriqeVNppn5RHbU+shN8RrLvHxTq+Nketd4anmTaWaeWi1P8RrLvHxTq+Nketc4mV7mqeVNkjuqavsMuuvdvT7kjUMa1YWZh1pPkqrauN4ZoA5KVa1396W9R47D1PImk808qVqe6DWWecmmVsfJJK/xpPImk808qVqe6DWWecmmVsfJ9K5xMr3MU8ub3FgtD3n096Ukp7Ztn5zv23VMVR3LrNv/2vUEGrn/ftAB9mlqeZNpZp6aKV5jmdnN1K7x1PIm08w8NVO8xjKzmyle46llnlreG7LnOqrzxvMrSe7OrCF9Osn/293PbRvz40ne293/X1Xdl+T/6e6P7PG5k/tLEexGLXMYqGMOC7XMYaCOOSxupJb3fPR3/p3TB5M8nuSmJJ/q7ueq6uEkG919OcnvJvlMVW0m+fsk9w0496Bnk2EC1DKHgTrmsFDLHAbqmMPiumt5zzuqAAAAsEpDvqMKAAAAK6NRBQAAYFSW3qhW1bmqeqGqNqvqoV2O31JVfzg//lRV3bbsTLBfA+r4p6rq+ap6tqo+X1XfcxA5YS971fK2cT9UVV1VJvNgdIbUcVV9ZP5z+bmq+oNVZ4QhBvx+cbqqnqiqL85/x7j3IHLC26mqT1XVy1X1pbc4XlX1G/M6f7aq3jfkc5faqFbVTUkeSXJPkrNJ7q+qszuGPZDkte7+3iS/nuSTy8wE+zWwjr+YZK27vz/J55L88mpTwt4G1nKq6l1JfjLJU6tNCHsbUsdVdSbJzyb5QHf/6yT/YdU5YS8Dfyb/fJLHuvvOzCYr/c3VpoRBPp3k3NscvyfJmfm/S0l+a8iHLvuO6l1JNrv7ane/nuTRJBd3jLmY5Pfmrz+X5O6qqiXngv3Ys467+4nu/vp888nM1huGsRnyMzlJfimzPxp+Y5XhYKAhdfzxJI9092tJ0t0vrzgjDDGkljvJd8xfH0/ytyvMB4N09xcyW/nlrVxM8vs982SS76yq797rc5fdqN6a5MVt29fm+3Yd091vJNlK8l1LzgX7MaSOt3sgyZ8tNRFcnz1ref44zqnu/tNVBoN9GPIz+T1J3lNVf1VVT1bV2/2lHw7KkFr+xSQfraprSa4k+YnVRIOF2u/v0kkGrKMKDFdVH02yluSDB50F9quq3pHk15J87ICjwI06ltkjZh/K7AmXL1TVe7v7Hw4yFFyH+5N8urt/tap+MMlnquqO7v6ngw4Gy7bsO6ovJTm1bfvkfN+uY6rqWGaPNXxtyblgP4bUcarqw0l+LsmF7v7mirLBfuxVy+9KckeSv6yq/53k/Ukum1CJkRnyM/laksvd/Y/d/TdJvpJZ4wpjMqSWH0jyWJJ0918neWeSEytJB4sz6HfpnZbdqD6d5ExV3V5VN2f2JfDLO8ZcTvKj89c/nOQvuruXnAv2Y886rqo7k/xOZk2q70IxVm9by9291d0nuvu27r4ts+9bX+jujYOJC7sa8rvFH2d2NzVVdSKzR4GvrjAjDDGklr+a5O4kqarvy6xRfWWlKeHGXU7yI/PZf9+fZKu7/26vNy310d/ufqOqHkzyeJKbknyqu5+rqoeTbHT35SS/m9ljDJuZfQn3vmVmgv0aWMe/kuTbk/zRfC6wr3b3hQMLDbsYWMswagPr+PEk/66qnk/yf5L8THd7WotRGVjLP53kv1XVf8xsYqWPuaHD2FTVZzP74+CJ+fepfyHJtyVJd/92Zt+vvjfJZpKvJ/mxQZ+r1gEAABiTZT/6CwAAAPuiUQUAAGBUNKoAAACMikYVAACAUdGoAgAAMCp7NqpV9amqermqvvQWx6uqfqOqNqvq2ap63+JjAgAAcFQMuaP66STn3ub4PUnOzP9dSvJbQ05cVZeGjIOxU8scBuqYw0ItcxioYw6LG6nlPRvV7v5Ckr9/myEXk/x+zzyZ5Dur6rsHnHty/wGr6vxBZ9iPqeVNppk5E6vlKV5jmVdiUnWcTO8aTy1vMs3MmVgtT/Eay7wSk6rjZJLXeHKZp5Z3bnmN6gC3Jnlx2/a1+b7DaGrFMbW8yTQzT80Ur7HM7GZq13hqeZNpZp6aKV5jmdnNFK/x1DJPLe8Nqe7ee1DVbUn+pLvv2OXYnyT5L939P+bbn0/yn7p7Y5exlzLvqm+55ZZ/c8cdb/q4Udva2srx48cPOsZgU8ubTDPzK6+8kne/+90HHWOwKV5jmZdvanWcTO8aTy1vMs3MU6vlKV5jmZdvanWcTO8aJ9PLPLW8SfLMM898M8n2uY7Wu3t9yHuPLeD8LyU5tW375Hzfm8xDrSfJ2tpab2y8qZcFAADgEKiqL3X32vW8dxGP/l5O8iPz2X/fn2Sru/9uAZ8LAADAEbTnHdWq+mySDyU5UVXXkvxCkm9Lku7+7SRXktybZDPJ15P82LLCAgAAcPjt2ah29/17HO8kP76wRAAAABxpi3j0FwAAABZm5Y1qVZ2vqvWtra1VnxoAAIDVOV5V69ezBuyg5WmWway/AAAAh1dVPXOQs/4CAADAwmhUAQAAGBWNKgAAAKOiUQUAAGBUNKoAAACMiuVpAAAAWAbL0wAAADAelqcBAADg0NCoAgAAMCoaVQAAAEZFowoAAMComPUXAACAZTDrLwAAAONh1l8AAAAOjUGNalWdq6oXqmqzqh7a5fjpqnqiqr5YVc9W1b2LjwoAAMBRsGejWlU3JXkkyT1Jzia5v6rO7hj280ke6+47k9yX5DcXHRQAAICjYcgd1buSbHb31e5+PcmjSS7uGNNJvmP++niSv11cRAAAAI6SYwPG3JrkxW3b15L8wI4xv5jkz6vqJ5L8iyQf3u2DqupSkktJcvr06f1mBQAAYDpOVNX2GXTXu3t9yBuHNKpD3J/k0939q1X1g0k+U1V3dPc/bR80D7WezGb9XdC5AQAAGJ9Xlznr70tJTm3bPjnft90DSR5Lku7+6yTvTHLiegIBAABwtA1pVJ9Ocqaqbq+qmzObLOnyjjFfTXJ3klTV92XWqL6yyKAAAAAcDXs2qt39RpIHkzye5MuZze77XFU9XFUX5sN+OsnHq+p/Jflsko91t0d7AQAA2LdB31Ht7itJruzY94ltr59P8oHFRgMAAOAoGvLoLwAAAKzMyhvVqjpfVetbW1urPjUAAACrc7yq1qvq/H7fWAf1VdK1tbXe2NjYeyAAAACTU1XPLHN5GgAAAFgZjSoAAACjolEFAABgVDSqAAAAjIpZfwEAAFgGs/4CAAAwHmb9BQAA4NDQqAIAADAqGlUAAABGRaMKAADAqGhUAQAAGBXL0wAAALAMlqcBAABgPJa+PE1VnauqF6pqs6oeeosxH6mq56vquar6g+sJAwAAAMf2GlBVNyV5JMm/TXItydNVdbm7n9825kySn03yge5+rar+5bICAwAAcLgNuaN6V5LN7r7a3a8neTTJxR1jPp7kke5+LUm6++XFxgQAAOCoGNKo3prkxW3b1+b7tntPkvdU1V9V1ZNVdW5RAQEAADha9nz0dx+fcybJh5KcTPKFqnpvd//D9kFVdSnJpSQ5ffr0gk4NAADACJ2oqu0z6K539/qQNw5pVF9Kcmrb9sn5vu2uJXmqu/8xyd9U1Vcya1yf3j5oHmo9mc36OyQgAAAAk/TqMmf9fTrJmaq6vapuTnJfkss7xvxxZndTU1UnMnsU+Or1BAIAAOBo27NR7e43kjyY5PEkX07yWHc/V1UPV9WF+bDHk3ytqp5P8kSSn+nury0rNAAAAIdXdR/ME7hra2u9sbGx90AAAAAmp6qeWeajvwAAALAyK29Uq+p8Va1vbW2t+tQAAACszvGqWq+q8/t9o0d/AQAAWDiP/gIAAHBoaFQBAAAYFY0qAAAAo6JRBQAAYFQ0qgAAAIyK5WkAAABYBsvTAAAAMB6WpwEAAODQ0KgCAAAwKhpVAAAARkWjCgAAwKiY9RcAAIBlMOsvAAAA42HWXwAAAA6NQY1qVZ2rqheqarOqHnqbcT9UVV1V19U1AwAAwJ6NalXdlOSRJPckOZvk/qo6u8u4dyX5ySRPLTokAAAAR8eQO6p3Jdns7qvd/XqSR5Nc3GXcLyX5ZJJvLDAfAAAAR8yQRvXWJC9u27423/d/VdX7kpzq7j9dYDYAAACOoGM3+gFV9Y4kv5bkYwPGXkpyKUlOnz59o6cGAABgvE5U1falXta7e33IG4c0qi8lObVt++R837e8K8kdSf6yqpLkXyW5XFUXuvufrT8zD7WezJanGRIQAACASXp1mcvTPJ3kTFXdXlU3J7kvyeVvHezure4+0d23dfdtSZ5M8qYmFQAAAIbYs1Ht7jeSPJjk8SRfTvJYdz9XVQ9X1YVlBwQAAOBoGfQd1e6+kuTKjn2feIuxH7rxWAAAABxVQx79BQAAgJVZeaNaVeeran1ra2vVpwYAAGB1jlfVelWd3+8bq/tgJt9dW1vrjQ3zLQEAABxGVfXMMmf9BQAAgJXRqAIAADAqGlUAAABGRaMKAADAqJj1FwAAgGUw6y8AAADjYdZfAAAADg2NKgAAAKOiUQUAAGBUNKoAAACMikYVAACAUbE8DQAAAMtgeRoAAADGY+nL01TVuap6oao2q+qhXY7/VFU9X1XPVtXnq+p7ricMAAAA7NmoVtVNSR5Jck+Ss0nur6qzO4Z9Mclad39/ks8l+eVFBwUAAOBoGHJH9a4km919tbtfT/JokovbB3T3E9399fnmk0lOLjYmAAAAR8WQRvXWJC9u27423/dWHkjyZzcSCgAAgKPr2CI/rKo+mmQtyQff4vilJJeS5PTp04s8NQAAAONyoqq2z6C73t3rQ944pFF9Kcmpbdsn5/v+mar6cJKfS/LB7v7mbh80D7WezGb9HRIQAACASXp1mbP+Pp3kTFXdXlU3J7kvyeXtA6rqziS/k+RCd798PUEAAAAgGdCodvcbSR5M8niSLyd5rLufq6qHq+rCfNivJPn2JH9UVf+zqi6/xccBAADA2xr0HdXuvpLkyo59n9j2+sMLzgUAAMARNeTRXwAAAFgZjSoAAACjsvJGtarOV9X61tbWqk8NAADA6hyvqvWqOr/fN1b3wawSs7a21hsbG3sPBAAAYHKq6pllLk8DAAAAK6NRBQAAYFQ0qgAAAIyKRhUAAIBRMesvAAAAy2DWXwAAAMbDrL8AAAAcGhpVAAAARkWjCgAAwKhoVAEAABgVs/4CAACwDGb9BQAAYDzM+gsAAMChMahRrapzVfVCVW1W1UO7HL+lqv5wfvypqrpt4UkBAAA4EvZsVKvqpiSPJLknydkk91fV2R3DHkjyWnd/b5JfT/LJRQcFAADgaBhyR/WuJJvdfbW7X0/yaJKLO8ZcTPJ789efS3J3VdXiYgIAAHBUDGlUb03y4rbta/N9u47p7jeSbCX5rkUEBAAA4Gg5tsqTVdWlJJfmm9+sqi+t8vwLcDyzJnwqppY3mWbmE0lePegQ+zDFayzz8k2tjpPpXeOp5U2mmXlqtTzFayzz8k2tjpPpXeNkepmnljdJ7qiq7Uu9rHf3+pA3DmlUX0pyatv2yfm+3cZcq6pjmV3Er+38oHmo9SSpqo3rnar4oFTVendf2nvkOEwtbzLZzJOq5YleY5mXbGp1nEzyGk8qbzLZzJOq5YleY5mXbGp1nEzvGifTyzy1vMmN1fKQR3+fTnKmqm6vqpuT3Jfk8o4xl5P86Pz1Dyf5iz6oBVqX678fdIB9mlreZJqZp2aK11hmdjO1azy1vMk0M0/NFK+xzOxmitd4apmnlveG1JB+sqruTfJfk9yU5FPd/Z+r6uEkG919uaremeQzSe5M8vdJ7uvuq3t85uT+UgS7UcscBuqYw0ItcxioYw6LG6nlQd9R7e4rSa7s2PeJba+/keTf7/Pcg55NhglQyxwG6pjDQi1zGKhjDovrruVBd1QBAABgVYZ8RxUAAABWZumNalWdq6oXqmqzqh7a5fgtVfWH8+NPVdVty84E+zWgjn+qqp6vqmer6vNV9T0HkRP2slctbxv3Q1XVVeU7UozOkDquqo/Mfy4/V1V/sOqMMMSA3y9OV9UTVfXF+e8Y9x5ETng7VfWpqnr5rZYerZnfmNf5s1X1viGfu9RGtapuSvJIknuSnE1yf1Wd3THsgSSvdff3Jvn1JJ9cZibYr4F1/MUka939/Uk+l+SXV5sS9jawllNV70ryk0meWm1C2NuQOq6qM0l+NskHuvtfJ/kPq84Jexn4M/nnkzzW3XdmtvLGb642JQzy6STn3ub4PUnOzP9dSvJbQz502XdU70qy2d1Xu/v1JI8mubhjzMUkvzd//bkkd1dVLTkX7MeeddzdT3T31+ebT2a23jCMzZCfyUnyS5n90fAbqwwHAw2p448neaS7X0uS7n55xRlhiCG13Em+Y/76eJK/XWE+GKS7v5DZyi9v5WKS3++ZJ5N8Z1V9916fu+xG9dYkL27bvjbft+uY7n4jyVaS71pyLtiPIXW83QNJ/mypieD67FnL88dxTnX3n64yGOzDkJ/J70nynqr6q6p6sqre7i/9cFCG1PIvJvloVV3LbAWOn1hNNFio/f4unWTg8jTAMFX10SRrST540Flgv6rqHUl+LcnHDjgK3KhjmT1i9qHMnnD5QlW9t7v/4SBDwXW4P8mnu/tXq+oHk3ymqu7o7n866GCwbMu+o/pSklPbtk/O9+06pqqOZfZYw9eWnAv2Y0gdp6o+nOTnklzo7m+uKBvsx161/K4kdyT5y6r630nen+SyCZUYmSE/k68ludzd/9jdf5PkK5k1rjAmQ2r5gSSPJUl3/3WSdyY5sZJ0sDiDfpfeadmN6tNJzlTV7VV1c2ZfAr+8Y8zlJD86f/3DSf6iLe7KuOxZx1V1Z5LfyaxJ9V0oxupta7m7t7r7RHff1t23ZfZ96wvdvXEwcWFXQ363+OPM7qamqk5k9ijw1RVmhCGG1PJXk9ydJFX1fZk1qq+sNCXcuMtJfmQ+++/7k2x199/t9aalPvrb3W9U1YNJHk9yU5JPdfdzVfVwko3uvpzkdzN7jGEzsy/h3rfMTLBfA+v4V5J8e5I/ms8F9tXuvnBgoWEXA2sZRm1gHT+e5N9V1fNJ/k+Sn+luT2sxKgNr+aeT/Leq+o+ZTaz0MTd0GJuq+mxmfxw8Mf8+9S8k+bYk6e7fzuz71fcm2Uzy9SQ/Nuhz1ToAAABjsuxHfwEAAGBfNKoAAACMikYVAACAUdGoAgAAMCoaVQAAAEZFowoAAMCoaFQBAAAYFY0qAAAAo/L/AzxNZIvFjamgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1, axs = plt.subplots(3,1,figsize=(16,8))\n",
    "\n",
    "ii = 7000\n",
    "\n",
    "axs[0].scatter(t_set_plt,noise,s=5)\n",
    "axs[0].set_xlim(raster1_start_mjd[0],raster1_stop_mjd[-1])\n",
    "axs2 = axs[0].twinx()\n",
    "axs2.plot(t_set_plt,dec,color='k')\n",
    "axs[0].scatter(t_noise,LL_dnoise_arr[:,ii]/np.nanmean(LL_dnoise_arr[:,ii]),s=10,c='red')\n",
    "#axs2.set_xlim(raster1_start_mjd,raster1_stop_mjd)\n",
    "axs2.set_ylim(5,65)\n",
    "axs2.set_ylabel('dec.')\n",
    "\n",
    "axs[1].scatter(t_set_plt,noise,s=5)\n",
    "axs[1].set_xlim(scan_start_mjd[0],scan_stop_mjd[-1])\n",
    "axs2 = axs[1].twinx()\n",
    "axs2.scatter(t_set_plt,az,color='k',s=0.01)\n",
    "axs[1].scatter(t_noise,LL_dnoise_arr[:,ii]/np.nanmean(LL_dnoise_arr[:,ii]),s=10,c='red')\n",
    "#axs2.set_xlim(scan_start_mjd[0],scan_stop_mjd[-1])\n",
    "axs2.set_ylim(-20,380)\n",
    "axs2.set_ylabel('az.')\n",
    "\n",
    "axs[2].scatter(t_set_plt,noise,s=5)\n",
    "axs[2].set_xlim(raster2_start_mjd[0],raster2_stop_mjd[-1])\n",
    "axs2 = axs[2].twinx()\n",
    "axs2.plot(t_set_plt,dec,color='k')\n",
    "axs[2].scatter(t_noise,LL_dnoise_arr[:,ii]/np.nanmean(LL_dnoise_arr[:,ii]),s=10,c='red')\n",
    "#axs2.set_xlim(raster1_start_mjd,raster1_stop_mjd)\n",
    "axs2.set_ylim(5,65)\n",
    "axs2.set_ylabel('dec.')\n",
    "\n",
    "for i in range(0,3):\n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].set_ylim(-0.1,1.2)\n",
    "    axs[i].set_ylabel('noise (arbitrary units)')\n",
    "axs[2].set_xlabel('UTC')\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/noise_firing_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44242e98",
   "metadata": {},
   "source": [
    "## Check noise deflection vs frequency for sample noise source instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d51c7ba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, axs = plt.subplots(2,1,figsize=(16,8))\n",
    "\n",
    "t_noise_idx = 100\n",
    "\n",
    "axs[0].scatter(freq,10*np.log10(LL_dnoise_arr[t_noise_idx,:]),s=1,label='LL')\n",
    "axs[0].scatter(freq,10*np.log10(RR_dnoise_arr[t_noise_idx,:]),s=1,label='RR')\n",
    "axs[0].scatter(freq,10*np.log10(np.sqrt((reRL_dnoise_arr[t_noise_idx,:])**2 + \n",
    "                                        (imRL_dnoise_arr[t_noise_idx,:])**2)),s=1,label='PI')\n",
    "axs[0].grid()\n",
    "axs[0].set_ylim(55,75)\n",
    "axs[0].set_xlim(350,1030)\n",
    "axs[0].set_ylabel('Power deflection (dB)')\n",
    "axs[0].legend(markerscale=10)\n",
    "\n",
    "axs[1].scatter(freq,reRL_dnoise_arr[t_noise_idx,:],s=1,label='reRL')\n",
    "axs[1].scatter(freq,imRL_dnoise_arr[t_noise_idx,:],s=1,label='imRL')\n",
    "axs[1].grid()\n",
    "axs[1].set_ylim(-2e7,2e7)\n",
    "axs[1].set_xlim(350,1030)\n",
    "axs[1].set_xlabel('Frequency (MHz)')\n",
    "axs[1].set_ylabel('Power deflection (spectrometer units)')\n",
    "axs[1].legend(markerscale=10)\n",
    "\n",
    "plt.tight_layout()\n",
    "#plt.savefig('../DVA_PLOTS/noise_gain_corr/noise_power_deflection_vs_freq_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6dc60c7",
   "metadata": {},
   "source": [
    "## Remove outliers in noise deflection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02c44fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "LL_dnoise_fix = LL_dnoise_arr.copy()\n",
    "RR_dnoise_fix = RR_dnoise_arr.copy()\n",
    "reRL_dnoise_fix = reRL_dnoise_arr.copy()\n",
    "imRL_dnoise_fix = imRL_dnoise_arr.copy()\n",
    "std_mult = 1.0\n",
    "for i in range(0,len(freq)):\n",
    "    \n",
    "    print(i,freq[i])    \n",
    "    for j in range(0,len(t_noise)):\n",
    "        \n",
    "        window = 10\n",
    "        if ((j >= window) & (j < len(t_noise)-window)):\n",
    "            LL_near = np.nanmedian(LL_dnoise_arr[j-window:j+window,i])\n",
    "            RR_near = np.nanmedian(RR_dnoise_arr[j-window:j+window,i])\n",
    "        else:\n",
    "            if j < window:\n",
    "                LL_near = np.nanmedian(LL_dnoise_arr[0:j+window,i])\n",
    "                RR_near = np.nanmedian(RR_dnoise_arr[0:j+window,i])\n",
    "            if j >= len(t_noise)-window:\n",
    "                LL_near = np.nanmedian(LL_dnoise_arr[j-window:-1,i])\n",
    "                RR_near = np.nanmedian(RR_dnoise_arr[j-window:-1,i])\n",
    "        \n",
    "        if ( LL_dnoise_arr[j,i] > 1.1*LL_near):\n",
    "            LL_dnoise_fix[j,i] = np.nan\n",
    "            reRL_dnoise_fix[j,i] = np.nan\n",
    "            imRL_dnoise_fix[j,i] = np.nan\n",
    "        if ( LL_dnoise_arr[j,i] < 0.9*LL_near):\n",
    "            LL_dnoise_fix[j,i] = np.nan\n",
    "            reRL_dnoise_fix[j,i] = np.nan\n",
    "            imRL_dnoise_fix[j,i] = np.nan\n",
    "            \n",
    "        if ( RR_dnoise_arr[j,i] > 1.1*RR_near):\n",
    "            RR_dnoise_fix[j,i] = np.nan\n",
    "            reRL_dnoise_fix[j,i] = np.nan\n",
    "            imRL_dnoise_fix[j,i] = np.nan\n",
    "        if ( RR_dnoise_arr[j,i] < 0.9*RR_near):\n",
    "            RR_dnoise_fix[j,i] = np.nan\n",
    "            reRL_dnoise_fix[j,i] = np.nan\n",
    "            imRL_dnoise_fix[j,i] = np.nan\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dd21b50",
   "metadata": {},
   "source": [
    "## Calculate fits to noise deflection vs time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49078f5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "LL_fit = np.empty_like(LL)\n",
    "RR_fit = np.empty_like(RR)\n",
    "reRL_fit = np.empty_like(reRL)\n",
    "imRL_fit = np.empty_like(imRL)\n",
    "\n",
    "LL_fit_noise_pt = np.empty_like(LL_dnoise_arr)\n",
    "RR_fit_noise_pt = np.empty_like(RR_dnoise_arr)\n",
    "reRL_fit_noise_pt = np.empty_like(reRL_dnoise_arr)\n",
    "imRL_fit_noise_pt = np.empty_like(imRL_dnoise_arr)\n",
    "\n",
    "for i in range(0,len(freq)):\n",
    "    print(freq[i])\n",
    "    try:\n",
    "        wuse = np.where(np.isfinite(LL_dnoise_fix[:,i]))[0]\n",
    "        pL = np.polyfit(np.array(t_noise)[wuse]-t_mean,LL_dnoise_fix[wuse,i],9)\n",
    "        LL_fit[:,i] = np.polyval(pL, t_set_plt-t_mean)\n",
    "        LL_fit_noise_pt[:,i] = np.polyval(pL, t_noise-t_mean)\n",
    "\n",
    "        wuse = np.where(np.isfinite(RR_dnoise_fix[:,i]))[0]\n",
    "        pR = np.polyfit(np.array(t_noise)[wuse]-t_mean,RR_dnoise_fix[wuse,i],9)\n",
    "        RR_fit[:,i] = np.polyval(pR, t_set_plt-t_mean)\n",
    "        RR_fit_noise_pt[:,i] = np.polyval(pR, t_noise-t_mean)\n",
    "    \n",
    "        wuse = np.where(np.isfinite(reRL_dnoise_fix[:,i]))[0]\n",
    "        pre = np.polyfit(np.array(t_noise)[wuse]-t_mean,reRL_dnoise_fix[wuse,i],9)\n",
    "        reRL_fit[:,i] = np.polyval(pre, t_set_plt-t_mean)\n",
    "        reRL_fit_noise_pt[:,i] = np.polyval(pre, t_noise-t_mean)\n",
    "    \n",
    "        wuse = np.where(np.isfinite(imRL_dnoise_fix[:,i]))[0]\n",
    "        pim = np.polyfit(np.array(t_noise)[wuse]-t_mean,imRL_dnoise_fix[wuse,i],9)\n",
    "        imRL_fit[:,i] = np.polyval(pim, t_set_plt-t_mean)\n",
    "        imRL_fit_noise_pt[:,i] = np.polyval(pim, t_noise-t_mean)\n",
    "    except:\n",
    "        pass\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5052978f",
   "metadata": {},
   "source": [
    "## Read in persistent RFI mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e841ed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0\n",
    "RFI_mask_idx = []\n",
    "with open('/srv/data/dva/RFIpersist_mask/RFIpersist_mask.txt') as fp:\n",
    "    for line in fp:\n",
    "        if i>0: \n",
    "            #print(line)\n",
    "            RFI_mask_idx.append(int(line.split()[0]))\n",
    "        i=i+1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b10b7a62",
   "metadata": {},
   "source": [
    "## Plot power and gain for a single frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d24b995e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fplot = 800\n",
    "df = freq[1]-freq[0]\n",
    "wf = np.where(abs(freq-fplot)<df/2)[0][0]\n",
    "print(freq[wf])\n",
    "\n",
    "fig1, axs = plt.subplots(3,1,figsize=(16,12))\n",
    "\n",
    "axs[0].scatter(t_noise,10*np.log10(RR_dnoise_arr[:,wf]),s=20,color='red')\n",
    "axs[0].scatter(t_noise,10*np.log10(RR_dnoise_fix[:,wf]),s=20,color='C1',label='RR',zorder=10)\n",
    "axs[0].scatter(t_noise,10*np.log10(LL_dnoise_arr[:,wf]),s=20,color='blue')\n",
    "axs[0].scatter(t_noise,10*np.log10(LL_dnoise_fix[:,wf]),s=20,color='C0',label='LL',zorder=10)\n",
    "axs[0].plot(t_set_plt,10*np.log10(LL_fit[:,wf]),color='blue',linewidth=2,zorder=0)\n",
    "axs[0].plot(t_set_plt,10*np.log10(RR_fit[:,wf]),color='red',linewidth=2,zorder=0)\n",
    "axs[0].set_ylim(65,69)\n",
    "axs[0].set_ylabel('Noise source deflection (dB)')\n",
    "axs[0].legend(markerscale=2)\n",
    "\n",
    "ax2 = axs[0].twinx()\n",
    "ax2.plot(t_weath_plt,temp_C,color='purple',linewidth=2)\n",
    "ax2.set_ylim(-5,35)\n",
    "ax2.set_ylabel('Outside Temperature')\n",
    "\n",
    "axs[1].scatter(t_noise,(RR_dnoise_fix[:,wf]-RR_fit_noise_pt[:,wf])/RR_dnoise_fix[:,wf],\n",
    "               color='C1',zorder=10)\n",
    "axs[1].set_ylim(-0.15,0.15)\n",
    "axs[1].fill_between(t_set_plt,-0.05,0.05,zorder=0,color='grey',alpha=0.4)\n",
    "axs[1].fill_between(t_set_plt,-0.01,0.01,zorder=1,color='grey',alpha=0.4)\n",
    "axs[1].set_ylabel('Fractional residuals')\n",
    "\n",
    "axs[2].scatter(t_noise,(LL_dnoise_fix[:,wf]-LL_fit_noise_pt[:,wf])/LL_dnoise_fix[:,wf],\n",
    "               color='C0',zorder=10)\n",
    "axs[2].set_ylim(-0.15,0.15)\n",
    "axs[2].fill_between(t_set_plt,-0.05,0.05,zorder=0,color='grey',alpha=0.4)\n",
    "axs[2].fill_between(t_set_plt,-0.01,0.01,zorder=1,color='grey',alpha=0.4)\n",
    "axs[2].set_ylabel('Fractional residuals')\n",
    "\n",
    "\n",
    "for i in range(0,3):\n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].set_xlim(raster1_start_mjd[0],raster2_stop_mjd[-1])\n",
    "    axs[i].grid()\n",
    "axs[2].set_xlabel('UTC')\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/noise_deflection_fits_800MHz_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "765f489b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fplot = 800\n",
    "df = freq[1]-freq[0]\n",
    "wf = np.where(abs(freq-fplot)<df/2)[0][0]\n",
    "print(freq[wf])\n",
    "\n",
    "fig1, axs = plt.subplots(3,1,figsize=(16,12))\n",
    "\n",
    "axs[0].scatter(t_noise,imRL_dnoise_arr[:,wf],s=20,color='red')\n",
    "axs[0].scatter(t_noise,imRL_dnoise_fix[:,wf],s=20,color='C1',label='imRL',zorder=10)\n",
    "axs[0].scatter(t_noise,reRL_dnoise_arr[:,wf],s=20,color='blue')\n",
    "axs[0].scatter(t_noise,reRL_dnoise_fix[:,wf],s=20,color='C0',label='reRL',zorder=10)\n",
    "axs[0].plot(t_set_plt,reRL_fit[:,wf],color='blue',linewidth=2,zorder=0)\n",
    "axs[0].plot(t_set_plt,imRL_fit[:,wf],color='red',linewidth=2,zorder=0)\n",
    "#axs[0].set_ylim(65,69)\n",
    "axs[0].set_ylabel('Noise source deflection (raw spectrometer units)')\n",
    "axs[0].legend(markerscale=2)\n",
    "\n",
    "ax2 = axs[0].twinx()\n",
    "ax2.plot(t_weath_plt,temp_C,color='purple',linewidth=2)\n",
    "ax2.set_ylim(-5,35)\n",
    "ax2.set_ylabel('Outside Temperature')\n",
    "\n",
    "axs[1].scatter(t_noise,(imRL_dnoise_fix[:,wf]-imRL_fit_noise_pt[:,wf])/imRL_dnoise_fix[:,wf],\n",
    "               color='C1',zorder=10)\n",
    "axs[1].set_ylim(-0.15,0.15)\n",
    "axs[1].fill_between(t_set_plt,-0.05,0.05,zorder=0,color='grey',alpha=0.4)\n",
    "axs[1].fill_between(t_set_plt,-0.01,0.01,zorder=1,color='grey',alpha=0.4)\n",
    "axs[1].set_ylabel('Fractional residuals')\n",
    "\n",
    "axs[2].scatter(t_noise,(reRL_dnoise_fix[:,wf]-reRL_fit_noise_pt[:,wf])/reRL_dnoise_fix[:,wf],\n",
    "               color='C0',zorder=10)\n",
    "axs[2].set_ylim(-0.15,0.15)\n",
    "axs[2].fill_between(t_set_plt,-0.05,0.05,zorder=0,color='grey',alpha=0.4)\n",
    "axs[2].fill_between(t_set_plt,-0.01,0.01,zorder=1,color='grey',alpha=0.4)\n",
    "axs[2].set_ylabel('Fractional residuals')\n",
    "\n",
    "\n",
    "for i in range(0,3):\n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].set_xlim(raster1_start_mjd[0],raster2_stop_mjd[-1])\n",
    "    axs[i].grid()\n",
    "axs[2].set_xlabel('UTC')\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/noise_deflection_pol_fits_800MHz_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e697571",
   "metadata": {},
   "outputs": [],
   "source": [
    "fplot = 800\n",
    "df = freq[1]-freq[0]\n",
    "wf = np.where(abs(freq-fplot)<df/2)[0][0]\n",
    "print(freq[wf])\n",
    "\n",
    "\n",
    "lims_raw1 = [65,65,65]\n",
    "lims_raw2 = [75,75,75]\n",
    "lims_cor1 = [2,2,2]\n",
    "lims_cor2 = [8,8,8]\n",
    "\n",
    "sz=1\n",
    "fig1, axs = plt.subplots(3,1,figsize=(20,12))\n",
    "\n",
    "for i in range(0,3):\n",
    "    axs[i].scatter(t_set_plt,10*np.log10(LL[:,wf]),s=sz)\n",
    "    axs[i].set_ylim(lims_raw1[i],lims_raw2[i])\n",
    "        \n",
    "    ax2 = axs[i].twinx()\n",
    "    ax2.scatter(t_set_plt,LL[:,wf]/LL_fit[:,wf],s=sz,color='red')\n",
    "    ax2.set_ylim(lims_cor1[i],lims_cor2[i])\n",
    "    ax2.set_ylabel('Corrected power (noise src units; red)')\n",
    "    \n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].grid()\n",
    "    axs[i].set_ylabel('Raw power (dB; blue)')\n",
    "\n",
    "axs[0].set_xlim(raster1_start_mjd[0],raster1_stop_mjd[-1])\n",
    "#axs[0].set_xlim(raster1_start_mjd,raster1_start_mjd+0.005)\n",
    "#axs[0].set_ylim(2e7,4e7)\n",
    "axs[1].set_xlim(scan_start_mjd[0],scan_stop_mjd[-1])\n",
    "#axs[1].set_ylim(1e7,10e7)\n",
    "axs[2].set_xlim(raster2_start_mjd[0],raster2_stop_mjd[-1])\n",
    "#axs[2].set_ylim(1e7,10e7)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/correction_example_LL_800MHz_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f0affac",
   "metadata": {},
   "source": [
    "## Checking frequency dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a38482",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, axs = plt.subplots(2,1,figsize=(16,12))\n",
    "\n",
    "t_idx = 20000\n",
    "print(t[t_idx])\n",
    "s = 0.5\n",
    "\n",
    "axs[0].scatter(freq,10*np.log10(LL[t_idx,:]),s=s,color='C0',label='LL')\n",
    "axs[0].scatter(freq,10*np.log10(RR[t_idx,:]),s=s,color='C1',label='RR')\n",
    "ax2 = axs[0].twinx()\n",
    "ax2.scatter(freq,10*np.log10(np.sqrt((reRL[t_idx,:])**2+(imRL[t_idx])**2)),s=s,color='C2',label='PI')\n",
    "axs[0].grid()\n",
    "axs[0].set_ylim(65,80)\n",
    "ax2.set_ylim(20,80)\n",
    "axs[0].set_xlim(350,1030)\n",
    "axs[0].set_ylabel('Raw power (dB)')\n",
    "axs[0].legend(markerscale=10,loc='upper right')\n",
    "ax2.legend(markerscale=10,loc='lower right')\n",
    "\n",
    "axs[1].scatter(freq,LL[t_idx,:]/LL_fit[t_idx,:],s=s,color='blue',label='LL')\n",
    "axs[1].scatter(freq,RR[t_idx,:]/RR_fit[t_idx,:],s=s,color='red',label='RR')\n",
    "ax2 = axs[1].twinx()\n",
    "ax2.scatter(freq,np.sqrt((reRL[t_idx,:]/reRL_fit[t_idx,:])**2+(imRL[t_idx]/imRL_fit[t_idx])**2),s=s,color='green',label='PI')\n",
    "ax2.set_ylim(-1,2)\n",
    "axs[1].set_ylim(0,10)\n",
    "axs[1].grid()\n",
    "axs[1].set_xlim(350,1030)\n",
    "axs[1].set_ylabel('Corrected power (noise source units)')\n",
    "axs[1].set_xlabel('Frequency (MHz)')\n",
    "axs[1].legend(markerscale=10,loc='upper right')\n",
    "ax2.legend(markerscale=10,loc='lower right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/correction_example_spectrum_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8549418",
   "metadata": {},
   "source": [
    "## Write out gain-corrected files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7014ec7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_gain_corrected(file,RR_out,LL_out,reRL_out,imRL_out,infiles,outfiles,inname,outname,GR,GL,GreRL,GimRL):\n",
    "\n",
    "    cmd2 = 'cp '+infiles+inname+' '+outfiles+outname+'.h5'\n",
    "    print(cmd2)\n",
    "    os.system(cmd2)\n",
    "    print('got here')\n",
    "    file_new = h5py.File(outfiles+outname+'.h5','r+')\n",
    "    \n",
    "    print(inname)\n",
    "    print(RR_out.shape)\n",
    "    print(GR.shape)\n",
    "    print(file_new['data']['beam_0']['band_SB0']['scan_0']['data'][:,0,:].shape)\n",
    "    print('')\n",
    "    \n",
    "    file_new['data']['beam_0']['band_SB0']['scan_0']['data'][:,0,:] = RR_out/GR\n",
    "    file_new['data']['beam_0']['band_SB0']['scan_0']['data'][:,1,:] = LL_out/GL\n",
    "    file_new['data']['beam_0']['band_SB0']['scan_0']['data'][:,2,:] = reRL_out/GreRL\n",
    "    file_new['data']['beam_0']['band_SB0']['scan_0']['data'][:,3,:] = imRL_out/GimRL\n",
    "           \n",
    "    file_new.close()\n",
    "\n",
    "    \n",
    "# Raster scan 1:\n",
    "for i in range(0,Ne):   \n",
    "    if Ne == 1:\n",
    "        inname  = 'dva_survey_phase'+phase+'_day0'+day+'_raster1'+'.h5'\n",
    "        outname = 'dva_survey_phase'+phase+'_day0'+day+'_raster1_gain_corr'\n",
    "    else:\n",
    "        inname  = 'dva_survey_phase'+phase+'_day0'+day+'_raster1'+rast_list[i]+'.h5'\n",
    "        outname = 'dva_survey_phase'+phase+'_day0'+day+'_raster1'+rast_list[i]+'_gain_corr'\n",
    "        \n",
    "    file = h5py.File(dir_in_rast+inname,'r')\n",
    "\n",
    "    write_gain_corrected(file,RR[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                              LL[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                            reRL[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                            imRL[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                         dir_in_rast,dir_out_rast,inname,outname,\n",
    "                         RR_fit[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                         LL_fit[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                       reRL_fit[idx1_rast1[i]:idx2_rast1[i],:],\n",
    "                       imRL_fit[idx1_rast1[i]:idx2_rast1[i],:])\n",
    "\n",
    "# Raster scan 2:\n",
    "for i in range(0,Nm):   \n",
    "    if Nm == 1:\n",
    "        inname  = 'dva_survey_phase'+phase+'_day0'+day+'_raster2'+'.h5'\n",
    "        outname = 'dva_survey_phase'+phase+'_day0'+day+'_raster2_gain_corr'\n",
    "    else:\n",
    "        inname  = 'dva_survey_phase'+phase+'_day0'+day+'_raster2'+rast_list[i]+'.h5'\n",
    "        outname = 'dva_survey_phase'+phase+'_day0'+day+'_raster2'+rast_list[i]+'_gain_corr'\n",
    "        \n",
    "    file = h5py.File(dir_in_rast+inname,'r')\n",
    "\n",
    "    write_gain_corrected(file,RR[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                              LL[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                            reRL[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                            imRL[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                         dir_in_rast,dir_out_rast,inname,outname,\n",
    "                         RR_fit[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                         LL_fit[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                       reRL_fit[idx1_rast2[i]:idx2_rast2[i],:],\n",
    "                       imRL_fit[idx1_rast2[i]:idx2_rast2[i],:])\n",
    "\n",
    "\n",
    "\n",
    "#file = h5py.File(dir_in_rast+inname,'r')\n",
    "#w = np.where((t_set_plt >= rast2_file_t[0]) & (t_set_plt <= rast2_file_t[-1]))[0]\n",
    "#write_gain_corrected(file,RR[w,:],LL[w,:],reRL[w,:],imRL[w,:],dir_in_rast,dir_out_rast,inname,outname,\n",
    "#                     GL_fit[w],GR_fit[w])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d9da19",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in scan_id:\n",
    "    inname  = 'dva_survey_phase'+phase+'_raw_'+f\"{i:04}\"+'.h5'\n",
    "    outname = 'dva_survey_phase'+phase+'_raw_'+f\"{i:04}\"+'_gain_corr'\n",
    "    file = h5py.File(dir_in_az+inname,'r')\n",
    "    print(i,file)\n",
    "    t_az = Time(file['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][:],\n",
    "                        format='isot',scale='utc').mjd\n",
    "    w = np.where((t_set_plt >= t_az[0]) & (t_set_plt <= t_az[-1]))[0]\n",
    "    print(w.shape)\n",
    "    print('')\n",
    "    write_gain_corrected(file,RR[w,:],LL[w,:],reRL[w,:],imRL[w,:],\n",
    "                         dir_in_az,dir_out_az,inname,outname,\n",
    "                         RR_fit[w,:],LL_fit[w,:],reRL_fit[w,:],imRL_fit[w,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57e51f11",
   "metadata": {},
   "outputs": [],
   "source": [
    "del RR\n",
    "del reRL\n",
    "del imRL\n",
    "del RR_fit\n",
    "del reRL_fit\n",
    "del imRL_fit\n",
    "\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c19fbf05",
   "metadata": {},
   "source": [
    "## Waterfall plot for power, gain and temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43a7845a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#Gmin = 0.9\n",
    "#Gmax = 1.1\n",
    "noisemin = 55\n",
    "noisemax = 75\n",
    "dBmin = 65\n",
    "dBmax = 75\n",
    "\n",
    "fplot = 800\n",
    "df = freq[1]-freq[0]\n",
    "wf = np.where(abs(freq-fplot)<df/2)[0][0]\n",
    "print(freq[wf])\n",
    "\n",
    "fig1, axs = plt.subplots(4,1,figsize=(20,25),sharex=True,gridspec_kw={'height_ratios': [1,0.5,1,1]},\n",
    "                         constrained_layout=True)\n",
    "fs = 16\n",
    "\n",
    "data1 = LL.T\n",
    "data2 = LL_fit.T\n",
    "data1[RFI_mask_idx,:] = np.nan\n",
    "data2[RFI_mask_idx,:] = np.nan\n",
    "\n",
    "\n",
    "for i in range(0,len(scan_id)):\n",
    "    w = np.where((t_set_plt>=scan_start_mjd[i]) & (t_set_plt<=scan_stop_mjd[i]))[0]\n",
    "    extent = [scan_start_mjd[i],scan_stop_mjd[i],freq[0],freq[-1]]\n",
    "    \n",
    "    im2 = axs[2].imshow(10*np.log10(data1[:,w]),aspect='auto',vmin=dBmin,vmax=dBmax,origin='lower',\n",
    "                        extent=extent,cmap='viridis')    \n",
    "    im0 = axs[0].imshow(10*np.log10(data2[:,w]),aspect='auto',vmin=noisemin,vmax=noisemax,origin='lower',\n",
    "                        extent=extent,cmap='viridis')\n",
    "    im3 = axs[3].imshow(data1[:,w]/data2[:,w],aspect='auto',vmin=1,vmax=8,origin='lower',\n",
    "                        extent=extent,cmap='viridis')\n",
    "\n",
    "\n",
    "for i in range(0,Ne):\n",
    "    w = np.where((t_set_plt>=raster1_start_mjd[i]) & (t_set_plt<=raster1_stop_mjd[i]))[0]\n",
    "    extent = [raster1_start_mjd[i],raster1_stop_mjd[i],freq[0],freq[-1]]\n",
    "\n",
    "    im2 = axs[2].imshow(10*np.log10(data1[:,w]),aspect='auto',vmin=dBmin,vmax=dBmax,origin='lower',\n",
    "                        extent=extent,cmap='viridis')    \n",
    "    im0 = axs[0].imshow(10*np.log10(data2[:,w]),aspect='auto',vmin=noisemin,vmax=noisemax,origin='lower',\n",
    "                        extent=extent,cmap='viridis')\n",
    "    im3 = axs[3].imshow(data1[:,w]/data2[:,w],aspect='auto',vmin=1,vmax=8,origin='lower',\n",
    "                        extent=extent,cmap='viridis')\n",
    "\n",
    "for i in range(0,Nm):\n",
    "    w = np.where((t_set_plt>=raster2_start_mjd[i]) & (t_set_plt<=raster2_stop_mjd[i]))[0]\n",
    "    extent = [raster2_start_mjd[i],raster2_stop_mjd[i],freq[0],freq[-1]]\n",
    "\n",
    "    im2 = axs[2].imshow(10*np.log10(data1[:,w]),aspect='auto',vmin=dBmin,vmax=dBmax,origin='lower',\n",
    "                         extent=extent,cmap='viridis')    \n",
    "    im0 = axs[0].imshow(10*np.log10(data2[:,w]),aspect='auto',vmin=noisemin,vmax=noisemax,origin='lower',\n",
    "                         extent=extent,cmap='viridis')\n",
    "    im3 = axs[3].imshow(data1[:,w]/data2[:,w],aspect='auto',vmin=1,vmax=8,origin='lower',\n",
    "                        extent=extent,cmap='viridis')\n",
    "\n",
    "cbar0= fig1.colorbar(im0,ax=axs[0])\n",
    "cbar0.ax.tick_params(labelsize=fs) \n",
    "cbar0.set_label('Fitted noise source LL Power deflection (dB)', fontsize=fs)\n",
    "\n",
    "cbar2= fig1.colorbar(im2,ax=axs[2])\n",
    "cbar2.ax.tick_params(labelsize=fs) \n",
    "cbar2.set_label('LL Power (dB)', fontsize=fs)\n",
    "\n",
    "cbar3= fig1.colorbar(im3,ax=axs[3])\n",
    "cbar3.ax.tick_params(labelsize=fs) \n",
    "cbar3.set_label('LL Power (noise units)', fontsize=fs)\n",
    "\n",
    "axs[1].plot(t_weath_plt,temp_C,color='purple',linewidth=3)\n",
    "axs[1].set_ylim(-5,35)     \n",
    "axs[1].set_ylabel('Temperature (C)',fontsize=fs)\n",
    "axs[1].axvspan(raster1_start_mjd[0],raster1_stop_mjd[-1],color='C0',alpha=0.1)\n",
    "axs[1].axvspan(raster2_start_mjd[0],raster2_stop_mjd[-1],color='C0',alpha=0.1)\n",
    "ax2 = axs[1].twinx()\n",
    "ax2.scatter(t_noise,10*np.log10(LL_dnoise_arr[:,wf]),s=10,color='red')\n",
    "ax2.scatter(t_noise,10*np.log10(LL_dnoise_fix[:,wf]),s=10,color='k')\n",
    "ax2.set_ylim(noisemin,noisemax)\n",
    "ax2.set_ylabel('Noise source power (LL) at 800 MHz',fontsize=fs)\n",
    "ax2.tick_params(axis='y', labelsize=fs)\n",
    "for i in range(0,len(scan_id)):\n",
    "    elhere = el[abs(t_set_plt-scan_start_mjd[i])<1e-4][0]\n",
    "    if abs(elhere-49.32)<0.5:\n",
    "        clr = 'C1'\n",
    "    elif abs(elhere-20.0)<0.5:\n",
    "        clr = 'C2'\n",
    "    axs[1].axvspan(scan_start_mjd[i],scan_stop_mjd[i],color=clr,alpha=0.1)\n",
    "axs[1].grid() \n",
    "\n",
    "for i in range(0,4):\n",
    "    axs[i].set_xlim(t_set_plt[0],t_set_plt[-1])\n",
    "    axs[i].tick_params(axis='both', labelsize=fs,labelbottom=True)\n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].set_xlabel('Time (UTC)',fontsize=fs)\n",
    "    if i != 1:\n",
    "        axs[i].set_ylim(freq[0],freq[-1])\n",
    "        axs[i].set_ylabel('Frequency (MHz)',fontsize=fs)\n",
    "    \n",
    "#plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/waterfall_phase'+phase+'_day'+day+'.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfe19dfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "del LL\n",
    "del LL_fit\n",
    "\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "662881ff",
   "metadata": {},
   "source": [
    "## Test new file by comparing to old file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79488208",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_test = []\n",
    "az_test = []\n",
    "dec_test = []\n",
    "ra_test = []\n",
    "el_test = []\n",
    "noise_test = []\n",
    "\n",
    "# Use one of the scans to get the list of frequencies:\n",
    "scan0 = f\"{scan_id[0]:04}\"\n",
    "file_test = h5py.File(dir_out_az+'dva_survey_phase'+phase+'_raw_'+scan0+'_gain_corr.h5','r')\n",
    "freq_test = file_test['data']['beam_0']['band_SB0']['frequency'][:]/1e6\n",
    "\n",
    "# Create empty arrays for the power data:\n",
    "RR_test = np.empty([0,len(freq_test)])\n",
    "LL_test = np.empty([0,len(freq_test)])\n",
    "reRL_test = np.empty([0,len(freq_test)])\n",
    "imRL_test = np.empty([0,len(freq_test)])\n",
    "\n",
    "# Raster scan 1:\n",
    "for i in range(0,Ne):\n",
    "    if Ne == 1:\n",
    "        file_test = h5py.File(dir_out_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster1'+'_gain_corr.h5','r')\n",
    "    else:\n",
    "        file_test = h5py.File(dir_out_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster1'+rast_list[i]+'_gain_corr.h5','r')\n",
    "\n",
    "    rast1_file_t_test = [Time(file_test['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][0],\n",
    "                            format='isot',scale='utc').mjd,\n",
    "                        Time(file_test['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][-1],\n",
    "                            format='isot',scale='utc').mjd]\n",
    "    RR_test,LL_test,reRL_test,imRL_test,dec_test,ra_test,el_test,az_test,t_test,noise_test,idx1_test,idx2_test = concatenate_data(file_test,RR_test,LL_test,reRL_test,imRL_test,\n",
    "                                                            dec_test,ra_test,el_test,az_test,t_test,noise_test)\n",
    "\n",
    "\n",
    "# Loop through all the scans in the \"scan_num\" list:\n",
    "for i in scan_id:\n",
    "    file_test = h5py.File(dir_out_az+'dva_survey_phase'+phase+'_raw_'+f\"{i:04}\"+'_gain_corr.h5','r')\n",
    "    print(i,file_test)\n",
    "    RR_test,LL_test,reRL_test,imRL_test,dec_test,ra_test,el_test,az_test,t_test,noise_test,idx1_test,idx2_test = concatenate_data(file_test,RR_test,LL_test,reRL_test,imRL_test,\n",
    "                                                            dec_test,ra_test,el_test,az_test,t_test,noise_test)\n",
    "\n",
    "# Raster scan 2:\n",
    "for i in range(0,Nm):\n",
    "    if Nm == 1:\n",
    "        file_test = h5py.File(dir_out_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster2'+'_gain_corr.h5','r')\n",
    "    else:\n",
    "        file_test = h5py.File(dir_out_rast+'dva_survey_phase'+phase+'_day0'+day+'_raster2'+rast_list[i]+'_gain_corr.h5','r')\n",
    "\n",
    "    rast2_file_t_test = [Time(file_test['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][0],\n",
    "                            format='isot',scale='utc').mjd,\n",
    "                        Time(file_test['data']['beam_0']['band_SB0']['scan_0']['metadata']['utc'][-1],\n",
    "                            format='isot',scale='utc').mjd]\n",
    "    RR_test,LL_test,reRL_test,imRL_test,dec_test,ra_test,el_test,az_test,t_test,noise_test,idx1_test,idx2_test = concatenate_data(file_test,RR_test,LL_test,reRL_test,imRL_test,\n",
    "                                                            dec_test,ra_test,el_test,az_test,t_test,noise_test)\n",
    "\n",
    "\n",
    "t_test_plt = Time(t_test, format='isot',scale='utc').mjd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "489d20ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "fplot = 800\n",
    "df = freq_test[1]-freq_test[0]\n",
    "wf = np.where(abs(freq_test-fplot)<df/2)[0][0]\n",
    "print(freq_test[wf])\n",
    "\n",
    "lims_raw1 = [65,65,65]\n",
    "lims_raw2 = [75,75,75]\n",
    "lims_cor1 = [2,2,2]\n",
    "lims_cor2 = [8,8,8]\n",
    "\n",
    "sz=1\n",
    "fig1, axs = plt.subplots(3,1,figsize=(20,12))\n",
    "\n",
    "for i in range(0,3):\n",
    "    axs[i].scatter(t_test_plt,LL_test[:,wf],s=sz,color='red',alpha=0)\n",
    "    axs[i].set_ylim(lims_raw1[i],lims_raw2[i])\n",
    "        \n",
    "    ax2 = axs[i].twinx()\n",
    "    ax2.scatter(t_test_plt,LL_test[:,wf],s=sz,color='red')\n",
    "    ax2.set_ylim(lims_cor1[i],lims_cor2[i])\n",
    "    axs[i].set_ylabel('Raw power (dB) - not plotted')\n",
    "    \n",
    "    axs[i].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n",
    "    axs[i].fmt_xdata = mdates.DateFormatter('%H:%M:%S')\n",
    "    axs[i].grid()\n",
    "    ax2.set_ylabel('Corrected power (noise src units)')\n",
    "\n",
    "axs[0].set_xlim(raster1_start_mjd[0],raster1_stop_mjd[-1])\n",
    "#axs[0].set_xlim(raster1_start_mjd,raster1_start_mjd+0.005)\n",
    "#axs[0].set_ylim(2e7,4e7)\n",
    "axs[1].set_xlim(scan_start_mjd[0],scan_stop_mjd[-1])\n",
    "#axs[1].set_ylim(1e7,10e7)\n",
    "axs[2].set_xlim(raster2_start_mjd[0],raster2_stop_mjd[-1])\n",
    "#axs[2].set_ylim(1e7,10e7)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('../DVA_PLOTS/noise_gain_corr/correction_example_LL_800MHz_phase'+phase+'_day'+day+'_test.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4476f306",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a410220",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8e279e4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82ba30cd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3429ab0f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8603c511",
   "metadata": {},
   "outputs": [],
   "source": [
    "beam_width = []\n",
    "freq_beam = []\n",
    "\n",
    "with open('../DATA/beam_widths.txt') as fp:\n",
    "    for line in fp:\n",
    "        #print(line.split()[0],line.split()[1])\n",
    "        beam_width.append(float(line.split()[1]))\n",
    "        freq_beam.append(float(line.split()[0]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19218832",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(LL_dnoise_arr.shape)\n",
    "print(len(t_noise))\n",
    "print(LL_fit_noise_pt.shape)\n",
    "\n",
    "ii = 30\n",
    "\n",
    "#print(LL_dnoise_arr[0,:]/LL_fit_noise_pt[0,:])\n",
    "\n",
    "fig1, axs = plt.subplots(1,1,figsize=(20,12))\n",
    "\n",
    "axs.scatter(freq,LL_dnoise_arr[ii,:],s=5)\n",
    "axs.scatter(freq,LL[20000,:],s=5)\n",
    "ax2 = axs.twinx()\n",
    "ax2.scatter(freq_beam,beam_width,s=5,color='C2')\n",
    "axs.set_ylim(0,5e7)\n",
    "ax2.set_ylim(0.5,1.8)\n",
    "axs.set_xlim(350,700)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da6b8ef5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, axs = plt.subplots(1,1,figsize=(20,12))\n",
    "\n",
    "i1 = 30\n",
    "print(t_noise[i1])\n",
    "i2 = 30000\n",
    "print(t_set_plt[i2])\n",
    "\n",
    "axs.scatter(freq,LL_fit[i2,:],s=5)\n",
    "ax2 = axs.twinx()\n",
    "axs.scatter(freq,LL[i2,:],s=5,color='C2')\n",
    "\n",
    "ax2.scatter(freq,LL[i2,:]/LL_fit[i2,:],color='C3',s=5)\n",
    "ax2.scatter(freq_beam,5.0*np.array(beam_width)-3.0,s=5,color='C4')\n",
    "\n",
    "axs.set_ylim(0,4e7)\n",
    "ax2.set_ylim(0,5)\n",
    "axs.set_xlim(350,500)\n",
    "axs.axvline(447.5)\n",
    "axs.axvline(420)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}