photometry_final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:14:13.665007Z",
     "start_time": "2023-08-27T20:14:11.492986Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import csv\n",
    "import glob\n",
    "import sep\n",
    "import astropy.units as u\n",
    "import astropy.coordinates as coord\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import astroalign as aa\n",
    "from scipy import odr\n",
    "import matplotlib.colors as mcolors\n",
    "from astropy.io import fits\n",
    "from astroquery.ipac.irsa import Irsa\n",
    "from astroquery.sdss import SDSS\n",
    "from astropy.wcs import WCS\n",
    "from astropy.wcs.utils import pixel_to_skycoord\n",
    "from astropy.coordinates import SkyCoord\n",
    "from matplotlib.colors import LogNorm, TABLEAU_COLORS, CSS4_COLORS\n",
    "from photutils.aperture import aperture_photometry, CircularAperture, CircularAnnulus, ApertureStats\n",
    "from regions import CirclePixelRegion, PixCoord\n",
    "from photutils import centroids\n",
    "\n",
    "\n",
    "#imports\n",
    "#FOR ASTROQUERRY: If you encounter an error related to \"fields\" in \"SDSS.querry_crossid\", try installing the newest development version of astroquerry rather than the default. There is a bug in the older official release version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:14:13.680191Z",
     "start_time": "2023-08-27T20:14:13.676607Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Night: #Each observing night is initialized as a Night object\n",
    "    def __init__(self, file_path, night_number):\n",
    "        #each attribute is declared here first, even if value is assigned later\n",
    "        self.path = file_path #file path of obs night\n",
    "        self.obs_night = night_number #each night is assigned a number, the first night is night 0\n",
    "        self.image_data = None #the image data for the night\n",
    "        self.headers = None #the night's headers\n",
    "        self.wcs = None #world coordinate system transformation matrix\n",
    "        self.readout_noise = None #detector readout noise\n",
    "        self.aligned = None #aligned image data\n",
    "        self.template = None #template frame\n",
    "        self.references = None #reference sources\n",
    "        self.obs_filter = None #observation filter\n",
    "        self.no_obs_images = None #number of images in night\n",
    "\n",
    "\n",
    "    def initialize_frames(self):\n",
    "        science_dir = sorted(glob.glob(self.path + '/*')) #pulls data files using glob\n",
    "        hdus = [fits.open(image) for image in science_dir] #opens each fits image\n",
    "        self.image_data = [image[1].data for image in hdus] #pulls image data from file\n",
    "        self.headers = [image[1].header for image in hdus] #pulls header data from file.\n",
    "        self.wcs = WCS(self.headers[0]) #gets wcs transformation, each night should only have one unique transformation\n",
    "        self.readout_noise = self.headers[0]['RDNOISE'] #pulls readout noise from header. Readout noise is detector based and should be the same across nights taken with same equipment but do this just in case.\n",
    "        self.mjd_times = []\n",
    "        self.date_times = []\n",
    "        self.date = []\n",
    "        self.start_times = []\n",
    "        for header in self.headers:\n",
    "            self.mjd_times.append(header['MJD-OBS']) #Modified Julian Date (date+start time)\n",
    "            self.date_times.append(header['DATE-OBS']) #Date/Start time (YYYY-MM-DD:HH:MM:SS in UTC)\n",
    "            self.date.append(header['DAY-OBS']) #Date (YYYYMMDD)\n",
    "            self.start_times.append(header['UTSTART']) #UTC start time HH:MM:SS\n",
    "        try:\n",
    "            self.aligned_images = self.align()\n",
    "        except:\n",
    "            print(\"Unable to align using affine transform. WCS matrix may be unavailble or points out of image bounds\")\n",
    "            try:\n",
    "                self.aligned_images = [aa.register(image, self.image_data[0])[0] for image in self.image_data[0:]]\n",
    "            except:\n",
    "                print(\"Could not align images\")\n",
    "                self.aligned_images = self.image_data\n",
    "                self.is_aligned = False\n",
    "\n",
    "        self.template = np.median(self.aligned_images, axis = 0) #night template\n",
    "        background = sep.Background(self.template) #sep background subtraction for source extraction\n",
    "        self.references = sep.extract(self.template - background.back(),  background.globalrms*3, minarea =25, segmentation_map=False) #finds sources based on given parameters\n",
    "\n",
    "        self.obs_filter = self.headers[0]['filter'][0] #observation filter.\n",
    "        self.no_obs_images = len(self.aligned_images) #number of images in night.\n",
    "\n",
    "\n",
    "    def get_info(self): #function to grab night info, useful for debugging.\n",
    "        print(f\"path: {self.path}, night {self.obs_night}, n_frames: {len(self.image_data)}, n_aligned: {len(self.aligned_images)}, wcs: {self.wcs}, n_ref: {len(self.references)}, filter: {self.obs_filter}\")\n",
    "\n",
    "\n",
    "    def align(self):\n",
    "        reference_pixels = np.array([[100, 200], [1500, 1500], [1800, 3000]]).astype(np.float32) #these are arbitrary points spanning the image, should be done dynamically at some point\n",
    "        reference_skycoords = []\n",
    "        for point in reference_pixels:\n",
    "            skyval = pixel_to_skycoord(point[1], point[0], wcs = self.wcs[0])\n",
    "            reference_skycoords.append(skyval)\n",
    "        \n",
    "        target_points = []\n",
    "        for i in range(len(self.image_data)):\n",
    "            target_pixels = []\n",
    "            for point in reference_skycoords:\n",
    "                pixval = SkyCoord.to_pixel(point, wcs= self.wcs[i])\n",
    "                target_pixels.append(pixval)\n",
    "            target_points.append(np.array(target_pixels).astype(np.float32))\n",
    "\n",
    "\n",
    "        affine_aligned = []\n",
    "        for i, image in enumerate(self.image_data):\n",
    "            M = cv2.getAffineTransform(reference_pixels, target_points[i])\n",
    "            rows, cols = image.shape\n",
    "            dst = cv2.warpAffine(image, M, (rows, cols))\n",
    "            affine_aligned.append(dst.T)\n",
    "\n",
    "        return affine_aligned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:14:13.706323Z",
     "start_time": "2023-08-27T20:14:13.680870Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Source: #initialize source object\n",
    "    def __init__(self, source, count, WCS):\n",
    "        self.position = pixel_to_skycoord(source['x'], source['y'], wcs= WCS).transform_to('icrs') #since pixel locations are inconsistent, store position as RA/DEC\n",
    "        self.radius = (source['xmax'] - source['xmin'])/2 #source radius (size) provided by SEP\n",
    "        self.source_id = count #identifying number\n",
    "        self.is_reference = None #if star is reference\n",
    "        self.ref_mag = None #SDSS magnitude, if available\n",
    "        self.ref_mag_err = None #reference mag error\n",
    "        self.inst_mags = [] #instrumnetal (our) magnitudes\n",
    "        self.inst_mag_errs = [] #instrumental mag errors\n",
    "        self.calibrated_mags = [] #calibrated magnitudes\n",
    "        self.flagged = False #bad source flag. Will be flipped true if a source is not present in all observing nights or has negative aperture sum\n",
    "        self.weights = []\n",
    "        self.errors = []\n",
    "        self.chi2 = None\n",
    "\n",
    "\n",
    "    def query_source(self): #querry a source through the sdss database\n",
    "        #we want the search to return ra, dec, mags, and mag error. region = False, returns first result of search.\n",
    "        search = SDSS.query_crossid(self.position, fields = ['ra', 'dec', f'psfMag_{Nights[0].obs_filter}', f'psfMagErr_{Nights[0].obs_filter}'], radius = 15 * u.arcsec, region = False)\n",
    "        if search:\n",
    "            if search['type'] == \"STAR\": #want to make sure reference objects are Stars\n",
    "                self.is_reference = True #lets us know star is a reference star\n",
    "                self.ref_mag = search[f'psfMag_{Nights[0].obs_filter}'] #fill in mag and error fields:\n",
    "                self.ref_mag_err = search[f'psfMagErr_{Nights[0].obs_filter}']\n",
    "\n",
    "    def boundary_check(self, night): #checks if a star is within frame for a given night.\n",
    "        source_xy = SkyCoord.to_pixel(self.position, wcs= night.wcs)\n",
    "        if (night.headers[0]['NAXIS1'] - source_xy[0]) < 0 or source_xy[0] < 0 or (night.headers[0]['NAXIS2'] - source_xy[1]) < 0 or source_xy[1] < 0:\n",
    "            self.flagged = True #if star is out of bounds, flags star as bad\n",
    "\n",
    "\n",
    "    def aperture_photometry(self, img, nght): #does aperture photometry\n",
    "\n",
    "        coords = SkyCoord.to_pixel(self.position, wcs = nght.wcs) #gets pixel values of source from RA DEC\n",
    "        pcoords = PixCoord(coords[0], coords[1]) #another coord object needed for Regions\n",
    "\n",
    "        radius_i = self.radius #inner aperture radius\n",
    "        radius_o_0 = radius_i + 5 #inner annulus radius\n",
    "        radius_o_1 = radius_o_0 + 5 #outer annulus radius\n",
    "\n",
    "        source_circle = CirclePixelRegion(pcoords, radius_i) #makes region of source shape\n",
    "        source_circle_mask = source_circle.to_mask()\n",
    "        source_aperture = source_circle_mask.cutout(img) #gets data of source\n",
    "        source_sum_unsub = np.sum(source_aperture)\n",
    "\n",
    "        background_annulus = CircularAnnulus(coords, radius_o_0, radius_o_1) #makes annulus for background subtraction\n",
    "        #background_mean = ApertureStats(img, background_annulus).mean #takes mean of background annulus\n",
    "        background_sum = aperture_photometry(img, background_annulus)['aperture_sum'][0]\n",
    "\n",
    "        #source_flux_pix = source_aperture-((source_circle.area/background_annulus.area)*background_sum*source_circle_mask) #pixel wise background subtraction\n",
    "        source_flux_total = np.sum(source_aperture) - (source_circle.area/background_annulus.area)*background_sum  #total flux\n",
    "\n",
    "\n",
    "\n",
    "        readout_sum_source = source_circle.area*(nght.readout_noise**2)\n",
    "        readout_sum_annulus = background_annulus.area*(nght.readout_noise**2)\n",
    "\n",
    "        delta_n = (readout_sum_source + source_flux_total + ((source_circle.area/background_annulus.area)**2)*(readout_sum_annulus + background_sum))**(1/2)\n",
    "\n",
    "        if source_flux_total < 0:\n",
    "            print(self.source_id)\n",
    "            self.flagged = True #flags source if aperture sum turns out to be negative\n",
    "\n",
    "        else:\n",
    "            instrumental_mag = -2.5*np.log10(source_flux_total) #magnitude\n",
    "            instrumental_mag_error = 2.5*np.log10(np.e)*abs(delta_n/source_flux_total) #magntiude error\n",
    "            self.inst_mags.append(instrumental_mag)\n",
    "            self.inst_mag_errs.append(instrumental_mag_error)\n",
    "\n",
    "\n",
    "    def add_calibrated_mag(self, mag):\n",
    "        self.calibrated_mags.append(mag) #adds calibrated mag. For some reason, math comes out unexpectedly if calibration takes place in class function.\n",
    "\n",
    "    def add_chi(self, chi):\n",
    "        self.chi2 = chi\n",
    "\n",
    "    def add_error(self, err):\n",
    "        self.errors.append(err)\n",
    "        self.weights.append(1/(err**2))\n",
    "\n",
    "    def get_info(self): #prints out source info.\n",
    "        print(f\"ra_dec: {self.position}, Night_0_xy: {SkyCoord.to_pixel(self.position, wcs = Nights[0].wcs)} rad: {self.radius}, ref_status: {self.is_reference}, ref_mag: {self.ref_mag}, inst_mag_avg:{np.mean(self.inst_mags)}, cal_mag_avg: {np.mean(self.calibrated_mags)}, flagged: {self.flagged}, ID: {self.source_id}, Chi2: {self.chi2}\")\n",
    "\n",
    "    def __iter__(self): #for writing out csv files\n",
    "        return iter([self.position, self.is_reference, self.ref_mag, self.chi2, self.flagged, self.source_id, self.calibrated_mags])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:14:13.710152Z",
     "start_time": "2023-08-27T20:14:13.701561Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def lin_model(p, x): #define a standard linear model for ODR fitting. Part of calibration.\n",
    "    return p[0] * x + p[1]\n",
    "\n",
    "def ODR(x_data, y_data):\n",
    "    x_bar = np.mean(x_data)\n",
    "    y_bar = np.mean(y_data)\n",
    "\n",
    "    s_xx = 1/len(x_data) * np.sum((x_data - x_bar)**2)\n",
    "    s_yy = 1/len(y_data) * np.sum((y_data - y_bar)**2)\n",
    "    s_xy = 1/len(x_data) * np.sum((x_data - x_bar) * (y_data - y_bar))\n",
    "\n",
    "    b_0 = (s_yy - s_xx + np.sqrt((s_yy - s_xx)**2 + 4*s_xy**2))/(2 * s_xy)\n",
    "    b_1 = y_bar - b_0 * x_bar\n",
    "\n",
    "    return [b_0, b_1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:19:36.998838Z",
     "start_time": "2023-08-27T20:14:13.712323Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# primary_dir = sorted(glob.glob(\"/Users/lucaangeleri/Documents/LCO/SN2023ixf_all/SN2023ixf_r/*\")) #main file directory; directory containing folders for each night.\n",
    "primary_dir = sorted(glob.glob(r\"C:\\Users\\Sam Whitebook\\Documents\\Visual Studio 2010\\Projects\\Lubin Lab\\Data\\sn2023ixf\\sn2023ixf_r\\*\")) #main file directory; directory containing folders for each night.\n",
    "Nights = [Night(directory, dir_number) for dir_number, directory in enumerate(primary_dir)] #initializes each night as a Night object\n",
    "for night in Nights:\n",
    "    night.initialize_frames() #see initialize_frames() in Night class definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:25:34.592438Z",
     "start_time": "2023-08-27T20:25:34.277612Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Sources = [Source(source, count, Nights[0].wcs) for count, source in enumerate(Nights[0].references)] #initializes sources based off first night's list. This ensures proper source tracking\n",
    "for source in Sources:\n",
    "    source.query_source() #see query_source() in source class definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:26:23.635481Z",
     "start_time": "2023-08-27T20:25:36.778963Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "night_array = [] #this is to help organize plotting later.\n",
    "mag_thresh = 15 #magnitude threshold for calibrating sources.\n",
    "for night in Nights: #for each night, iterates through every source for each image.\n",
    "    for image in night.aligned_images:\n",
    "        for source in Sources:\n",
    "            source.boundary_check(night) #see boundary_check() in source class definition\n",
    "            if source.flagged == False:\n",
    "                source.aperture_photometry(image, night)  #see aperture_photometry() in source class definition\n",
    "        night_array.append(night.obs_night)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:27:05.610202Z",
     "start_time": "2023-08-27T20:27:05.059196Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "slopes = []\n",
    "zeros = []\n",
    "slope_errs = []\n",
    "zero_errs = []\n",
    "\n",
    "counter = 0\n",
    "\n",
    "for night in Nights:\n",
    "    for image in night.aligned_images:\n",
    "        instrumental_magnitudes = [s.inst_mags[counter] for s in Sources if s.is_reference == True and s.ref_mag < mag_thresh and s.flagged != True]\n",
    "        reference_magnitudes = [s.ref_mag[0] for s in Sources if s.is_reference == True and s.ref_mag < mag_thresh and s.flagged != True]\n",
    "        jk_params = np.zeros((len(instrumental_magnitudes), 2))\n",
    "        for i in range(len(instrumental_magnitudes)):\n",
    "            x_sample = np.append(instrumental_magnitudes[:i], instrumental_magnitudes[i+1:])\n",
    "            y_sample = np.append(reference_magnitudes[:i], reference_magnitudes[i+1:])\n",
    "            jk_params[i] = ODR(x_sample, y_sample)\n",
    "\n",
    "        mean_params = np.mean(jk_params, axis = 0)\n",
    "        sig_params = np.std(jk_params, axis = 0)\n",
    "\n",
    "        slopes.append(mean_params[0])\n",
    "        zeros.append(mean_params[1])\n",
    "        slope_errs.append(sig_params[0])\n",
    "        zero_errs.append(sig_params[1])\n",
    "\n",
    "        counter += 1\n",
    "\n",
    "\n",
    "\n",
    "print(len(slopes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:27:08.841377Z",
     "start_time": "2023-08-27T20:27:08.837066Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(f\"Slopes: {slopes[5]}\")\n",
    "print(f\"Zero: {zeros[5]}\")\n",
    "print(f\"Slope ERRs: {slope_errs[5]}\")\n",
    "print(f\"Zero errs: {zero_errs[5]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:27:11.732218Z",
     "start_time": "2023-08-27T20:27:11.584220Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for source in Sources:\n",
    "    if source.flagged != True:\n",
    "        for i in range(0, len(slopes)):\n",
    "            mag = (source.inst_mags[i]*slopes[i] + zeros[i])\n",
    "            # final_err = np.sqrt((slopes[i]* source.inst_mags[i])**2 * ((slope_errs[i]/slopes[i])**2 + (source.inst_mag_errs[i]/source.inst_mags[i])**2) + zero_errs[i]**2)\n",
    "            final_err = source.inst_mag_errs[i]\n",
    "            source.add_calibrated_mag(mag)\n",
    "            print(final_err)\n",
    "            source.add_error(final_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:27:34.886346Z",
     "start_time": "2023-08-27T20:27:34.658592Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "color_arr = np.array(sorted(CSS4_COLORS, key=lambda c: tuple(mcolors.rgb_to_hsv(mcolors.to_rgb(c)))))\n",
    "l = np.arange(0, len(slopes))\n",
    "\n",
    "median_curves = []\n",
    "for night in Nights:\n",
    "    median_mags = []\n",
    "    index_l = len(night.image_data)*night.obs_night\n",
    "    index_h = len(night.image_data)*(1+night.obs_night)\n",
    "    for source in Sources:\n",
    "        if source.flagged != True:\n",
    "            median_mags.append(source.calibrated_mags[index_l:index_h])\n",
    "    median_curve = np.median(median_mags, axis = 0)/np.median(median_mags)\n",
    "    median_curves.append(median_curve)\n",
    "med_curve = np.concatenate(median_curves)\n",
    "plt.scatter(l, med_curve, c = color_arr[93::5][np.array(night_array)])\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()\n",
    "\n",
    "print(len(med_curve))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:27:57.289054Z",
     "start_time": "2023-08-27T20:27:35.645664Z"
    }
   },
   "outputs": [],
   "source": [
    "for source in Sources:\n",
    "    if source.flagged != True:\n",
    "        if len(Nights) > 1:\n",
    "            print(source.get_info())\n",
    "            night_array = np.array(night_array)\n",
    "            source_mags = np.array(source.calibrated_mags/med_curve)\n",
    "            source_errs = np.array(source.errors)\n",
    "            median_mags = []\n",
    "            median_errs = []\n",
    "            mjd_times = []\n",
    "            utc_times = []\n",
    "            for night in range(np.max(night_array) + 1): # +1 avoids fencepost error in the night because range is exclusive.\n",
    "                #Magnitudes\n",
    "                index_array = np.where(night_array == night)[0]\n",
    "                night_mags = source_mags[index_array]\n",
    "                night_errs = source_errs[index_array]\n",
    "                median_mags.append(np.median(night_mags))\n",
    "                median_errs.append(np.median(night_errs))\n",
    "                #Times\n",
    "                mjd_array = np.array(Nights[night].mjd_times)\n",
    "                utc_array = np.array(Nights[night].date)\n",
    "                mjd_times.append(np.median(mjd_array))\n",
    "                utc_times.append(utc_array[int(len(utc_array)/2)])\n",
    "            print(median_errs)\n",
    "            avg_mag = np.average(source.calibrated_mags/med_curve, weights= source.weights)\n",
    "            Chis = []\n",
    "            for i, m in enumerate(source.calibrated_mags):\n",
    "                chi_i = (((m/med_curve[i]) - avg_mag)**2)/(source.errors[i]**2)\n",
    "                Chis.append(chi_i)\n",
    "            dof = len(source.calibrated_mags) - 1\n",
    "            chi_dof = np.sum(Chis)/dof\n",
    "            source.add_chi(chi_dof)\n",
    "            dof_string =  \"%.2f\" % chi_dof\n",
    "            plt.figure(figsize=(10,7.5))\n",
    "            plt.title(f\"Source: {source.source_id}, Location: {source.position}, Chi2/Dof: {dof_string}\")\n",
    "            plt.errorbar(mjd_times, median_mags, yerr=median_errs, elinewidth=2, capsize=5, linestyle=\"\", marker=\"o\", color=\"black\")\n",
    "            plt.xticks(mjd_times, utc_times, rotation=45, ha='right')\n",
    "            plt.xlabel(\"Observation Date (YYYYMMDD)\")\n",
    "            plt.ylabel(f\"Magnitude, {Nights[0].obs_filter}-band\")\n",
    "            plt.gca().invert_yaxis()\n",
    "            plt.tight_layout()\n",
    "            plt.plot(mjd_times, np.ones(len(mjd_times))*avg_mag, linestyle = '--', color = 'black', label = \"TRIPP Average Mag: {}\".format(\"%.3f\" % avg_mag))\n",
    "            if source.is_reference:\n",
    "                plt.plot(mjd_times, np.ones(len(mjd_times))*source.ref_mag, linestyle = 'dashdot', color = f\"{Nights[0].obs_filter}\", label = \"SDSS Mag: {}\".format(\"%.3f\" % source.ref_mag))\n",
    "            plt.legend()\n",
    "            #plt.savefig(f\"/Users/lucaangeleri/Desktop/test/source{source.source_id}.png\", dpi = 1000)\n",
    "            plt.show()\n",
    "\n",
    "        else: #single night\n",
    "            print(source.get_info())\n",
    "            r = np.arange(0, len(source.calibrated_mags))\n",
    "            avg_mag = np.average(source.calibrated_mags/med_curve, weights=source.weights)\n",
    "            #avg_mag = np.mean(source.calibrated_mags/med_curve)\n",
    "            Chis = []\n",
    "            for i, m in enumerate(source.calibrated_mags):\n",
    "                chi_i = (((m / med_curve[i]) - avg_mag) ** 2) / (source.errors[i] ** 2)\n",
    "                Chis.append(chi_i)\n",
    "            dof = len(source.calibrated_mags) - 1\n",
    "            chi_dof = np.sum(Chis) / dof\n",
    "            dof_string = \"%.2f\" % chi_dof\n",
    "            plt.figure(figsize=(12, 9))\n",
    "            plt.errorbar(r, source.calibrated_mags/med_curve, yerr=source.errors, elinewidth=1, capsize=2, markersize = 3, linestyle = 'none', marker = 'o', c = 'black')\n",
    "            plt.plot(r, np.ones(len(r))*avg_mag, linestyle = '--', color = 'black', label = f\"TRIPP Avg Mag:{avg_mag}\")\n",
    "            if source.is_reference:\n",
    "                plt.plot(r, np.ones(len(r))*source.ref_mag, linestyle = 'dashdot', color = f\"{Nights[0].obs_filter}\", label = \"SDSS Mag: {}\".format(\"%.3f\" % source.ref_mag))\n",
    "            x_locs= []\n",
    "            times = []\n",
    "            for i in range(len(Nights[0].image_data))[::10]:\n",
    "                x_locs.append(i)\n",
    "                times.append(Nights[0].start_times[i])\n",
    "            plt.xlabel(\"Observation Start Time, UTC (HH:MM:SS)\")\n",
    "            plt.ylabel(f\"Magnitude {Nights[0].obs_filter}-band\")\n",
    "            plt.title(f\"Source: {source.source_id}, Location: {source.position}, Chi2/Dof: {chi_dof}\")\n",
    "            plt.gca().invert_yaxis()\n",
    "            plt.xticks(x_locs, times, rotation = 45)\n",
    "            plt.legend()\n",
    "            #plt.savefig(\"/Users/lucaangeleri/Desktop/test/sec32source{}.png\".format(source.source_id), format = 'png', dpi = 500,  bbox_inches=\"tight\")\n",
    "            plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:20:46.104829Z",
     "start_time": "2023-08-27T20:20:46.094713Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "transient_candidates = []\n",
    "for source in Sources:\n",
    "    if source.flagged != True:\n",
    "        if source.chi2 > 25:\n",
    "            transient_candidates.append(source.get_info())\n",
    "\n",
    "print(len(transient_candidates))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:20:46.113198Z",
     "start_time": "2023-08-27T20:20:46.107876Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "header_row = [\"Source_Position\", \"Reference_Source\", \"Reference_Magnitude\", \"Chi2/DoF\", \"Flagged\", \"Source_ID\", \"TRIPP_Calibrated_Mags\" ]\n",
    "with open('/Users/lucaangeleri/Desktop/output.csv', 'w') as stream:\n",
    "    writer = csv.writer(stream)\n",
    "    writer.writerow(header_row)\n",
    "    writer.writerows([source for source in Sources if source.flagged != True and source.chi2>25])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:20:46.120161Z",
     "start_time": "2023-08-27T20:20:46.115639Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for night in Nights:\n",
    "    print(night.headers[0]['SITE'], night.headers[0]['TELID'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-27T20:20:46.222297Z",
     "start_time": "2023-08-27T20:20:46.122362Z"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "chis = []\n",
    "for source in Sources:\n",
    "    if not source.flagged:\n",
    "        chis.append(source.chi2)\n",
    "print(np.mean(chis), np.std(chis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "inst_mags = [source.inst_mags[counter] for source in Sources if source.is_reference == True and source.ref_mag < mag_thresh and source.flagged == False]\n",
    "        #inst_errs = [source.inst_mag_errs[counter] for source in Sources if source.is_reference == True and source.ref_mag < mag_thresh and source.flagged == False]\n",
    "        ref_mags = [source.ref_mag[0] for source in Sources if source.is_reference == True and source.ref_mag < mag_thresh and source.flagged == False]\n",
    "        #ref_mag_errs = [source.ref_mag_err[0] for source in Sources if source.is_reference == True and source.ref_mag < mag_thresh and source.flagged == False]\n",
    "        p0, cov0 = np.polyfit(inst_mags, ref_mags, deg = 1, cov = True)\n",
    "        diffs = [s.ref_mag[0] - (s.inst_mags[counter]*p0[0] + p0[1]) for s in Sources if s.is_reference == True and s.ref_mag < mag_thresh and source.flagged == False]\n",
    "        stdv = np.std(diffs)\n",
    "        inst_mags_final = []\n",
    "        sky_mags_final = []\n",
    "        outlier_inst =[]\n",
    "        outlier_sky =[]\n",
    "        for diff in diffs: #rudementary sigma clipping to remove outliers from calibration model.\n",
    "            if diff < stdv:\n",
    "                i = diffs.index(diff)\n",
    "                inst_mags_final.append(inst_mags[i])\n",
    "                sky_mags_final.append(ref_mags[i])\n",
    "            else:\n",
    "                i = diffs.index(diff)\n",
    "                outlier_inst.append(inst_mags[i])\n",
    "                outlier_sky.append(ref_mags[i])\n",
    "        if len(sky_mags_final) > 0:\n",
    "            p1, cov1 = np.polyfit(inst_mags_final, sky_mags_final, deg = 1, cov = True)\n",
    "            if abs(p1[0]-1) < abs(p0[1]-1):\n",
    "                print(\"using second\")\n",
    "                slopes.append(p1[0])\n",
    "                zeros.append(p1[1])\n",
    "                slope_errs.append(np.sqrt(np.diag(cov1))[0])\n",
    "                zero_errs.append(np.sqrt(np.diag(cov1))[1])\n",
    "        else:\n",
    "            slopes.append(p0[0])\n",
    "            zeros.append(p0[1])\n",
    "            slope_errs.append(np.sqrt(np.diag(cov0))[0])\n",
    "            zero_errs.append(np.sqrt(np.diag(cov0))[1])\n",
    "        counter += 1"
   ]
  }
 ],
 "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}