Merge pull request #176 from FRBs/alopeke

Code related to the Alopeke analysis by Kilpatrick+2024

papers/Kilpatrick2024_Alopeke/Analysis/py/alopeke_defs.py

@@ -0,0 +1,88 @@
+import numpy as np
+from astropy.time import Time
+from astropy import units
+import json
+
+from frb.frb import FRB
+frb180916 = FRB.by_name('FRB20180916B')
+
+# CHIME
+# Tendulkar (Feb 11, 2021): 
+# This pulse arrival time is topocentric at 400 MHz.
+# In the database I see 2020-10-23 07:48:30.777667 UTC+00:00 
+# The uncertainty is 1 ms. 
+
+chime_time_arrival = {
+    '20201023':Time('2020-10-23T07:48:30.777667'),
+    '20220908':Time('2022-09-08T10:53:26.888596'),
+}
+ata_chime = 1*units.s/1000.  # see comment above
+
+#  Alopeke absolute time accuracy (from email Feb 8, 2021 Nic Scott)
+#  Nic: The major contributor in this uncertainty is thought to be the variable lag
+#  between the computer receipt from the NTP server and the triggering of the
+#  cameras.
+ata_alopeke = 163*units.s/1000.
+
+# final absolute time accuracy (1 sigma)
+ata = np.sqrt(ata_alopeke**2 + ata_chime**2)  # quadrature of absolute uncertainties
+
+# According to the calculation of Kilpatrick, from 400 MHz to optical
+# frequencies, a delay of 9.1s should be applied to the radio.
+# Also account for light-travel time for CHIME to Alopeke (14.8ms)
+chime_time_arrival_optical = {}
+FRB_time = {}
+mjd_low = {}
+mjd_high = {}
+for key in chime_time_arrival.keys():
+    chime_time_arrival_optical[key]=chime_time_arrival[key] - 9.08*units.s - 0.0148*units.s
+    FRB_time[key]=chime_time_arrival_optical[key]
+    mjd_low[key] = (FRB_time[key]-ata).mjd
+    mjd_high[key] = (FRB_time[key]+ata).mjd
+
+# Gains 
+gain_red = 5.61
+gain_blue = 5.54
+EM = 1000
+
+# Center of band
+XSDSS_r = 620. # nm
+XSDSS_i = 765. # nm
+
+# Exposure time
+#dt_alopeke = 11.6 * (1e-3 * units.s)
+dt_alopeke = 10.419 * (1e-3 * units.s)  # on-sky exposure time
+
+# Photometry Star-1
+# From Charlie Kilpatrick #frb180916-speckle channel (5 deb 2021)
+# r=15.7481+/-0.00017 (PS1)
+# i=15.1387+/-0.00237 (PS1)
+r_1 = 15.7481  # Panstarrs-1 magnitude
+i_1 = 15.1387  # Panstarrs-1 magnitude
+
+r_2 = 17.1016  # Panstarrs-1 magnitude
+i_2 = 16.6247  # Panstarrs-1 magnitude
+
+redshift = 0.0337
+
+r_nu = (2.998e18 * units.angstrom/units.second) / (6231 * units.angstrom)
+i_nu = (2.998e18 * units.angstrom/units.second) / (7625 * units.angstrom)
+
+distance = 150.0e6 * units.parsec
+
+# Radio data from CHIME (email from Emmanuel Fonseca, 2023-04-21)
+burst1_file = '../../Data/results_R3/results_fitburst_139459007.json'
+burst2_file = '../../Data/results_R3/results_fitburst_244202260.json'
+
+burst1_data = {}
+burst2_data = {}
+
+with open(burst1_file, 'r') as f:
+    burst1_data = json.load(f)
+with open(burst2_file, 'r') as f:
+    burst2_data = json.load(f)
+
+radio_data = {
+    '20201023':burst1_data,
+    '20220908':burst2_data,
+}

papers/Kilpatrick2024_Alopeke/Analysis/py/alopeke_utils2.py

@@ -0,0 +1,130 @@
+import sys, os
+import numpy as np
+from scipy.stats import norm
+
+from astropy.table import Table
+
+sys.path.append(os.path.abspath("../Analysis/py"))
+import alopeke_defs
+
+def calc_zeropoint(data_dict, camera:str):
+
+    # Instrumental magnitude for star 1 and star 2
+    m_star1 = -2.5 * np.log10(data_dict['C_star1_full']-data_dict['C_bkg_full'])
+    m_star2 = -2.5 * np.log10(data_dict['C_star2_full']-data_dict['C_bkg_full'])
+
+    # Uncertainties on instrumental magnitudes - assume Poisson
+    merr_star1 = 1.086 * np.sqrt(data_dict['C_star1_full'])/data_dict['C_star1_full']
+    merr_star2 = 1.086 * np.sqrt(data_dict['C_star1_full'])/data_dict['C_star1_full']
+
+    # Now get zero points
+    zpt_star1 = 0.
+    zpt_star2 = 0.
+    if camera=='red':
+        zpt_star1 = alopeke_defs.i_1 - m_star1
+        zpt_star2 = alopeke_defs.i_2 - m_star2
+    elif camera=='blue':
+        zpt_star1 = alopeke_defs.r_1 - m_star1
+        zpt_star2 = alopeke_defs.r_2 - m_star2
+    else:
+        raise Exception(f'ERROR: Unrecognized camera {camera}')
+
+    zpt = (merr_star1 * zpt_star1 + merr_star2 * zpt_star2)/(merr_star1 + merr_star2)
+    zpt_err = 1./np.sqrt(2) * np.sqrt(merr_star1**2 + merr_star2**2)
+
+    return(zpt, zpt_err)
+
+def get_star_flux(date:str, data_dict):
+    mask = np.abs(data_dict['MJD_star1_full']-alopeke_defs.FRB_time[date].mjd)*24*3600<alopeke_defs.ata.to('s').value
+
+    data_dict['mean_bkg'] = np.mean(data_dict['C_bkg_full'][mask])
+    data_dict['mean_star1'] = np.mean(data_dict['C_star1_full'][mask])
+    data_dict['mean_star2'] = np.mean(data_dict['C_star2_full'][mask])
+    data_dict['mean_star3'] = np.mean(data_dict['C_star3_full'][mask])
+
+    return(data_dict)
+
+def load_camera(camera:str, date:str, cut=True):
+
+    if camera == 'red':
+        data = Table.read(f'../../Data/master_table_{date}_r.fits')
+    elif camera == 'blue':
+        data = Table.read(f'../../Data/master_table_{date}_b.fits')
+    else:
+        raise IOError("Bad camera")
+
+    # Cut down to Galaxy
+    i_FRB = (data['MJD']>=alopeke_defs.mjd_low[date]) & (
+        data['MJD']<=alopeke_defs.mjd_high[date])
+    i_gal = np.invert(i_FRB)
+
+    # Expunge first reads
+    if cut:
+        gd_data = data['frame'] > 1
+    else:
+        gd_data = np.ones(len(data), dtype='bool')
+
+    # ADUs
+    C_gal = data['flux_2FWHM'][i_gal & gd_data]
+    C_gal_full = data['flux_2FWHM'][gd_data]
+    C_frb = data['flux_2FWHM'][i_FRB & gd_data]
+
+    C_star1 = data['flux_star1_2FWHM'][i_gal & gd_data]
+    C_star1_full = data['flux_star1_2FWHM'][gd_data]
+
+    C_star2 = data['flux_star2_2FWHM'][i_gal & gd_data]
+    C_star2_full = data['flux_star2_2FWHM'][gd_data]
+
+    C_star3 = data['flux_star3_2FWHM'][i_gal & gd_data]
+    C_star3_full = data['flux_star3_2FWHM'][gd_data]
+
+    C_bkg = data['flux_bkg_2FWHM'][i_gal & gd_data]
+    C_bkg_full = data['flux_bkg_2FWHM'][gd_data]
+
+
+    out_dict = {}
+
+    # Time
+    out_dict['MJD_gal'] = data['MJD'][i_gal & gd_data]
+    out_dict['MJD_FRB'] = data['MJD'][i_FRB & gd_data]
+    out_dict['MJD_star1'] = data['MJD'][i_gal & gd_data]
+    out_dict['MJD_star1_full'] = data['MJD'][gd_data]
+    out_dict['MJD_star2'] = data['MJD'][i_gal & gd_data]
+    out_dict['MJD_star2_full'] = data['MJD'][gd_data]
+    out_dict['MJD_star3'] = data['MJD'][i_gal & gd_data]
+    out_dict['MJD_star3_full'] = data['MJD'][gd_data]
+    out_dict['dt'] = (data['MJD'][2]-data['MJD'][1])*24*3600 # seconds
+
+    # Counts
+    gain = alopeke_defs.gain_red if camera == 'red' else alopeke_defs.gain_blue
+    out_dict['C_FRB'] = C_frb * gain / alopeke_defs.EM
+    out_dict['C_gal'] = C_gal* gain / alopeke_defs.EM
+    out_dict['C_gal_full'] = C_gal_full * gain / alopeke_defs.EM
+
+    out_dict['C_star1'] = C_star1 * gain / alopeke_defs.EM
+    out_dict['C_star1_full'] = C_star1_full * gain / alopeke_defs.EM
+
+    out_dict['C_star2'] = C_star2 * gain / alopeke_defs.EM
+    out_dict['C_star2_full'] = C_star2_full * gain / alopeke_defs.EM
+
+    out_dict['C_star3'] = C_star3 * gain / alopeke_defs.EM
+    out_dict['C_star3_full'] = C_star3_full * gain / alopeke_defs.EM
+
+    out_dict['C_bkg'] = C_bkg * gain / alopeke_defs.EM
+    out_dict['C_bkg_full'] = C_bkg_full * gain / alopeke_defs.EM
+
+
+    # Gauss
+    out_dict['Gauss_gal'] = norm.fit(out_dict['C_gal'])
+    out_dict['Gauss_star1'] = norm.fit(out_dict['C_star1'])
+    out_dict['Gauss_star2'] = norm.fit(out_dict['C_star2'])
+    out_dict['Gauss_star3'] = norm.fit(out_dict['C_star3'])
+
+    out_dict = get_star_flux(date, out_dict)
+
+    zpt, zpt_err = calc_zeropoint(out_dict, camera)
+    out_dict['zeropoint'] = zpt
+    out_dict['zeropoint_error'] = zpt_err
+
+    # Return
+    return out_dict

papers/Kilpatrick2024_Alopeke/Analysis/py/analysis_data_20201023.txt

@@ -0,0 +1,46 @@
+This camera = red
+Time of start (UTC) 2020-10-23T07:41:40.265
+We started 6.841878353152424 minutes before
+And ended 12.09422416635789 minutes after
+N values (+/-0.16300306745579976 s): 28
+Max flux: 72.33394791687047
+Mean flux: 47.82070429999975
+Standard deviation: 9.000776067207367
+Zeropoint: 22.136190941626683 AB mag
+Zeropoint error: 0.035881878355318225 mag
+The upper limit is 51.61662613981765
+Extinction (A_filter): 1.5465383884105393
+AB mag (no dust correction): 18.18696094759817
+AB mag (MW dust correction): 16.64042255918763
+fluence (no dust correction): 2.0092855408094445 s uJy
+fluence (MW dust correction): 8.349433582310564 s uJy
+fluence ratio (opt/radio): 0.003211320608580986
+Equivalent isotropic energy: 8.837761298901163e+40 erg
+Intercept is: 908.6946900880495 electrons/s
+Slope is: -3.99180150411254 electrons/s
+The upper limit is 51.66037543041208
+Fraction exceeding: 0.0678
+Total observing = 1136.1661511706188s
+This camera = blue
+Time of start (UTC) 2020-10-23T07:41:40.265
+We started 6.84188335086219 minutes before
+And ended 12.095297552878037 minutes after
+N values (+/-0.16300306745579976 s): 28
+Max flux: 56.42290357139614
+Mean flux: 43.515621034091886
+Standard deviation: 7.255609825974984
+Zeropoint: 22.151109428099474 AB mag
+Zeropoint error: 0.04506495957178212 mag
+The upper limit is 38.949428571428534
+Extinction (A_filter): 2.205121094427004
+AB mag (no dust correction): 18.59129404177557
+AB mag (MW dust correction): 16.386172947348566
+fluence (no dust correction): 1.3845492902092602 s uJy
+fluence (MW dust correction): 10.552536404938659 s uJy
+fluence ratio (opt/radio): 0.00405866784805333
+Equivalent isotropic energy: 1.366860549177359e+41 erg
+Intercept is: 576.6310231475677 electrons/s
+Slope is: 4.910245350656469 electrons/s
+The upper limit is 38.94579361689137
+Fraction exceeding: 0.6364
+Total observing = 1136.2308542244136s

papers/Kilpatrick2024_Alopeke/Analysis/py/analysis_data_20220908.txt

@@ -0,0 +1,46 @@
+This camera = red
+Time of start (UTC) 2022-09-08T10:39:04.629
+We started 14.370993955526501 minutes before
+And ended 7.545006988802925 minutes after
+N values (+/-0.16300306745579976 s): 28
+Max flux: 99.15800454694453
+Mean flux: 88.44920933888261
+Standard deviation: 8.789814983629809
+Zeropoint: 22.118906012299913 AB mag
+Zeropoint error: 0.03531093383954578 mag
+The upper limit is 41.997860785932886
+Extinction (A_filter): 1.5465383884105393
+AB mag (no dust correction): 18.279726116280635
+AB mag (MW dust correction): 16.733187727870096
+fluence (no dust correction): 1.8447418125932844 s uJy
+fluence (MW dust correction): 7.665684606755237 s uJy
+fluence ratio (opt/radio): 0.002190195601930068
+Equivalent isotropic energy: 8.114022356042966e+40 erg
+Intercept is: 825.3673651910062 electrons/s
+Slope is: 7.763991414138952 electrons/s
+The upper limit is 41.988571154215386
+Fraction exceeding: 0.9929
+Total observing = 1314.9600566597655s
+This camera = blue
+Time of start (UTC) 2022-09-08T10:39:04.629
+We started 14.370995422359556 minutes before
+And ended 7.546799898846075 minutes after
+N values (+/-0.16300306745579976 s): 28
+Max flux: 134.17818092168486
+Mean flux: 112.8543244232961
+Standard deviation: 16.212168141402877
+Zeropoint: 21.993747416781208 AB mag
+Zeropoint error: 0.042791267583366197 mag
+The upper limit is 76.1798582600194
+Extinction (A_filter): 2.205121094427004
+AB mag (no dust correction): 17.659275130073198
+AB mag (MW dust correction): 15.454154035646194
+fluence (no dust correction): 3.266750648610154 s uJy
+fluence (MW dust correction): 24.89799777377808 s uJy
+fluence ratio (opt/radio): 0.007113713649650881
+Equivalent isotropic energy: 3.2250152574271826e+41 erg
+Intercept is: 555.8676614392847 electrons/s
+Slope is: -1.9068632522228652 electrons/s
+The upper limit is 76.23242942144554
+Fraction exceeding: 0.9591
+Total observing = 1315.0677192723379s

papers/Kilpatrick2024_Alopeke/Data/README

@@ -0,0 +1,3 @@
+Grab the files from the Drive.
+
+Don't add to the Repo