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 [1]:
#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 [2]:
#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 [250]:
#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, magnitude):
fig=plt.figure()
ax=fig.add_subplot(111)
abs_plotter1(merger)
abs_plot_interp(merger)
LSST_plot(merger)
ax.axis([0,35,rounded_up(correct_max(merger))+5,(rounded_down(correct_min(merger))-5)])
#major_xticks = np.arange(0, 36, 5)
#minor_xticks = np.arange(0, 36, 1)
#ax.set_xticks(major_xticks)
#ax.set_xticks(minor_xticks, minor=True)
#plt.grid(which='minor')
plt.title("Light Curve of " + name, fontsize=18)
plt.xlabel("Time (Days)", fontsize=14)
plt.ylabel('Magnitude ' + magnitude, fontsize=14)
plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
return
#abs_plotter1 Merger -> Image
#Plots the data from the ascii.reads of the merger.
def abs_plotter1(merger):
plot_band(merger, 'u', '.', 'b')
plot_band(merger, 'g', '.', 'g')
plot_band(merger, 'r', '.', 'y')
plot_band(merger, 'i', '.', 'orange')
plot_band(merger, 'z', '.', '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, label=name)
return
#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):
approx='cubic'
u=interp1d(merger['day'][np.isfinite(merger['u'])], merger['u'][np.isfinite(merger['u'])], kind=approx)
g=interp1d(merger['day'][np.isfinite(merger['g'])], merger['g'][np.isfinite(merger['g'])], kind=approx)
r=interp1d(merger['day'][np.isfinite(merger['r'])], merger['r'][np.isfinite(merger['r'])], kind=approx)
i=interp1d(merger['day'][np.isfinite(merger['i'])], merger['i'][np.isfinite(merger['i'])], kind=approx)
z=interp1d(merger['day'][np.isfinite(merger['z'])], merger['z'][np.isfinite(merger['z'])], kind=approx)
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 [251]:
#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", "(Apparent)")
full_plot(APR1314_app, "APR4-1314", "(Apparent)")
full_plot(H41215_app, "H4-1215", "(Apparent)")
full_plot(H41314_app, "H4-1314", "(Apparent)")
full_plot(Sly135_app, "Sly-135", "(Apparent)")
In [249]:
#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', "(Apparent)")
full_plot(H4Q3a75_app, 'H4Q3a75', "(Apparent)")
full_plot(MS1Q3a75_app, 'MS1Q3a75', "(Apparent)")
full_plot(MS1Q7a75_app, 'MS1Q7a75', "(Apparent)")
In [253]:
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=10000
bands=['u','g','r', 'i', 'z']
In [254]:
#rand_plot: Merger String String String -> Image
#Plots a randomly chosen band from the merger.
def rand_plot(merger, name, magnitude):
band=random.choice(bands)
fig=plt.figure()
ax=fig.add_subplot(111)
if band == 'u':
plot_band(merger, 'u', '.', 'b')
plot_interp(merger, 'u', '--', 'b', 0)
elif band == 'g':
plot_band(merger, 'g', '.', 'g')
plot_interp(merger, 'g', '--', 'g', 1)
elif band == 'r':
plot_band(merger, 'r', '.', 'y')
plot_interp(merger, 'r', '--', 'y', 2)
elif band == 'i':
plot_band(merger, 'i', '.', 'orange')
plot_interp(merger, 'i', '--', 'orange', 3)
elif band == 'z':
plot_band(merger, 'z', '.', 'r')
plot_interp(merger, 'z', '--', 'r', 4)
ax.axis([0,35,rounded_up(correct_max(merger))+5,(rounded_down(correct_min(merger))-5)])
plt.title("Light Curve of " + name)
plt.xlabel("Time (Days)")
plt.ylabel('Magnitude ' + magnitude)
plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
return
In [256]:
#Analyze: Merger -> Table
#Randomly runs 10,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):
band_list=[]
dist_list=[]
detect_list=[]
for i in range(trial_num):
band=random.choice(bands)
band_list.append(band)
dist_list.append(rand_dist(merger))
if rand_detect(merger, band) == True:
detect_list.append(1)
else:
detect_list.append(0)
return Table([band_list, dist_list, detect_list], names=('Band', 'Distance', '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):
t_0=random.random()+1.06
if detect(merger, band, t_0) and detect(merger, band, t_0+1) and not detect(merger, band, t_0+6):
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
In [10]:
#NS-NS Mergers
APR1215_table = Analyze(APR1215_app)
APR1314_table = Analyze(APR1314_app)
H41215_table = Analyze(H41215_app)
H41314_table = Analyze(H41314_app)
Sly135_table = Analyze(Sly135_app)
In [11]:
#NS-BH Mergers
APR4Q3a75_table = Analyze(APR4Q3a75_app)
H4Q3a75_table = Analyze(H4Q3a75_app)
MS1Q3a75_table = Analyze(MS1Q3a75_app)
MS1Q7a75_table = Analyze(MS1Q7a75_app)
In [257]:
#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):
f, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, sharex=True, sharey=True)
ax1.hist(np.array(filter_dist(filter_band(filter_hits(table), 'u'))),
bins=49, range=(0,200), facecolor='blue', edgecolor='black', label='u')
ax2.hist(np.array(filter_dist(filter_band(filter_hits(table), 'g'))),
bins=49, range=(0,200), facecolor='green', edgecolor='black', label='g')
ax3.hist(np.array(filter_dist(filter_band(filter_hits(table), 'r'))),
bins=49, range=(0,200), facecolor='yellow', edgecolor='black', label='r')
ax4.hist(np.array(filter_dist(filter_band(filter_hits(table), 'i'))),
bins=49, range=(0,200), facecolor='orange', edgecolor='black', label='i')
ax5.hist(np.array(filter_dist(filter_band(filter_hits(table), 'z'))),
bins=49,facecolor='red', edgecolor='black', label='z')
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
plt.suptitle(name, fontsize=18, y=0.91)
plt.xlabel('Distance (Mpc)', fontsize=14)
ax1.set_ylabel('Detections')
ax2.set_ylabel('Detections')
ax3.set_ylabel('Detections')
ax4.set_ylabel('Detections')
ax5.set_ylabel('Detections')
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]
band_hits = table['Band'][hits]
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 [258]:
plot_hist(APR1215_table, "APR4-1215")
plot_hist(APR1314_table, "APR4-1314")
plot_hist(H41215_table, "H4-1215")
plot_hist(H41314_table, "H4-1314")
plot_hist(Sly135_table, "Sly-135")
In [172]:
plot_hist(APR4Q3a75_table, "APR4Q3a75")
plot_hist(H4Q3a75_table, "H4Q3a75")
plot_hist(MS1Q3a75_table, "MS1Q3a75")
plot_hist(MS1Q7a75_table, "MS1Q7a75")
In [259]:
#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
#Constants:
bar_width=3
distance=np.linspace(0,200,50)
#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):
f, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, sharex=True, sharey=True)
ax1.bar(distance, efficiency(table, 'u'),
width=bar_width, align='center', facecolor='blue', edgecolor='black', label='u')
ax2.bar(distance, efficiency(table, 'g'),
width=bar_width, align='center', facecolor='green', edgecolor='black', label='g')
ax3.bar(distance, efficiency(table, 'r'),
width=bar_width, align='center', facecolor='yellow', edgecolor='black', label='r')
ax4.bar(distance, efficiency(table, 'i'),
width=bar_width, align='center', facecolor='orange', edgecolor='black', label='i')
ax5.bar(distance, efficiency(table, 'z'),
width=bar_width, align='center', facecolor='red', edgecolor='black', label='z')
ax1.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax2.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax3.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax4.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
ax5.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)
plt.suptitle(name, fontsize=18, y=0.91)
plt.xlabel('Distance (Mpc)', fontsize=14)
ax1.set_ylabel('Efficiency(%)')
ax2.set_ylabel('Efficiency(%)')
ax3.set_ylabel('Efficiency(%)')
ax4.set_ylabel('Efficiency(%)')
ax5.set_ylabel('Efficiency(%)')
return
In [260]:
plot_bar(APR1215_table, "APR4-1215")
plot_bar(APR1314_table, "APR4-1314")
plot_bar(H41215_table, "H4-1215")
plot_bar(H41314_table, "H4-1314")
plot_bar(Sly135_table, "Sly-135")
In [261]:
plot_bar(APR4Q3a75_table, "APR4Q3a75")
plot_bar(H4Q3a75_table, "H4Q3a75")
plot_bar(MS1Q3a75_table, "MS1Q3a75")
plot_bar(MS1Q7a75_table, "MS1Q7a75")
In [ ]: