Goal: We need to define some criteria for a "detection," which is going to be related to how bright the source is relative to the detection limit. For now, let's be conservative and say a source is detected only if it is brighter than the limits you found by 0.5 mag. Let's also start by observing in only a single filter. So, to simulate observations I'd like for you to do the following - pick a rondom number between 0 and 24 (this is the time in hours after which LIGO detected the event). Then observe the light curve (for now just use the z = 0.015 model, though we will vary this eventually) in one of the filters, do this again 24 hr later, and 6 d after that. If we "detect" the source in each of the first 2 observations, but we do not detect the source in the last observation, consider the source "discovered". What we are now trying to make is a histogram as a function of 0-24 hr for each filter for how likely we are to detect the kilonova.
In [30]:
#Necessary files needed for plotting the curves and manipulating the data.
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import astropy as astro
from astropy.table import Table,Column
from astropy.io import ascii
from astropy.cosmology import WMAP9 as cosmo
from scipy.interpolate import interp1d
plt.rcParams["figure.figsize"] = (10,8)
#We will need to grab the data from where it is saved on the computer.
#NS-NS
APR1215 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/APR4-1215.txt')
APR1314 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/APR4-1314.txt')
H41215 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/H4-1215.txt')
H41314 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/H4-1314.txt')
Sly135 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/Sly-135.txt')
#NS-BH
APR4Q3a75 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/APR4Q3a75.txt')
H4Q3a75 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/H4Q3a75.txt')
MS1Q3a75 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/MS1Q3a75.txt')
MS1Q7a75 = ascii.read('/Users/kristophermortensen/Desktop/Kilonova Light Curve Data/MS1Q7a75.txt')
In [3]:
#Convert Absolute Magnitude to Apparent Magnitude
#app_mag: Merger Number -> Merger
#Converts the merger data from absolute magnitude to apparent magnitude.
def app_mag(merger, redshift):
merger['u']=app_mag_band(merger['u'],redshift)
merger['g']=app_mag_band(merger['g'],redshift)
merger['r']=app_mag_band(merger['r'],redshift)
merger['i']=app_mag_band(merger['i'],redshift)
merger['z']=app_mag_band(merger['z'],redshift)
return merger
#app_mag_band: ListofNumbers Number -> ListofNumbers
#converts all the absolute magnitudes into apparent magnitudes.
def app_mag_band(data, redshift):
return 5*np.log10(lumo_dist(redshift)/10)+data
#lumo_dist: Number -> Number
#converts redshift to luminosity distance
def lumo_dist(redshift):
return cosmo.luminosity_distance(redshift).to(astro.units.pc).value
#NS-NS Conversions:
APR1215_app=app_mag(APR1215, 0.015)
APR1314_app=app_mag(APR1314, 0.015)
H41215_app=app_mag(H41215, 0.015)
H41314_app=app_mag(H41314, 0.015)
Sly135_app=app_mag(Sly135, 0.015)
#NS-BH Conversions:
APR4Q3a75_app=app_mag(APR4Q3a75, 0.015)
H4Q3a75_app=app_mag(H4Q3a75, 0.015)
MS1Q3a75_app=app_mag(MS1Q3a75, 0.015)
MS1Q7a75_app=app_mag(MS1Q7a75, 0.015)
In [48]:
#Constants:
title = 28
subtitle = 20
axis = 24
ticksize = 16
legend = 16
#full_plot: Merger String String -> Image
#Takes the data from the Merger and plots the magnitudes and interpolations of each associated
# photometric band as a function of time measured in days.
def full_plot(merger,name):
fig=plt.figure(figsize=(12,11))
ax=fig.add_subplot(111)
abs_plotter1(merger)
abs_plot_interp(merger)
#LSST_plot(merger)
ax.axis([0,35,rounded_down(correct_max(merger)+3),(rounded_down(correct_min(merger))-2)])
ax.tick_params(labelsize=ticksize)
#plt.title("Light Curve of " + name, fontsize=title)
plt.xlabel("Time (Days)", fontsize=axis)
plt.ylabel('Magnitude', fontsize=axis)
plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
return
#abs_plotter1 Merger -> Image
#Plots the data from the ascii.reads of the merger.
def abs_plotter1(merger):
plot_band(merger, 'u', 'o', 'b')
plot_band(merger, 'g', 'o', 'g')
plot_band(merger, 'r', 'o', 'y')
plot_band(merger, 'i', 'o', 'orange')
plot_band(merger, 'z', 'o', 'r')
return
#plot_band: Merger String String String -> Image
#Plots the magnitudes of a photometric band in the merger.
def plot_band(merger, name, style, color):
plt.plot(merger['day'], merger[name], style, color=color, markeredgecolor='black', label=name)
#abs_plotter1 Merger -> Image
#Plots the interpolations calculated from data from the ascii.reads of the merger.
def abs_plot_interp(merger):
plot_interp(merger, 'u', '-', 'b', 0)
plot_interp(merger, 'g', '-', 'g', 1)
plot_interp(merger, 'r', '-', 'y', 2)
plot_interp(merger, 'i', '-', 'orange', 3)
plot_interp(merger, 'z', '-', 'r', 4)
return
#plot_band: Merger String String String Number -> Image
#Plots the interpolation of a photometric band in the merger.
def plot_interp(merger, name, style, color, number):
plt.plot(day_new(merger, name), interpolate(merger)[number](day_new(merger, name)),
style, color=color, label=name+' interp.')
return
#interpolate: Merger-> List of Interpolations
#Provides the necessary interpolation functions of the mergers in each associated photometric bands.
def interpolate(merger):
LSST_u=23.9
LSST_g=25.0
LSST_r=24.7
LSST_i=24.0
LSST_z=23.3
approx='cubic'
u=interp1d(merger['day'][np.isfinite(merger['u'])], merger['u'][np.isfinite(merger['u'])],
kind=approx, bounds_error=False, fill_value=LSST_u)
g=interp1d(merger['day'][np.isfinite(merger['g'])], merger['g'][np.isfinite(merger['g'])],
kind=approx, bounds_error=False, fill_value=LSST_g)
r=interp1d(merger['day'][np.isfinite(merger['r'])], merger['r'][np.isfinite(merger['r'])],
kind=approx, bounds_error=False, fill_value=LSST_r)
i=interp1d(merger['day'][np.isfinite(merger['i'])], merger['i'][np.isfinite(merger['i'])],
kind=approx, bounds_error=False, fill_value=LSST_i)
z=interp1d(merger['day'][np.isfinite(merger['z'])], merger['z'][np.isfinite(merger['z'])],
kind=approx, bounds_error=False, fill_value=LSST_z)
return [u, g, r, i, z]
#day_new Merger String -> List of Numbers
#Produces a new set of points between the range of the 'day' column in the merger data.
# This is useful for plotting the interpolations.
def day_new(merger, band):
set_length=100
day_new = np.linspace(min(merger['day'][np.isfinite(merger[band])]),
max(merger['day'][np.isfinite(merger[band])]),set_length)
return day_new
#LSST_plot: Merger -> Image
#Plots graphs of the LSST's 5-sigma depths in all five bands.
def LSST_plot(merger):
LSST_u=(merger['day']/merger['day'])-1+23.9
LSST_g=(merger['day']/merger['day'])-1+25.0
LSST_r=(merger['day']/merger['day'])-1+24.7
LSST_i=(merger['day']/merger['day'])-1+24.0
LSST_z=(merger['day']/merger['day'])-1+23.3
plt.plot(merger['day'], LSST_u, 'b', label='LSST \n u-band')
plt.plot(merger['day'], LSST_g, 'g', label='LSST \n g-band')
plt.plot(merger['day'], LSST_r, 'y', label='LSST \n r-band')
plt.plot(merger['day'], LSST_i, 'orange', label='LSST \n i-band')
plt.plot(merger['day'], LSST_z, 'r', label='LSST \n z-band')
return
#correct_max: Merger -> Number
#Outputs the largest magnitude in the entire merger data.
def correct_max(merger):
max_list = [max(merger['u'][np.isfinite(merger['u'])]),
max(merger['g'][np.isfinite(merger['g'])]),
max(merger['r'][np.isfinite(merger['r'])]),
max(merger['i'][np.isfinite(merger['i'])]),
max(merger['z'][np.isfinite(merger['z'])])]
return max(max_list)
#correct_min: Merger -> Number
#Outputs the smallest magnitude in the entire merger data.
def correct_min(merger):
min_list = [min(merger['u'][np.isfinite(merger['u'])]),
min(merger['g'][np.isfinite(merger['g'])]),
min(merger['r'][np.isfinite(merger['r'])]),
min(merger['i'][np.isfinite(merger['i'])]),
min(merger['z'][np.isfinite(merger['z'])])]
return min(min_list)
#round_up: number base -> number
#Rounds the number up to the nearest multiple of the base number.
def rounded_up(x, base=5):
return int(base * (round(float(x)/base)+1))
#round_down: number base -> number
#Rounds the number down to the nearest multiple of the base number.
def rounded_down(x, base=5):
return int(base * (round(float(x)/base)))
In [5]:
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
majorLocator = MultipleLocator(5)
majorFormatter = FormatStrFormatter('%d')
minorLocator = MultipleLocator(1)
bands=['u','g','r', 'i', 'z']
#rand_plot: Merger String String String -> Image
#Plots a randomly chosen band from the merger.
def band_plot(merger, name):
LSST_u=(merger['day']/merger['day'])-1+23.9
LSST_g=(merger['day']/merger['day'])-1+25.0
LSST_r=(merger['day']/merger['day'])-1+24.7
LSST_i=(merger['day']/merger['day'])-1+24.0
LSST_z=(merger['day']/merger['day'])-1+23.3
fig=plt.figure(figsize=(10,25))
for band in bands:
if band == 'u':
band_graph(merger, 'u', 'b', 511, 0, LSST_u)
elif band == 'g':
band_graph(merger, 'g', 'g', 512, 1, LSST_g)
elif band == 'r':
band_graph(merger, 'r', 'y', 513, 2, LSST_r)
elif band == 'i':
band_graph(merger, 'i', 'orange', 514, 3, LSST_i)
elif band == 'z':
band_graph(merger, 'z', 'r', 515, 4, LSST_z)
plt.suptitle("Light Curve Bands of " + name, fontsize=title, y=1.01)
plt.tight_layout()
return
def band_graph(merger, band, color, sub_num, num, LSST_band):
ax = plt.subplot(sub_num)
plot_band(merger, band, '.', color)
plot_interp(merger, band, '--', color, num)
plt.plot(merger['day'], LSST_band, color, label='LSST (' + band + ')')
plt.axis([0,35,rounded_up(correct_max(merger))+5,(rounded_down(correct_min(merger))-5)])
ax.xaxis.set_major_locator(majorLocator)
ax.xaxis.set_major_formatter(majorFormatter)
ax.xaxis.set_minor_locator(minorLocator)
plt.xlabel("Time (Days)", fontsize=axis)
plt.ylabel('Magnitude', fontsize=axis)
ax.tick_params(labelsize=ticksize)
plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
return
In [6]:
#The following are plots of only the data and the interpolations retrived from the data files in the ascii.reads,
# converted into apparent magnitude.
full_plot(APR1215_app, "APR-1215")
plt.close()
band_plot(APR1215_app, "APR-1215")
plt.close()
In [6]:
full_plot(APR1314_app, "APR4-1314")
plt.close()
band_plot(APR1314_app, "APR4-1314")
plt.close()
In [8]:
full_plot(H41215_app, "H4-1215")
plt.close()
band_plot(H41215_app, "H4-1215")
plt.close()
In [10]:
full_plot(H41314_app, "H4-1314")
plt.close()
band_plot(H41314_app, "H4-1314")
plt.close()
In [49]:
full_plot(Sly135, "Sly-135")
band_plot(Sly135_app, "Sly-135")
plt.close()
In [10]:
#The following are plots of only the data and the interpolations retrived from the data files in the ascii.reads,
# converted into apparent magnitude.
full_plot(APR4Q3a75_app, 'APR4Q3a75')
plt.close()
band_plot(APR4Q3a75_app, 'APR4Q3a75')
plt.close()
In [11]:
full_plot(H4Q3a75_app, 'H4Q3a75')
plt.close()
band_plot(H4Q3a75_app, 'H4Q3a75')
plt.close()
In [12]:
full_plot(MS1Q3a75_app, 'MS1Q3a75')
plt.close()
band_plot(MS1Q3a75_app, 'MS1Q3a75')
plt.close()
In [40]:
full_plot(MS1Q7a75, 'MS1Q7a75')
band_plot(MS1Q7a75_app, 'MS1Q7a75')
plt.close()
In [14]:
import random
#Constants:
LSST_u=23.9
LSST_g=25.0
LSST_r=24.7
LSST_i=24.0
LSST_z=23.3
trial_num=1000
epoch2=1
epoch3=2
bands=['u','g','r', 'i', 'z']
In [15]:
#Analyze: Merger -> Table
#Randomly runs 1,000 trials can places the information of the type of band, the distance, and the kilonova detection
# (0=fail; 1=success) into a table.
def Analyze(merger, epoch2, epoch3, trial_num, bands):
t_0 = [random.random()+1.06 for x in range(trial_num)]
band_list=[]
dist_list=[]
detect_list=[]
m0=[]
m1=[]
m7=[]
epoch_2 = [t+epoch2 for t in t_0]
epoch_3 = [t+epoch3 for t in t_0]
for t in t_0:
band=random.choice(bands)
band_list.append(band)
dist_list.append(rand_dist(merger))
m0.append(get_mag(merger, band, t))
m1.append(get_mag(merger, band, t+epoch2))
m7.append(get_mag(merger, band, t+epoch3))
if rand_detect(merger, band, t, epoch2, epoch3) == True:
detect_list.append(1)
else:
detect_list.append(0)
return Table([band_list, dist_list, t_0, epoch_2, epoch_3, m0, m1, m7, detect_list],
names=('Band', 'Distance', 'Epoch 1 (Days)', 'Epoch 2 (Days)', 'Epoch 3 (Days)',
'Magnitude\n(Epoch 1)', 'Magnitude\n(Epoch 2)', 'Magnitude\n(Epoch 3)','Kilonova Hit'))
#rand_dist: Merger -> Number
#Picks a random distance from 0 to 200Mpc.
def rand_dist(merger):
x=random.uniform(-1,1)
y=random.uniform(-1,1)
z=random.uniform(-1,1)
dist=np.sqrt(x**2+y**2+z**2)
if dist < 1:
return dist*200
else:
return rand_dist(merger)
return
#rand_detect: Merger String -> Boolean
#Determines whether a Kilonova is detected or not.
def rand_detect(merger, band, time, epoch2, epoch3):
if detect(merger, band, time) and detect(merger, band, time+epoch2) and not detect(merger, band, time+epoch3):
return True
else:
return False
return
#detect: Merger String Number -> Boolean
#Determines if LSST can detect the light curve of a specific band at a specific time.
def detect(merger, band, time):
if band == 'u':
mag_u = interpolate(merger)[0](time)+0.5
return in_view(mag_u, LSST_u)
elif band == 'g':
mag_g = interpolate(merger)[1](time)+0.5
return in_view(mag_g, LSST_g)
elif band == 'r':
mag_r = interpolate(merger)[2](time)+0.5
return in_view(mag_r, LSST_r)
elif band == 'i':
mag_i = interpolate(merger)[3](time)+0.5
return in_view(mag_i, LSST_i)
elif band == 'z':
mag_z = interpolate(merger)[4](time)+0.5
return in_view(mag_z, LSST_z)
return
#in_view: Number Number -> Boolean
#Checks if LSST can detect the light curve based on magnitude.
def in_view(mag, LSST):
if mag < LSST:
return True
else:
return False
return
def get_mag(merger, band, time):
if detect(merger, band, time)==True:
if band == 'u':
mag = interpolate(merger)[0](time)+0.5
elif band == 'g':
mag = interpolate(merger)[1](time)+0.5
elif band == 'r':
mag = interpolate(merger)[2](time)+0.5
elif band == 'i':
mag = interpolate(merger)[3](time)+0.5
elif band == 'z':
mag = interpolate(merger)[4](time)+0.5
elif detect(merger, band, time)==False:
if band == 'u':
mag = LSST_u
elif band == 'g':
mag = LSST_g
elif band == 'r':
mag = LSST_r
elif band == 'i':
mag = LSST_i
elif band == 'z':
mag = LSST_z
return mag
In [16]:
#NS-NS Mergers
APR1215_table =Analyze(APR1215_app, epoch2, epoch3, trial_num, bands)
In [17]:
APR1314_table = Analyze(APR1314_app, epoch2, epoch3, trial_num, bands)
In [18]:
H41215_table = Analyze(H41215_app, epoch2, epoch3, trial_num, bands)
In [19]:
H41314_table = Analyze(H41314_app, epoch2, epoch3, trial_num, bands)
In [20]:
Sly135_table = Analyze(Sly135_app, epoch2, epoch3, trial_num, bands)
In [21]:
#NS-BH Mergers
APR4Q3a75_table = Analyze(APR4Q3a75_app, epoch2, epoch3, trial_num, bands)
In [22]:
H4Q3a75_table = Analyze(H4Q3a75_app, epoch2, epoch3, trial_num, bands)
In [23]:
MS1Q3a75_table = Analyze(MS1Q3a75_app, epoch2, epoch3, trial_num, bands)
In [24]:
MS1Q7a75_table = Analyze(MS1Q7a75_app, epoch2, epoch3, trial_num, bands)
In [25]:
def total_efficiency(table, merger_base, epoch2, epoch3):
table[np.where(table['Band']=='u')] = band_efficiency(table[np.where(table['Band']=='u')], 'u',
merger_base, epoch2, epoch3, LSST_u)
table[np.where(table['Band']=='g')] = band_efficiency(table[np.where(table['Band']=='g')], 'g',
merger_base, epoch2, epoch3, LSST_g)
table[np.where(table['Band']=='r')] = band_efficiency(table[np.where(table['Band']=='r')], 'r',
merger_base, epoch2, epoch3, LSST_r)
table[np.where(table['Band']=='i')] = band_efficiency(table[np.where(table['Band']=='i')], 'i',
merger_base, epoch2, epoch3, LSST_i)
table[np.where(table['Band']=='z')] = band_efficiency(table[np.where(table['Band']=='z')], 'z',
merger_base, epoch2, epoch3, LSST_z)
return table
def band_efficiency(table, band, merger_base, epoch2, epoch3, LSST):
if all([item==LSST for item in table['Magnitude\n(Epoch 2)']]):
return table
elif any([item==LSST for item in table['Magnitude\n(Epoch 2)']]):
return table
else:
new_table=efficiency(table, band)[0]
if all([item==0 for item in new_table]):
epoch3 +=1
return band_efficiency(
Analyze(merger_base, epoch2, epoch3, len(table), [band]), band, merger_base, epoch2, epoch3, LSST)
elif all([item==100 or item==0 for item in new_table]):
return table
else:
epoch3 +=1
return band_efficiency(
Analyze(merger_base, epoch2, epoch3, len(table), [band]), band, merger_base, epoch2, epoch3, LSST)
#efficiency: Table String -> ListofNumbers
#Converts the number of detections into a fraction: number of succesful detections over the total number of
# observations at that distance.
def efficiency(merger_table, band):
countx, binx=np.histogram(np.array(filter_dist(filter_band(filter_hits(merger_table), band))),
bins=50, range=(0,200))
county, biny=np.histogram(np.array(filter_dist(filter_band(merger_table, band))), bins=50, range=(0,200))
np.seterr(divide='ignore', invalid='ignore')
counts=countx/county
for i in range(len(counts)):
if np.isnan(counts[i]):
counts[i]=0
return counts*100, binx
#filter_hits: Table -> Table
#Filters the table to only have data containing kilonova hits.
def filter_hits(table):
hits = np.where(table['Kilonova Hit'] == 1)[0]
return table[hits]
#filter_band: Table String -> Table
#Filters the table to only have data containing the designated band.
def filter_band(table, band):
band_filter=np.where(table['Band']==band)[0]
return table[band_filter]
#filter_hits: Table -> Table
#Filters the table to only have data containing the distance (Mpc).
def filter_dist(table):
return table['Distance']
In [26]:
#NS-NS Mergers
APR1215_eff = total_efficiency(APR1215_table, APR1215_app, epoch2, epoch3)
In [27]:
APR1314_eff = total_efficiency(APR1314_table, APR1314_app, epoch2, epoch3)
In [28]:
H41215_eff = total_efficiency(H41215_table, H41215_app, epoch2, epoch3)
In [29]:
H41314_eff = total_efficiency(H41314_table, H41314_app, epoch2, epoch3)
In [30]:
Sly135_eff = total_efficiency(Sly135_table, Sly135_app, epoch2, epoch3)
In [31]:
#NS-BH Mergers
APR4Q3a75_eff = total_efficiency(APR4Q3a75_table, APR4Q3a75_app, epoch2, epoch3)
In [32]:
H4Q3a75_eff = total_efficiency(H4Q3a75_table, H4Q3a75_app, epoch2, epoch3)
In [33]:
MS1Q3a75_eff = total_efficiency(MS1Q3a75_table, MS1Q3a75_app, epoch2, epoch3)
In [34]:
MS1Q7a75_eff = total_efficiency(MS1Q7a75_table, MS1Q7a75_app, epoch2, epoch3)
In [35]:
def max_day(table):
u=max_day_band(filter_band(table, 'u'))
g=max_day_band(filter_band(table, 'g'))
r=max_day_band(filter_band(table, 'r'))
i=max_day_band(filter_band(table, 'i'))
z=max_day_band(filter_band(table, 'z'))
return [u, g, r, i, z]
def max_day_band(table):
day = np.subtract(table['Epoch 3 (Days)'], table['Epoch 1 (Days)'])
return day[0]
In [111]:
#plot_hist: Table String -> Image
#Plots a histogram of the number of Kilonova detections for each range of distance (Mpc).
def plot_hist(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(511)
ax2 = plt.subplot(512)
ax3 = plt.subplot(513)
ax4 = plt.subplot(514)
ax5 = plt.subplot(515)
#_______________________________________________________________________________________________________
ax1.hist(np.array(filter_dist(filter_band(filter_hits(table), 'u'))),
bins=50, range=(0,200), facecolor='blue', edgecolor='black', label='u')
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax1.set_xlabel('Distance (Mpc)', fontsize=axis)
ax1.set_ylabel('Discoveries', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax2.hist(np.array(filter_dist(filter_band(filter_hits(table), 'g'))),
bins=50, range=(0,200), facecolor='green', edgecolor='black', label='g')
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax2.set_xlabel('Distance (Mpc)', fontsize=axis)
ax2.set_ylabel('Discoveries', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax3.hist(np.array(filter_dist(filter_band(filter_hits(table), 'r'))),
bins=50, range=(0,200), facecolor='yellow', edgecolor='black', label='r')
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax3.set_xlabel('Distance (Mpc)', fontsize=axis)
ax3.set_ylabel('Discoveries', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax4.hist(np.array(filter_dist(filter_band(filter_hits(table), 'i'))),
bins=50, range=(0,200), facecolor='orange', edgecolor='black', label='i')
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax4.set_xlabel('Distance (Mpc)', fontsize=axis)
ax4.set_ylabel('Discoveries', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax5.hist(np.array(filter_dist(filter_band(filter_hits(table), 'z'))),
bins=50, range=(0,200), facecolor='red', edgecolor='black', label='z')
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax5.set_xlabel('Distance (Mpc)', fontsize=axis)
ax5.set_ylabel('Discoveries', fontsize=axis)
ax5.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05)
plt.tight_layout()
return
#filter_hits: Table -> Table
#Filters the table to only have data containing kilonova hits.
def filter_hits(table):
hits = np.where(table['Kilonova Hit'] == 1)[0]
return table[hits]
#filter_band: Table String -> Table
#Filters the table to only have data containing the designated band.
def filter_band(table, band):
band_filter=np.where(table['Band']==band)[0]
return table[band_filter]
#filter_hits: Table -> Table
#Filters the table to only have data containing the distance (Mpc).
def filter_dist(table):
return table['Distance']
from matplotlib import gridspec
def plot_dimension(table, band, number):
gs = gridspec.GridSpec(5, 5)
new_table=efficiency(table, band)
if all([item==0 for item in new_table]):
ax = plt.subplot(gs[number, 0])
else:
ax = plt.subplot(gs[number, :])
return ax
def epoch_table(table):
return Table([bands, max_day(table)], names=('Bands', 'Epoch 3 \n (Days)'))
In [112]:
plot_hist(APR1215_eff, "APR4-1215")
epoch_table(APR1215_eff)
Out[112]:
In [113]:
plot_hist(APR1314_eff, "APR4-1314")
epoch_table(APR1314_eff)
plt.close()
In [114]:
plot_hist(H41215_eff, "H4-1215")
epoch_table(H41215_eff)
plt.close()
In [115]:
plot_hist(H41314_eff, "H4-1314")
epoch_table(H41314_eff)
plt.close()
In [116]:
plot_hist(Sly135_eff, "Sly-135")
epoch_table(Sly135_eff)
plt.close()
In [117]:
plot_hist(APR4Q3a75_eff, "APR4Q3a75")
plt.close()
epoch_table(APR4Q3a75_eff)
plt.close()
In [118]:
plot_hist(H4Q3a75_eff, "H4Q3a75")
epoch_table(H4Q3a75_eff)
plt.close()
In [119]:
plot_hist(MS1Q3a75_eff, "MS1Q3a75")
epoch_table(MS1Q3a75_eff)
plt.close()
In [120]:
plot_hist(MS1Q7a75_eff, "MS1Q7a75")
epoch_table(MS1Q7a75_eff)
plt.close()
In [46]:
def box_plot(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(511)
ax2 = plt.subplot(512)
ax3 = plt.subplot(513)
ax4 = plt.subplot(514)
ax5 = plt.subplot(515)
#_______________________________________________________________________________________________________
if len(np.array(filter_dist(filter_band(filter_hits(table), 'u')))) == 0:
ax1.hist(np.array(filter_dist(filter_band(filter_hits(table), 'u'))))
else:
ax1.boxplot(np.array(filter_dist(filter_band(filter_hits(table), 'u'))), 0, 'rs', 0)
ax1.set_title('u-band', fontsize=subtitle)
ax1.set_xlabel('Distance (Mpc)', fontsize=axis)
ax1.set_ylabel('Discoveries', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
if len(np.array(filter_dist(filter_band(filter_hits(table), 'g')))) == 0:
ax2.hist(np.array(filter_dist(filter_band(filter_hits(table), 'g'))))
else:
ax2.boxplot(np.array(filter_dist(filter_band(filter_hits(table), 'g'))), 0, 'rs', 0)
ax2.set_title('g-band', fontsize=subtitle)
ax2.set_xlabel('Distance (Mpc)', fontsize=axis)
ax2.set_ylabel('Discoveries', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
if len(np.array(filter_dist(filter_band(filter_hits(table), 'r')))) == 0:
ax3.hist(np.array(filter_dist(filter_band(filter_hits(table), 'r'))))
else:
ax3.boxplot(np.array(filter_dist(filter_band(filter_hits(table), 'r'))), 0, 'rs', 0)
ax3.set_title('r-band', fontsize=subtitle)
ax3.set_xlabel('Distance (Mpc)', fontsize=axis)
ax3.set_ylabel('Discoveries', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
if len(np.array(filter_dist(filter_band(filter_hits(table), 'i')))) == 0:
ax4.hist(np.array(filter_dist(filter_band(filter_hits(table), 'i'))))
else:
ax4.boxplot(np.array(filter_dist(filter_band(filter_hits(table), 'i'))), 0, 'rs', 0)
ax4.set_title('i-band', fontsize=subtitle)
ax4.set_xlabel('Distance (Mpc)', fontsize=axis)
ax4.set_ylabel('Discoveries', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
if len(np.array(filter_dist(filter_band(filter_hits(table), 'z')))) == 0:
ax5.hist(np.array(filter_dist(filter_band(filter_hits(table), 'z'))))
else:
ax5.boxplot(np.array(filter_dist(filter_band(filter_hits(table), 'z'))), 0, 'rs', 0)
ax5.set_title('z-band', fontsize=subtitle)
ax5.set_xlabel('Distance (Mpc)', fontsize=axis)
ax5.set_ylabel('Discoveries', fontsize=axis)
ax5.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05)
plt.tight_layout()
return
In [122]:
#plot_bar: Table String -> Image
#Plots a bar graph of the efficiency of the LSST as a function of distance(Mpc).
def plot_bar(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(511)
ax2 = plt.subplot(512)
ax3 = plt.subplot(513)
ax4 = plt.subplot(514)
ax5 = plt.subplot(515)
#_______________________________________________________________________________________________________
center1 = (efficiency(table, 'u')[1][:-1] + efficiency(table, 'u')[1][1:]) / 2
width1 = (efficiency(table, 'u')[1][1] - efficiency(table, 'u')[1][0])
ax1.bar(center1, efficiency(table, 'u')[0],
width=width1, align='center', facecolor='blue', edgecolor='black', label='u')
ax1.tick_params(labelsize=ticksize)
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax1.set_xlabel('Distance (Mpc)', fontsize=axis)
ax1.set_ylabel('Efficiency(%)', fontsize=axis)
ax1.set_ylim(0,100)
#_______________________________________________________________________________________________________
center2 = (efficiency(table, 'g')[1][:-1] + efficiency(table, 'g')[1][1:]) / 2
width2 = (efficiency(table, 'g')[1][1] - efficiency(table, 'g')[1][0])
ax2.bar(center2, efficiency(table, 'g')[0],
width=width2, align='center', facecolor='green', edgecolor='black', label='g')
ax2.tick_params(labelsize=ticksize)
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax2.set_xlabel('Distance (Mpc)', fontsize=axis)
ax2.set_ylabel('Efficiency(%)', fontsize=axis)
ax2.set_ylim(0,100)
#_______________________________________________________________________________________________________
center3 = (efficiency(table, 'r')[1][:-1] + efficiency(table, 'r')[1][1:]) / 2
width3 = (efficiency(table, 'r')[1][1] - efficiency(table, 'r')[1][0])
ax3.bar(center3, efficiency(table, 'r')[0],
width=width3, align='center', facecolor='yellow', edgecolor='black', label='r')
ax3.tick_params(labelsize=ticksize)
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax3.set_xlabel('Distance (Mpc)', fontsize=axis)
ax3.set_ylabel('Efficiency(%)', fontsize=axis)
ax3.set_ylim(0,100)
#_______________________________________________________________________________________________________
center4 = (efficiency(table, 'i')[1][:-1] + efficiency(table, 'i')[1][1:]) / 2
width4 = (efficiency(table, 'i')[1][1] - efficiency(table, 'i')[1][0])
ax4.bar(center4, efficiency(table, 'i')[0],
width=width4, align='center', facecolor='orange', edgecolor='black', label='i')
ax4.tick_params(labelsize=ticksize)
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax4.set_xlabel('Distance (Mpc)', fontsize=axis)
ax4.set_ylabel('Efficiency(%)', fontsize=axis)
ax4.set_ylim(0,100)
#_______________________________________________________________________________________________________
center5 = (efficiency(table, 'z')[1][:-1] + efficiency(table, 'z')[1][1:]) / 2
width5 = (efficiency(table, 'z')[1][1] - efficiency(table, 'z')[1][0])
ax5.bar(center5, efficiency(table, 'z')[0],
width=width5, align='center', facecolor='blue', edgecolor='black', label='z')
ax5.tick_params(labelsize=ticksize)
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax5.set_xlabel('Distance (Mpc)', fontsize=axis)
ax5.set_ylabel('Efficiency(%)', fontsize=axis)
ax5.set_ylim(0,100)
#_______________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05)
plt.tight_layout()
return
In [125]:
plot_bar(APR1215_eff, "APR4-1215")
epoch_table(APR1215_eff)
Out[125]:
In [51]:
plot_bar(APR1314_eff, "APR4-1314")
epoch_table(APR1314_eff)
plt.close()
In [52]:
plot_bar(H41215_eff, "H4-1215")
epoch_table(H41215_eff)
plt.close()
In [53]:
plot_bar(H41314_eff, "H4-1314")
epoch_table(H41314_eff)
plt.close()
In [54]:
plot_bar(Sly135_eff, "Sly-135")
epoch_table(Sly135_eff)
plt.close()
In [55]:
plot_bar(APR4Q3a75_eff, "APR4Q3a75")
plt.close()
epoch_table(APR4Q3a75_eff)
plt.close()
In [56]:
plot_bar(H4Q3a75_eff, "H4Q3a75")
epoch_table(H4Q3a75_eff)
plt.close()
In [57]:
plot_bar(MS1Q3a75_eff, "MS1Q3a75")
epoch_table(MS1Q3a75_eff)
plt.close()
In [58]:
plot_bar(MS1Q7a75_eff, "MS1Q7a75")
epoch_table(MS1Q7a75_eff)
plt.close()
In [59]:
#filter_hits: Table -> Table
#Filters the table to only have data containing kilonova hits.
def filter_hits(table):
hits = np.where(table['Kilonova Hit'] == 1)[0]
return table[hits]
#filter_band: Table String -> Table
#Filters the table to only have data containing the designated band.
def filter_band(table, band):
band_filter=np.where(table['Band']==band)[0]
return table[band_filter]
def filter_mags(table):
return table['Magnitude\n(Epoch 1)', 'Magnitude\n(Epoch 2)', 'Magnitude\n(Epoch 3)']
def delta_m1_band(table, band):
m1 = np.subtract(filter_mags(filter_band(table, band))['Magnitude\n(Epoch 2)'],
filter_mags(filter_band(table, band))['Magnitude\n(Epoch 1)'])
return m1
def plot_m1(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(511)
ax2 = plt.subplot(512)
ax3 = plt.subplot(513)
ax4 = plt.subplot(514)
ax5 = plt.subplot(515)
#_______________________________________________________________________________________________________
ax1.hist(np.array(delta_m1_band(table, 'u')),
bins=75, range=(np.floor(min_delta(table)), np.ceil(max_delta(table))),
facecolor='blue', edgecolor='black', label='u')
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax1.set_ylabel('Frequency', fontsize=axis)
ax1.set_xlabel('$Δm_{1}$', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax2.hist(np.array(delta_m1_band(table, 'g')),
bins=75, range=(np.floor(min_delta(table)), np.ceil(max_delta(table))),
facecolor='green', edgecolor='black', label='g')
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax2.set_ylabel('Frequency', fontsize=axis)
ax2.set_xlabel('$Δm_{1}$', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax3.hist(np.array(delta_m1_band(table, 'r')),
bins=75, range=(np.floor(min_delta(table)), np.ceil(max_delta(table))),
facecolor='yellow', edgecolor='black', label='r')
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax3.set_ylabel('Frequency', fontsize=axis)
ax3.set_xlabel('$Δm_{1}$', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax4.hist(np.array(delta_m1_band(table, 'i')),
bins=75, range=(np.floor(min_delta(table)), np.ceil(max_delta(table))),
facecolor='orange', edgecolor='black', label='i')
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax4.set_ylabel('Frequency', fontsize=axis)
ax4.set_xlabel('$Δm_{1}$', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax5.hist(np.array(delta_m1_band(table, 'z')),
bins=75, range=(np.floor(min_delta(table)), np.ceil(max_delta(table))),
facecolor='red', edgecolor='black', label='z')
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax5.set_ylabel('Frequency', fontsize=axis)
ax5.set_xlabel('$Δm_{1}$', fontsize=axis)
ax5.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05)
plt.tight_layout()
return
def max_delta(table):
u = check_empty(np.array(delta_m1_band(table, 'u')))
g = check_empty(np.array(delta_m1_band(table, 'g')))
r = check_empty(np.array(delta_m1_band(table, 'r')))
i = check_empty(np.array(delta_m1_band(table, 'i')))
z = check_empty(np.array(delta_m1_band(table, 'z')))
return max([u, g, r, i, z])
def check_empty(array):
if len(array) == 0:
max_num = 0
else:
max_num = max(array)
return max_num
def min_delta(table):
u = check_empty_min(np.array(delta_m1_band(table, 'u')))
g = check_empty_min(np.array(delta_m1_band(table, 'g')))
r = check_empty_min(np.array(delta_m1_band(table, 'r')))
i = check_empty_min(np.array(delta_m1_band(table, 'i')))
z = check_empty_min(np.array(delta_m1_band(table, 'z')))
return min([u, g, r, i, z])
def check_empty_min(array):
if len(array) == 0:
min_num = 0
else:
min_num = min(array)
return min_num
In [60]:
plot_m1(APR1215_eff, 'APR4-1215')
plt.close()
In [61]:
plot_m1(APR1314_eff, 'APR4-1314')
plt.close()
In [62]:
plot_m1(H41215_eff, "H4-1215")
plt.close()
In [63]:
plot_m1(H41314_eff, "H4-1314")
plt.close()
In [64]:
plot_m1(Sly135_eff, "Sly-135")
plt.close()
In [65]:
plot_m1(APR4Q3a75_eff, "APR4Q3a75")
plt.close()
In [66]:
plot_m1(H4Q3a75_eff, "H4Q3a75")
plt.close()
In [67]:
plot_m1(MS1Q3a75_eff, "MS1Q3a75")
plt.close()
In [68]:
plot_m1(MS1Q7a75_table, "MS1Q7a75")
plt.close()
In [69]:
#filter_hits: Table -> Table
#Filters the table to only have data containing kilonova hits.
def filter_hits(table):
hits = np.where(table['Kilonova Hit'] == 1)[0]
return table[hits]
#filter_band: Table String -> Table
#Filters the table to only have data containing the designated band.
def filter_band(table, band):
band_filter=np.where(table['Band']==band)[0]
return table[band_filter]
def filter_mags(table):
return table['Magnitude\n(Epoch 1)', 'Magnitude\n(Epoch 2)', 'Magnitude\n(Epoch 3)']
def delta_m_epoch_band(table, band):
m1 = np.subtract(filter_mags(filter_band(table, band))['Magnitude\n(Epoch 3)'],
filter_mags(filter_band(table, band))['Magnitude\n(Epoch 1)'])
return m1
def max_delta7(table):
u = check_empty(np.array(delta_m_epoch_band(table, 'u')))
g = check_empty(np.array(delta_m_epoch_band(table, 'g')))
r = check_empty(np.array(delta_m_epoch_band(table, 'r')))
i = check_empty(np.array(delta_m_epoch_band(table, 'i')))
z = check_empty(np.array(delta_m_epoch_band(table, 'z')))
return max([u, g, r, i, z])
def check_empty(array):
if len(array) == 0:
max_num = 0
else:
max_num = max(array)
return max_num
def plot_m_epoch(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(511)
ax2 = plt.subplot(512)
ax3 = plt.subplot(513)
ax4 = plt.subplot(514)
ax5 = plt.subplot(515)
#_______________________________________________________________________________________________________
ax1.hist(np.array(delta_m_epoch_band(table, 'u')),
bins=75, range=(0, np.ceil(max_delta7(table))), facecolor='blue', edgecolor='black', label='u')
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax1.set_ylabel('Frequency', fontsize=14)
ax1.set_xlabel('$Δm_{'+str(int(round(max_day(table)[0])))+'}$', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax2.hist(np.array(delta_m_epoch_band(table, 'g')),
bins=75, range=(0, np.ceil(max_delta7(table))), facecolor='green', edgecolor='black', label='g')
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax2.set_ylabel('Frequency', fontsize=14)
ax2.set_xlabel('$Δm_{'+str(int(round(max_day(table)[1])))+'}$', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax3.hist(np.array(delta_m_epoch_band(table, 'r')),
bins=75, range=(0, np.ceil(max_delta7(table))), facecolor='yellow', edgecolor='black', label='r')
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax3.set_ylabel('Frequency', fontsize=14)
ax3.set_xlabel('$Δm_{'+str(int(round(max_day(table)[2])))+'}$', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax4.hist(np.array(delta_m_epoch_band(table, 'i')),
bins=75, range=(0, np.ceil(max_delta7(table))), facecolor='orange', edgecolor='black', label='i')
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax4.set_ylabel('Frequency', fontsize=14)
ax4.set_xlabel('$Δm_{'+str(int(round(max_day(table)[3])))+'}$', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
ax5.hist(np.array(delta_m_epoch_band(table, 'z')),
bins=75, range=(0, np.ceil(max_delta7(table))), facecolor='red', edgecolor='black', label='z')
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax5.set_ylabel('Frequency', fontsize=14)
ax5.set_xlabel('$Δm_{'+str(int(round(max_day(table)[4])))+'}$', fontsize=axis)
ax5.tick_params(labelsize=ticksize)
#_______________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05, ha = 'center')
plt.tight_layout()
return
In [70]:
plot_m_epoch(APR1215_eff, 'APR4-1215')
plt.close()
In [71]:
plot_m_epoch(APR1314_eff, 'APR4-1314')
plt.close()
In [72]:
plot_m_epoch(H41215_eff, "H4-1215")
plt.close()
In [73]:
plot_m_epoch(H41314_eff, "H4-1314")
plt.close()
In [74]:
plot_m_epoch(Sly135_eff, "Sly-135")
plt.close()
In [75]:
plot_m_epoch(APR4Q3a75_eff, "APR4Q3a75")
plt.close()
In [76]:
plot_m_epoch(H4Q3a75_eff, "H4Q3a75")
plt.close()
In [77]:
plot_m_epoch(MS1Q3a75_eff, "MS1Q3a75")
plt.close()
In [78]:
plot_m_epoch(MS1Q7a75_eff, "MS1Q7a75")
plt.close()
In [79]:
trial_num = 1000
def Analyze2(merger):
time = weather_output(merger)
colors = measure_color(merger, time)
discoveries = []
distance = []
for i in range(len(colors['Color (Epoch 1)'])):
distance.append(rand_dist(merger))
if colors['Color (Epoch 1)'][i] > 0.5 and colors['Color (Epoch 2)'][i] > 0.5 and \
colors['$\Delta m_{3}$ (z-band)'][i] > 1.0:
discoveries.append(1)
else:
discoveries.append(0)
return Table([colors['Band'], distance, discoveries], names=('Band', 'Distance (Mpc)', 'Discovery'))
def measure_color(merger, time0):
bands = ['u', 'g', 'r', 'i']
time = time0
color_list1 = []
color_list2 = []
band_list = []
delta_z = []
for t in time['Epoch 1']:
band = random.choice(bands)
color1 = np.subtract(get_mag(merger, band, t), get_mag(merger, 'z', t))
color_list1.append(color1)
band_list.append(band+'-z')
color2 = np.subtract(get_mag(merger, band, float(time['Epoch 2'][np.where(time['Epoch 1']==t)])),
get_mag(merger, 'z', float(time['Epoch 2'][np.where(time['Epoch 1']==t)])))
color_list2.append(color2)
z = np.subtract(get_mag(merger, 'z', float(time['Epoch 2'][np.where(time['Epoch 1']==t)])),
get_mag(merger, 'z', t))
delta_z.append(z)
return Table([band_list, color_list1, color_list2, delta_z],
names=('Band', 'Color (Epoch 1)', 'Color (Epoch 2)', '$\Delta m_{3}$ (z-band)'))
def weather_output(merger):
time = np.empty((trial_num,4))
for i in range(len(time[:,1])):
t = random.random()+1.06
t1 = weather_simulator(t)
t2 = weather_simulator(t1+3)
t3 = weather_simulator(t2+3)
t4 = weather_simulator(t3+3)
time[i,0] = t1
time[i,1] = t2
time[i,2] = t3
time[i,3] = t4
return Table([time[:,0], time[:,1], time[:,2], time[:,3]], names=('Epoch 1', 'Epoch 2', 'Epoch 3', 'Epoch 4'))
def weather_simulator(time):
weather = random.random()
if weather < 0.8:
time = time
else:
time = weather_simulator(time+1)
return time
#rand_detect: Merger String -> Boolean
#Determines whether a Kilonova is detected or not.
def rand_detect2(merger, band, time1, time2):
if detect(merger, band, time1) and detect(merger, band, time2):
return True
else:
return False
return
#detect: Merger String Number -> Boolean
#Determines if LSST can detect the light curve of a specific band at a specific time.
def detect(merger, band, time):
LSST_u=23.9
LSST_g=25.0
LSST_r=24.7
LSST_i=24.0
LSST_z=23.3
if band == 'u':
mag_u = interpolate(merger)[0](time)+0.5
return in_view(mag_u, LSST_u)
elif band == 'g':
mag_g = interpolate(merger)[1](time)+0.5
return in_view(mag_g, LSST_g)
elif band == 'r':
mag_r = interpolate(merger)[2](time)+0.5
return in_view(mag_r, LSST_r)
elif band == 'i':
mag_i = interpolate(merger)[3](time)+0.5
return in_view(mag_i, LSST_i)
elif band == 'z':
mag_z = interpolate(merger)[4](time)+0.5
return in_view(mag_z, LSST_z)
return
#in_view: Number Number -> Boolean
#Checks if LSST can detect the light curve based on magnitude.
def in_view(mag, LSST):
if mag < LSST:
return True
else:
return False
return
In [80]:
#NS-NS Mergers
APR1215_table2 =Analyze2(APR1215_app)
In [81]:
APR1314_table2 = Analyze2(APR1314_app)
In [82]:
H41215_table2 = Analyze2(H41215_app)
In [83]:
H41314_table2 = Analyze2(H41314_app)
In [84]:
Sly135_table2 = Analyze2(Sly135_app)
In [85]:
#NS-BH Mergers
APR4Q3a75_table2 = Analyze2(APR4Q3a75_app)
In [86]:
H4Q3a75_table2 = Analyze2(H4Q3a75_app)
In [87]:
MS1Q3a75_table2 = Analyze2(MS1Q3a75_app)
In [88]:
MS1Q7a75_table2 = Analyze2(MS1Q7a75_app)
In [89]:
#filter_hits: Table -> Table
#Filters the table to only have data containing kilonova hits.
def filter_hits2(table):
hits = np.where(table['Discovery'] == 1)[0]
return table[hits]
#filter_band: Table String -> Table
#Filters the table to only have data containing the designated band.
def filter_band2(table, band):
band_filter=np.where(table['Band']==band)[0]
return table[band_filter]
#filter_hits: Table -> Table
#Filters the table to only have data containing the distance (Mpc).
def filter_dist2(table):
return table['Distance (Mpc)']
#plot_hist: Table String -> Image
#Plots a histogram of the number of Kilonova detections for each range of distance (Mpc).
def plot_hist2(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(411)
ax2 = plt.subplot(412)
ax3 = plt.subplot(413)
ax4 = plt.subplot(414)
#________________________________________________________________________________________________
n1, bins1 = np.histogram(np.array(filter_dist2(filter_band2(filter_hits2(table), 'u-z'))),
bins = 50, range=(0,200))
n1b, bins1b = np.histogram(np.array(filter_dist2(filter_band2(table, 'u-z'))),
bins = 50, range=(0,200))
center1 = (bins1[:-1] + bins1[1:]) / 2
width1 = (bins1[1] - bins1[0])
ax1.bar(center1, n1, facecolor='blue', edgecolor='black',
align='center', width=width1, label='Discoveries')
ax1.plot(center1, n1b, 'o-', color='orange', markeredgecolor='black', label='Total Observations')
handles1,labels1 = ax1.get_legend_handles_labels()
handles1 =handles1[::-1]
labels1 = labels1[::-1]
ax1.legend(handles1, labels1, bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax1.set_title('u $-$ z', fontsize=subtitle)
ax1.set_xlabel('Distance (Mpc)', fontsize=axis)
ax1.set_ylabel('Discoveries', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
ax1.set_ylim(0,25)
#________________________________________________________________________________________________
n2, bins2 = np.histogram(np.array(filter_dist2(filter_band2(filter_hits2(table), 'g-z'))),
bins = 50, range=(0,200))
n2b, bins2b = np.histogram(np.array(filter_dist2(filter_band2(table, 'g-z'))),
bins = 50, range=(0,200))
center2 = (bins2[:-1] + bins2[1:]) / 2
width2 = (bins2[1] - bins2[0])
ax2.bar(center2, n2, facecolor='green', edgecolor='black',
align='center', width=width2, label='Discoveries')
ax2.plot(center2, n2b, 'o-', color='red', markeredgecolor='black', label='Total Observations')
handles2,labels2 = ax2.get_legend_handles_labels()
handles2 =handles2[::-1]
labels2 = labels2[::-1]
ax2.legend(handles2, labels2, bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax2.set_title('g $-$ z', fontsize=subtitle)
ax2.set_xlabel('Distance (Mpc)', fontsize=axis)
ax2.set_ylabel('Discoveries', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
ax2.set_ylim(0,25)
#________________________________________________________________________________________________
n3, bins3 = np.histogram(np.array(filter_dist2(filter_band2(filter_hits2(table), 'r-z'))),
bins = 50, range=(0,200))
n3b, bins3b = np.histogram(np.array(filter_dist2(filter_band2(table, 'r-z'))),
bins = 50, range=(0,200))
center3 = (bins3[:-1] + bins3[1:]) / 2
width3 = (bins3[1] - bins3[0])
ax3.bar(center3, n3, facecolor='yellow', edgecolor='black',
align='center', width=width3, label='Discoveries')
ax3.plot(center3, n3b, 'o-', color='violet', markeredgecolor='black', label='Total Observations')
handles3,labels3 = ax3.get_legend_handles_labels()
handles3 =handles3[::-1]
labels3 = labels3[::-1]
ax3.legend(handles3, labels3, bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax3.set_title('r $-$ z', fontsize=subtitle)
ax3.set_xlabel('Distance (Mpc)', fontsize=axis)
ax3.set_ylabel('Discoveries', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
ax3.set_ylim(0,25)
#________________________________________________________________________________________________
n4, bins4 = np.histogram(np.array(filter_dist2(filter_band2(filter_hits2(table), 'i-z'))),
bins = 50, range=(0,200))
n4b, bins4b = np.histogram(np.array(filter_dist2(filter_band2(table, 'i-z'))),
bins = 50, range=(0,200))
center4 = (bins4[:-1] + bins4[1:]) / 2
width4 = (bins4[1] - bins4[0])
ax4.bar(center4, n4, facecolor='orange', edgecolor='black',
align='center', width=width4, label='Discoveries')
ax4.plot(center4, n4b, 'o-', color='turquoise', markeredgecolor='black', label='Total Observations')
handles4,labels4 = ax4.get_legend_handles_labels()
handles4 =handles4[::-1]
labels4 = labels4[::-1]
ax4.legend(handles4, labels4, bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0., prop={'size':legend})
ax4.set_title('i $-$ z', fontsize=subtitle)
ax4.set_xlabel('Distance (Mpc)', fontsize=axis)
ax4.set_ylabel('Discoveries', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
ax4.set_ylim(0,25)
#________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05, horizontalalignment='center')
plt.tight_layout()
return
In [90]:
plot_hist2(APR1215_table2, 'APR4-1215')
In [91]:
plot_hist2(APR1314_table2, 'APR4-1314')
plt.close()
In [92]:
plot_hist2(H41215_table2, 'H4-1215')
plt.close()
In [93]:
plot_hist2(H41314_table2, 'H4-1314')
plt.close()
In [94]:
plot_hist2(Sly135_table2, 'Sly-135')
In [95]:
plot_hist2(APR4Q3a75_table2, 'APR4Q3a75')
In [96]:
plot_hist2(H4Q3a75_table2, 'H4Q3a75')
plt.close()
In [97]:
plot_hist2(MS1Q3a75_table2, 'MS1Q3a75')
plt.close()
In [98]:
plot_hist2(MS1Q7a75_table2, 'MS1Q7a75')
plt.close()
In [99]:
#efficiency: Table String -> ListofNumbers
#Converts the number of detections into a fraction: number of succesful detections over the total number of
# observations at that distance.
def efficiency2(merger_table, band):
countx, binx=np.histogram(np.array(filter_dist2(filter_band2(filter_hits2(merger_table), band))),
bins=50, range=(0,200))
county, biny=np.histogram(np.array(filter_dist2(filter_band2(merger_table, band))), bins=50, range=(0,200))
np.seterr(divide='ignore', invalid='ignore')
counts=countx/county
for i in range(len(counts)):
if np.isnan(counts[i]):
counts[i]=0
return counts*100, binx
#plot_bar: Table String -> Image
#Plots a bar graph of the efficiency of the LSST as a function of distance(Mpc).
def plot_bar2(table, name):
plt.figure(figsize=(10,11))
ax1 = plt.subplot(411)
ax1b = ax1.twinx()
ax2 = plt.subplot(412)
ax2b = ax2.twinx()
ax3 = plt.subplot(413)
ax3b = ax3.twinx()
ax4 = plt.subplot(414)
ax4b = ax4.twinx()
#________________________________________________________________________________________________
n1b, bins1b = np.histogram(np.array(filter_dist2(filter_band2(table, 'u-z'))),
bins = 50, range=(0,200))
center1 = (efficiency2(table, 'u-z')[1][:-1] + efficiency2(table, 'u-z')[1][1:]) / 2
width1 = (efficiency2(table, 'u-z')[1][1] - efficiency2(table, 'u-z')[1][0])
ax1.bar(center1, efficiency2(table, 'u-z')[0],
width=width1, align='center', facecolor='blue', edgecolor='black', label='Efficiency')
ax1b.plot(center1, n1b, 'o-', color='orange', markeredgecolor='black', label='Total Observations')
handles1,labels1 = ax1b.get_legend_handles_labels()
handles1 =handles1[::-1]
labels1 = labels1[::-1]
ax1.legend(handles1, labels1, bbox_to_anchor=(1.05, 1), loc=3, borderaxespad=0., prop={'size':12})
ax1.set_title('u $-$ z', fontsize=subtitle)
ax1.set_xlabel('Distance (Mpc)', fontsize=axis)
ax1.set_ylabel('Efficiency ($\%$)', fontsize=axis)
ax1.tick_params(labelsize=ticksize)
ax1b.tick_params(labelsize=ticksize)
ax1.set_ylim(0,100)
ax1b.set_ylim(0,25)
ax1b.set_ylabel('Discoveries', fontsize=axis)
#________________________________________________________________________________________________
n2b, bins2b = np.histogram(np.array(filter_dist2(filter_band2(table, 'g-z'))),
bins = 50, range=(0,200))
center2 = (efficiency2(table, 'g-z')[1][:-1] + efficiency2(table, 'g-z')[1][1:]) / 2
width2 = (efficiency2(table, 'g-z')[1][1] - efficiency2(table, 'g-z')[1][0])
ax2.bar(center2, efficiency2(table, 'g-z')[0],
width=width2, align='center', facecolor='green', edgecolor='black', label='Efficiency')
ax2b.plot(center2, n2b, 'o-', color='red', markeredgecolor='black', label='Total Observations')
handles2,labels2 = ax2b.get_legend_handles_labels()
handles2 =handles2[::-1]
labels2 = labels2[::-1]
ax2.legend(handles2, labels2, bbox_to_anchor=(1.05, 1), loc=3, borderaxespad=0., prop={'size':12})
ax2.set_title('g $-$ z', fontsize=subtitle)
ax2.set_xlabel('Distance (Mpc)', fontsize=axis)
ax2.set_ylabel('Efficiency ($\%$)', fontsize=axis)
ax2.tick_params(labelsize=ticksize)
ax2b.tick_params(labelsize=ticksize)
ax2.set_ylim(0,100)
ax2b.set_ylim(0,25)
ax2b.set_ylabel('Discoveries', fontsize=axis)
#________________________________________________________________________________________________
n3b, bins3b = np.histogram(np.array(filter_dist2(filter_band2(table, 'r-z'))),
bins = 50, range=(0,200))
center3 = (efficiency2(table, 'r-z')[1][:-1] + efficiency2(table, 'r-z')[1][1:]) / 2
width3 = (efficiency2(table, 'r-z')[1][1] - efficiency2(table, 'r-z')[1][0])
ax3.bar(center3, efficiency2(table, 'r-z')[0],
width=width3, align='center', facecolor='yellow', edgecolor='black', label='Efficiency')
ax3b.plot(center3, n3b, 'o-', color='violet', markeredgecolor='black', label='Total Observations')
handles3,labels3 = ax3b.get_legend_handles_labels()
handles3 =handles3[::-1]
labels3 = labels3[::-1]
ax3.legend(handles3, labels3, bbox_to_anchor=(1.05, 1), loc=3, borderaxespad=0., prop={'size':12})
ax3.set_title('r $-$ z', fontsize=subtitle)
ax3.set_xlabel('Distance (Mpc)', fontsize=axis)
ax3.set_ylabel('Efficiency ($\%$)', fontsize=axis)
ax3.tick_params(labelsize=ticksize)
ax3b.tick_params(labelsize=ticksize)
ax3.set_ylim(0,100)
ax3b.set_ylim(0,25)
ax3b.set_ylabel('Discoveries', fontsize=axis)
#________________________________________________________________________________________________
n4b, bins4b = np.histogram(np.array(filter_dist2(filter_band2(table, 'i-z'))),
bins = 50, range=(0,200))
center4 = (efficiency2(table, 'i-z')[1][:-1] + efficiency2(table, 'i-z')[1][1:]) / 2
width4 = (efficiency2(table, 'i-z')[1][1] - efficiency2(table, 'i-z')[1][0])
ax4.bar(center4, efficiency2(table, 'i-z')[0],
width=width4, align='center', facecolor='orange', edgecolor='black', label='Efficiency')
ax4b.plot(center4, n4b, 'o-', color='turquoise', markeredgecolor='black', label='Total Observations')
handles4,labels4 = ax4b.get_legend_handles_labels()
handles4 =handles4[::-1]
labels4 = labels4[::-1]
ax4.legend(handles4, labels4, bbox_to_anchor=(1.05, 1), loc=3, borderaxespad=0., prop={'size':12})
ax4.set_title('i $-$ z', fontsize=subtitle)
ax4.set_xlabel('Distance (Mpc)', fontsize=axis)
ax4.set_ylabel('Efficiency ($\%$)', fontsize=axis)
ax4.tick_params(labelsize=ticksize)
ax4b.tick_params(labelsize=ticksize)
ax4.set_ylim(0,100)
ax4b.set_ylim(0,25)
ax4b.set_ylabel('Discoveries', fontsize=axis)
#________________________________________________________________________________________________
plt.suptitle(name, fontsize=title, y=1.05, horizontalalignment='center')
plt.tight_layout()
return
In [100]:
plot_bar2(APR1215_table2, 'APR4-1215')
In [101]:
plot_bar2(APR1314_table2, 'APR4-1314')
plt.close()
In [102]:
plot_bar2(H41215_table2, 'H4-1215')
plt.close()
In [103]:
plot_bar2(H41314_table2, 'H4-1314')
plt.close()
In [104]:
plot_bar2(Sly135_table2, 'Sly-135')
In [105]:
plot_bar2(APR4Q3a75_table2, 'APR4Q3a75')
plt.close()
In [106]:
plot_bar2(H4Q3a75_table2, 'H4Q3a75')
plt.close()
In [107]:
plot_bar2(MS1Q3a75_table2, 'MS1Q3a75')
plt.close()
In [108]:
plot_bar2(MS1Q7a75_table2, 'MS1Q7a75')
In [ ]: