In [22]:
import sqlite3
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
from pyMagnetar import Magnetar
from numba import jit
%matplotlib inline
In [23]:
conn = sqlite3.connect('SLSN.db', detect_types=sqlite3.PARSE_DECLTYPES|sqlite3.PARSE_COLNAMES)
cur = conn.cursor()
In [24]:
# query = """
# DROP TABLE IF EXISTS four_ops
# """
# cur.execute(query)
# conn.commit()
In [25]:
query = """
CREATE TABLE IF NOT EXISTS four_ops (
RUN_ID INT,
T_M REAL,
B REAL,
P REAL,
p0 REAL,
t0 REAL,
M_400_peak REAL,
M_520_peak REAL,
M_400_p30 REAL,
M_520_p30 REAL,
dM_400_p30 REAL,
dM_520_p30 REAL,
C_peak REAL,
C_p30 REAL,
in_A BOOL,
in_B BOOL,
in_C BOOL,
in_D BOOL,
pass BOOL
)
"""
cur.execute(query)
conn.commit()
insert_query = """
INSERT INTO four_ops (
RUN_ID,
T_M, B, P, p0, t0,
M_400_peak, M_520_peak, M_400_p30, M_520_p30,
dM_400_p30, dM_520_p30, C_peak, C_p30,
in_A, in_B, in_C, in_D,
pass
) VALUES (
{},
{}, {}, {}, {}, {},
{}, {}, {}, {},
{}, {}, {}, {},
{}, {}, {}, {},
{}
)
"""
In [26]:
# p = {"A": [-22.62, 0.75, 0.32], "B": [-22.02, 1.14, 0.29], "C": [-0.30, 0.16, 0.14], "D": [-0.22, 0.35, 0.08]}
# xl = {"A": [0, 3.0], "B": [-0.5, 1.5], "C": [0, 3.0], "D": [-0.5, 1.5]}
# yl = {"A": [-23, -20.0], "B": [-23, -20.0], "C": [-0.6, 0.5], "D": [-0.6, 0.5]}
p = {"A": [-21.62, 0.75, 0.62], "B": [-21.02, 1.14, 0.59], "C": [-0.30, 0.16, 0.14], "D": [-0.22, 0.35, 0.08]}
xl = {"A": [0, 3.0], "B": [-0.5, 1.5], "C": [0, 3.0], "D": [-0.5, 1.5]}
yl = {"A": [-23, -19.0], "B": [-23, -19.0], "C": [-1.0, 0.5], "D": [-1.0, 0.5]}
In [27]:
@jit
def line(x, p):
return p[0] + p[1] * x
@jit
def check(x, y, panel):
if not (x >= xl[panel][0] and x < xl[panel][1]):
return False
if not (y >= yl[panel][0] and y < yl[panel][1]):
return False
if np.abs(y - line(x, p[panel])) > 3.0*p[panel][2]:
return False
return True
In [28]:
m = Magnetar(b"filters")
t = np.arange(1, 200, 1.0)
zp400 = 20.39
zp520 = 20.96
In [29]:
def sim(T, B, P, t0, run):
# Set up the magnetar for given params
m.setup(T, B, P, t0, 0.0)
# calculate LC for 400nm band
flux400 = m.flux(t, b"B400")
mag400 = -2.5 * np.log10(flux400) - zp400
p0 = np.argmax(flux400)
if p0 > 150:
return False
# calculate LC for 520nm band
flux520 = m.flux(t, b"B520")
mag520 = -2.5 * np.log10(flux520) - zp520
if np.isnan(flux400[p0]) or np.isnan(flux520[p0]) or np.isnan(flux400[p0+30]) or np.isnan(flux520[p0+30]):
return False
# 4OPS params
M_400_peak = mag400[p0]
M_520_peak = mag520[p0]
M_400_p30 = mag400[p0 + 30]
M_520_p30 = mag520[p0 + 30]
dM_400_p30 = M_400_p30 - M_400_peak
dM_520_p30 = M_520_p30 - M_520_peak
C_peak = M_400_peak - M_520_peak
C_p30 = M_400_p30 - M_520_p30
# 4OPS checks
check_A = check(dM_400_p30, M_400_peak, "A")
check_B = check(C_p30, M_400_peak, "B")
check_C = check(dM_400_p30, C_peak, "C")
check_D = check(C_p30, C_peak, "D")
check_all = int((check_A & check_D) | (check_B & check_C))
insert = insert_query.format(run,
T, B, P, p0, t0,
M_400_peak, M_520_peak, M_400_p30, M_520_p30,
dM_400_p30, dM_520_p30, C_peak, C_p30,
int(check_A), int(check_B), int(check_C), int(check_D),
check_all
)
cur.execute(insert)
conn.commit()
return check_all
In [30]:
# for T in np.arange(10.0, 150.0, 5.0):
# print(T, '-', datetime.now())
# for B in np.arange(0.1, 20, 0.1):
# for P in np.arange(0.1, 10, 0.1):
# for t0 in np.arange(0.0, 20.0, 5.0):
# sim(T, B, P, t0, 7)
In [31]:
# query = """
# SELECT * FROM four_ops WHERE run_id=4 AND pass=1
# """
# cur.execute(query)
# df = pd.DataFrame(cur.fetchall(), columns=np.array(cur.description)[:,0])
In [32]:
# for index, row in df.iterrows():
# sim(row['T_M'], row['B'], row['P'], 15.0, 7)
In [33]:
count = 0
bPrint = True
while count < 100:
T = np.random.random() * 150.0 + 10.0
B = np.random.random() * 20.0 + 0.01
P = np.random.random() * 10.0 + 0.01
if count % 1000 == 0 and bPrint == True:
print(count, ': ', datetime.now())
bPrint = False
if sim(T, B, P, 0.0, 99):
bPrint = True
count += 1
In [21]:
conn.close()
In [ ]: