<jupyter_start><jupyter_text>Adding SNe to CoCo___This notebook will go through adding a SN to the CoCo template framework.<jupyter_code>%matplotlib inline
# %matplotlib notebook

try:
    reload  # Python 2.7
except NameError:
    try:
        from importlib import reload  # Python 3.4+
    except ImportError:
        from imp import reload  # Python 3.0 - 3.3
        
from matplotlib import pyplot as plt
from matplotlib import rc
rc('text', usetex=True)

import os
import numpy as np
from astropy.table import Table, Column
from collections import OrderedDict

import pycoco as pcc<jupyter_output><empty_output><jupyter_text>Download SN---The example that I will go through is probably the simplest one - from a JSON that I have downloaded from http://sne.space contaning the information for SN2010as.<jupyter_code>import json
snname = "SN2010as"
path = "/Users/berto/data/CoreCollapse/SNe/SN2010as.json"

sndict = json.load(open(path, 'r'))[snname]
len(sndict["spectra"])
S1 = sndict["spectra"][0]
S1.keys()
w = np.array(S1["data"], dtype=np.float64).T[0]
f = np.array(S1["data"], dtype=np.float64).T[1]
f
from astropy.table import Table, Column, Row
import astropy.units as u

wavelength = Column(w, name = ("wavelength"), unit=u.Angstrom)
flux = Column(f, name = ("flux"), unit=u.cgs.erg / u.si.cm ** 2 / u.si.s / u.angstrom)
Table([wavelength, flux])
spec_dict_list = sndict["spectra"]
for s in spec_dict_list:
    print(s["time"], s["filename"])
print(len(spec_dict_list))
to_combine = ["SN10as_100325_VLT_XS_", 
              "SN10as_100326_VLT_XS_", 
              "SN10as_100328_VLT_XS_", 
              "SN10as_100329_VLT_XS_"]

to_combine_phases = [55280.0,
                     55281.0,
                     55283.0,
                     55284.0]

spec_filenames = [spec["filename"] for spec in sndict["spectra"]]
# spec_filenames = [spec["filename"][:-7] for spec in sndict["spectra"]]
mask = np.ones(len(spec_dict_list), dtype=bool)

combined = np.array([])

for prefix in to_combine:
    matches = [f for f in spec_filenames if prefix in f]
    print(matches)
    
    spec_dir = "/Users/berto/data/CoreCollapse/data/spectra/SN2010as/"
    
    tables = []
    out_dict = OrderedDict()
    for file in matches:
        w = np.where(np.array(spec_filenames) == file)[0]
        
        wavl = np.array(spec_dict_list[w]["data"], dtype=np.float64).T[0]
        flux = np.array(spec_dict_list[w]["data"], dtype=np.float64).T[1]
        
        wavelength = Column(wavl, name = ("wavelength"), unit=u.Angstrom)
        flux = Column(flux, name = ("flux"), unit=u.cgs.erg / u.si.cm ** 2 / u.si.s / u.angstrom)
     
        tables.append(Table([wavelength, flux]))
        
        for i in w:
            print(i, spec_dict_list[i]["filename"])
            mask[i] = False
#             spec_dict_list.pop(i)
            
    tables = np.array(tables)[np.argsort([i["wavelength"][0] for i in tables])]
    
    S1 = pcc.SpectrumClass()
    S2 = pcc.SpectrumClass()
    S3 = pcc.SpectrumClass()
    
    S1.load_table(tables[0], path="")
    S2.load_table(tables[1], path="")
    S3.load_table(tables[2], path="")
    
    sr=pcc.functions.combine_spectra(S2, S3, 9975, 10150, scale=True)
    spec = pcc.functions.combine_spectra(S1,sr, 5502.5, 5525)
    
    out_spec = pcc.SpectrumClass()
    out_spec.load_table(spec.data, path="", trim_wavelength=True, wmin=0, wmax=15000)
    
    filename = prefix+"combined.dat"
    
    out_dict["filename"] = os.path.join(spec_dir, filename)
    out_dict["time"] = np.float64(spec_dict_list[w]["time"])
    out_dict["redshift"] = spec_dict_list[w]["redshift"]
    out_dict["snname"] = snname
    
    combined = np.append(combined, out_dict)
    
#     spec.plot()
#     spec.save(filename, path=spec_dir, format="ascii.commented_header", verbose=True)
    out_spec.plot()
    out_spec.save(filename, path=spec_dir, format="ascii.commented_header", verbose=True)
out_dict
uncombined=np.array([])

for spec_file in np.array(spec_filenames)[mask]:
    print(spec_file)
    
    w = np.where(np.array(spec_filenames) == spec_file)[0]
    print(w)
    
    wavl = np.array(spec_dict_list[w]["data"], dtype=np.float64).T[0]
    flux = np.array(spec_dict_list[w]["data"], dtype=np.float64).T[1]

    wavelength = Column(wavl, name = ("wavelength"), unit=u.Angstrom)
    flux = Column(flux, name = ("flux"), unit=u.cgs.erg / u.si.cm ** 2 / u.si.s / u.angstrom)

    spec_table = Table([wavelength, flux])

    s = pcc.SpectrumClass()

    s.load_table(spec_table, path="")
    print(os.path.join(spec_dir, spec_file))
    s.save(spec_file, path=spec_dir, )
    out_dict = OrderedDict()
    out_dict["filename"] = os.path.join(spec_dir, spec_file)
    out_dict["time"] = np.float64(spec_dict_list[w]["time"])
    out_dict["redshift"] = spec_dict_list[w]["redshift"]
    out_dict["snname"] = snname
    print(spec_dict_list[w]["time"])
    uncombined = np.append(uncombined, out_dict)
    
    s.plot()
out_array = np.append(combined, uncombined)
out_array
# [i["filename"] for i in out_array]
# [i["redshift"] for i in out_array]
# [i["snname"] for i in out_array]
[i["time"] for i in out_array]<jupyter_output><empty_output><jupyter_text>We need a (rough) estimate of the time of maximum light (only to ensure the spectra adhere to our naming convention), so from <NAME> al. 2014 (http://iopscience.iop.org/article/10.1088/0004-637X/792/1/7/pdf) we have:<jupyter_code>mjdmax = 55286.4
# keys = out_array[0].keys()

keys = ["filename",
        "snname",
        "time",
        "redshift"]
columns = []
for key in keys:
    columns.append([i[key] for i in out_array])
keys[2] = "MJD"
list_table = Table(columns, names = keys)
list_table.add_index("MJD")
list_table.sort()
list_table
list_out_path = os.path.join(pcc.defaults._default_list_dir_path, snname+".list")
list_table.write(list_out_path, format="ascii.no_header")
info = pcc.classes.InfoClass()
info.load()
info._data
small_table = info._data[("snname", "z_obs", "mu", "RA", "Dec", "Type", "z_distmod")]
small_table
sndict.keys()
sndict["lumdist"]
from astropy.cosmology import LambdaCDM
from astropy.coordinates import Distance

mu = 32.16 ## from NED -- https://ned.ipac.caltech.edu/ui/?q=objsearch&omegav=0.73&obj_sort=RA%20or%20Longitude&omegam=0.27&extend=no&hconst=73&out_csys=Equatorial&of=pre_text&img_stamp=YES&list_limit=5&out_equinox=J2000.0&corr_z=1&objname=NGC%206000

d = Distance(distmod=mu)

z_distmod = d.z
(snname, sndict["redshift"][0]["value"], mu, sndict["ra"][0]["value"], sndict["dec"][0]["value"], sndict["claimedtype"][0]["value"].replace("/", ""), z_distmod)
small_table.add_row((snname, sndict["redshift"][0]["value"], mu, sndict["ra"][0]["value"], sndict["dec"][0]["value"], sndict["claimedtype"][0]["value"].replace("/", ""), z_distmod))
small_table
outfile = pcc.defaults._default_info_path.replace(".dat", "_new.dat")
small_table.write(outfile, format="ascii.commented_header")

SN = pcc.classes.SNClass(snname)
SN.load_phot()
SN.plot_lc()
list_file_path = "/Users/berto/Code/CoCo/lists/"+snname+".list"
# list_file_path = "/Users/berto/data/CoreCollapse/SNe/SN2010as.list"
SN.load_list(path=list_file_path)
SN.load_spec()
SN.plot_spec()
SN.spec
fit_path = "/Users/berto/Code/CoCo/recon/"+snname+".dat"
SN.get_lcfit(fit_path)

SN.plot_lc()
# pcc.coco.run_specfit(list_file_path)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16029391753980132, 'fraction_numerical': 0.030072118655599403, 'mean_word_length': 4.104166666666667, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '433871', 'license': 'MIT', 'n_tokens_mistral': 2943, 'n_tokens_neox': 2777, 'n_words': 506, 'path': 'notebooks/Tutorials/pycoco_adding_SNe.ipynb', 'repo_name': 'RobFirth/pycoco', 'size': '7355'}