<commit_before>import tables as tb
import numpy as np
from astropy import units as u
from mag2flux import mag2flux
from fbol import integrate_fqbol as fqbol_trapezoidal
from fbol import ir_correction, uv_correction_linear, uv_correction_blackbody
from fit_blackbody import bb_fit_parameters
from fit_blackbody import bb_flux_nounits
from luminosity import calc_Lbol
from specutils import extinction


class SN(object):
    """A supernova is the explosion that ends the life of a star

    The SN needs to be conatained within the HDF5 database before it is used
    by SNoBoL. Once there, simply create a supernova by calling the constructor
    with the name of the SN as a string of the form "sn[YEAR][Letter(s)]"

    For example:
    sn1987A = SN('sn1987a')
    sn1999em = SN('sn1999em')

    Attributes
    ----------
    name : Name of the supernova, "sn" followed by the year of first observation
           along with a letter designating the order of observation in that
           year. "sn1987a" was the first SN observed in 1987. "sn2000cb" was the
           eightieth SN observed in 2000.
    """

    def __init__(self, name):
        """Initializes the SN with supplied value for [name]"""
        self.name = name
        self.min_num_obs = 4

        self.read_hdf5()

    def read_hdf5(self):
        """Reads the hdf5 file and returns data on supernova matching [name]"""
        hdf5_filename = './hdf5/sn_data.h5'
        h5file = tb.open_file(hdf5_filename, 'r')
        
        self.filter_table = h5file.root.filters
        
        self.sn_node = h5file.get_node('/sn', self.name)
        
        self.phot_table = self.sn_node.phot
        self.parameter_table = self.sn_node.parameters

    def lbol_direct_bh09(self):
        """Calculate the bolometric lightcurve using the direct integration
           method published in Bersten & Hamuy 2009 (2009ApJ...701..200B)"""
        self.convert_magnitudes_to_fluxes()
        self.deredden_fluxes()
        self.get_lbol_epochs()
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()
        
        self.lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for jd in self.lbol_epochs:
            names = np.array([x['name'] for x in self.converted_obs 
                              if x['jd'] == jd and x['name'] != 'z'])
            wavelengths = np.array([x['wavelength'] for x in self.converted_obs
                                    if x['jd'] == jd and x['name'] != 'z'])
            fluxes = np.array([x['flux'] for x in self.converted_obs
                               if x['jd'] == jd and x['name'] != 'z'])
            flux_errs = np.array([x['uncertainty'] for x in self.converted_obs
                                  if x['jd'] == jd and x['name'] != 'z'])

            sort_indices = np.argsort(wavelengths)
            wavelengths = wavelengths[sort_indices]
            fluxes = fluxes[sort_indices]
            flux_errs = flux_errs[sort_indices]

            fqbol, fqbol_err = fqbol_trapezoidal(wavelengths, fluxes, flux_errs)
            temperature, angular_radius, perr = bb_fit_parameters(wavelengths,
                                                                   fluxes,
                                                                   flux_errs)

            temperature_err = perr[0]
            angular_radius_err = perr[1]

            shortest_wl = np.amin(wavelengths)
            shortest_flux = np.amin(fluxes)
            shortest_flux_err = np.amin(flux_errs)
            longest_wl = np.amax(wavelengths)

            ir_corr, ir_corr_err = ir_correction(temperature,
                                                 temperature_err,
                                                 angular_radius,
                                                 angular_radius_err,
                                                 longest_wl)
            if 'U' in names:
                idx = np.nonzero(names == 'U')[0][0]
                U_flux = fluxes[idx]
                U_wl = wavelengths[idx]
                if U_flux < bb_flux_nounits(U_wl,
                                            temperature, 
                                            angular_radius):
                    uv_corr, uv_corr_err = uv_correction_linear(shortest_wl, 
                                                                shortest_flux, 
                                                                shortest_flux_err) 
                else:
                    uv_corr, uv_corr_err = uv_correction_blackbody(temperature,
                                                             temperature_err,
                                                             angular_radius,
                                                             angular_radius_err,
                                                             shortest_wl)
            else:
                uv_corr, uv_corr_err = uv_correction_blackbody(temperature,
                                                               temperature_err,
                                                               angular_radius,
                                                               angular_radius_err,
                                                               shortest_wl)

            fbol = fqbol + ir_corr + uv_corr
            fbol_err = np.sqrt(np.sum(x*x for x in [fqbol_err, ir_corr_err, uv_corr_err]))
            lum = fbol * 4.0 * np.pi * self.distance_cm**2.0
            lum_err = np.sqrt((4.0 * np.pi * self.distance_cm**2 * fbol_err)**2
                              +(8.0*np.pi * fbol * self.distance_cm * self.distance_cm_err)**2)
            phase = jd - self.parameter_table.cols.explosion_JD[0]
            phase_err = self.parameter_table.cols.explosion_JD_err[0]
            self.lc = np.append(self.lc, [[jd, phase, phase_err, lum, lum_err, temperature, angular_radius]], axis=0)

        self.lc = np.delete(self.lc, (0), axis=0)

        self.write_lbol_plaintext(self.lc, 'direct')

    def lqbol(self):
        """Calculate the quasi-bolometric lightcurve using direct integration
           with trapezoidal integration of the fluxes"""
        self.convert_magnitudes_to_fluxes()
        self.deredden_fluxes()
        self.get_lbol_epochs()
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()
        
        self.qbol_lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for jd in self.lbol_epochs:
            wavelengths = np.array([x['wavelength'] for x in self.converted_obs
                                    if x['jd'] == jd and x['name'] != 'z'])
            fluxes = np.array([x['flux'] for x in self.converted_obs
                               if x['jd'] == jd and x['name'] != 'z'])
            flux_errs = np.array([x['uncertainty'] for x in self.converted_obs
                                  if x['jd'] == jd and x['name'] != 'z'])

            sort_indices = np.argsort(wavelengths)
            wavelengths = wavelengths[sort_indices]
            fluxes = fluxes[sort_indices]
            flux_errs = flux_errs[sort_indices]

            fqbol, fqbol_err = fqbol_trapezoidal(wavelengths, fluxes, flux_errs)

            lqbol = fqbol * 4.0 * np.pi * self.distance_cm**2.0
            lqbol_err = np.sqrt((4.0 * np.pi * self.distance_cm**2 * fqbol_err)**2
                              +(8.0*np.pi * fqbol * self.distance_cm * self.distance_cm_err)**2)
            phase = jd - self.parameter_table.cols.explosion_JD[0]
            phase_err = self.parameter_table.cols.explosion_JD_err[0]
            self.qbol_lc = np.append(self.qbol_lc, [[jd, phase, phase_err, lqbol, lqbol_err]], axis=0)
        self.qbol_lc = np.delete(self.qbol_lc, (0), axis=0)
        self.write_lbol_plaintext(self.qbol_lc, 'qbol')

    def lbol_bc_bh09(self, filter1, filter2):
        """Calculate the bolometric lightcurve using the bolometric corrections
           found in Bersten & Hamuy 2009 (2009ApJ...701..200B). These require 
           specifying a color, taken to be filter1 - filter2"""
        self.get_magnitudes()
        self.deredden_UBVRI_magnitudes()
        self.get_bc_epochs(filter1, filter2)
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()

        self.bc_lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for i in range(len(self.bc_epochs)):
            jd = self.bc_epochs[i]
            color = self.get_bc_color(jd, filter1, filter2)
            color_err = self.get_bc_color_uncertainty(jd, filter1, filter2) 
            v_mag = np.array([x['magnitude'] for x in self.photometry 
                               if x['jd'] == jd and x['name'] == 'V'])
            v_mag_err = np.array([x['uncertainty'] for x in self.photometry <commit_msg>Fix IR-only nights<commit_after>import tables as tb
import numpy as np
from astropy import units as u
from mag2flux import mag2flux
from fbol import integrate_fqbol as fqbol_trapezoidal
from fbol import ir_correction, uv_correction_linear, uv_correction_blackbody
from fit_blackbody import bb_fit_parameters
from fit_blackbody import bb_flux_nounits
from luminosity import calc_Lbol
from specutils import extinction


class SN(object):
    """A supernova is the explosion that ends the life of a star

    The SN needs to be conatained within the HDF5 database before it is used
    by SNoBoL. Once there, simply create a supernova by calling the constructor
    with the name of the SN as a string of the form "sn[YEAR][Letter(s)]"

    For example:
    sn1987A = SN('sn1987a')
    sn1999em = SN('sn1999em')

    Attributes
    ----------
    name : Name of the supernova, "sn" followed by the year of first observation
           along with a letter designating the order of observation in that
           year. "sn1987a" was the first SN observed in 1987. "sn2000cb" was the
           eightieth SN observed in 2000.
    """

    def __init__(self, name):
        """Initializes the SN with supplied value for [name]"""
        self.name = name
        self.min_num_obs = 4

        self.read_hdf5()

    def read_hdf5(self):
        """Reads the hdf5 file and returns data on supernova matching [name]"""
        hdf5_filename = './hdf5/sn_data.h5'
        h5file = tb.open_file(hdf5_filename, 'r')
        
        self.filter_table = h5file.root.filters
        
        self.sn_node = h5file.get_node('/sn', self.name)
        
        self.phot_table = self.sn_node.phot
        self.parameter_table = self.sn_node.parameters

    def lbol_direct_bh09(self):
        """Calculate the bolometric lightcurve using the direct integration
           method published in Bersten & Hamuy 2009 (2009ApJ...701..200B)"""
        self.convert_magnitudes_to_fluxes()
        self.deredden_fluxes()
        self.get_lbol_epochs()
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()
        
        self.lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for jd in self.lbol_epochs:
            names = np.array([x['name'] for x in self.converted_obs 
                              if x['jd'] == jd and x['name'] != 'z'])
            wavelengths = np.array([x['wavelength'] for x in self.converted_obs
                                    if x['jd'] == jd and x['name'] != 'z'])
            fluxes = np.array([x['flux'] for x in self.converted_obs
                               if x['jd'] == jd and x['name'] != 'z'])
            flux_errs = np.array([x['uncertainty'] for x in self.converted_obs
                                  if x['jd'] == jd and x['name'] != 'z'])

            sort_indices = np.argsort(wavelengths)
            wavelengths = wavelengths[sort_indices]
            fluxes = fluxes[sort_indices]
            flux_errs = flux_errs[sort_indices]

            fqbol, fqbol_err = fqbol_trapezoidal(wavelengths, fluxes, flux_errs)
            temperature, angular_radius, perr = bb_fit_parameters(wavelengths,
                                                                   fluxes,
                                                                   flux_errs)

            temperature_err = perr[0]
            angular_radius_err = perr[1]

            shortest_wl = np.amin(wavelengths)
            shortest_flux = np.amin(fluxes)
            shortest_flux_err = np.amin(flux_errs)
            longest_wl = np.amax(wavelengths)

            ir_corr, ir_corr_err = ir_correction(temperature,
                                                 temperature_err,
                                                 angular_radius,
                                                 angular_radius_err,
                                                 longest_wl)
            if 'U' in names:
                idx = np.nonzero(names == 'U')[0][0]
                U_flux = fluxes[idx]
                U_wl = wavelengths[idx]
                if U_flux < bb_flux_nounits(U_wl,
                                            temperature, 
                                            angular_radius):
                    uv_corr, uv_corr_err = uv_correction_linear(shortest_wl, 
                                                                shortest_flux, 
                                                                shortest_flux_err) 
                else:
                    uv_corr, uv_corr_err = uv_correction_blackbody(temperature,
                                                             temperature_err,
                                                             angular_radius,
                                                             angular_radius_err,
                                                             shortest_wl)
            else:
                uv_corr, uv_corr_err = uv_correction_blackbody(temperature,
                                                               temperature_err,
                                                               angular_radius,
                                                               angular_radius_err,
                                                               shortest_wl)

            fbol = fqbol + ir_corr + uv_corr
            fbol_err = np.sqrt(np.sum(x*x for x in [fqbol_err, ir_corr_err, uv_corr_err]))
            lum = fbol * 4.0 * np.pi * self.distance_cm**2.0
            lum_err = np.sqrt((4.0 * np.pi * self.distance_cm**2 * fbol_err)**2
                              +(8.0*np.pi * fbol * self.distance_cm * self.distance_cm_err)**2)
            phase = jd - self.parameter_table.cols.explosion_JD[0]
            phase_err = self.parameter_table.cols.explosion_JD_err[0]
            self.lc = np.append(self.lc, [[jd, phase, phase_err, lum, lum_err, temperature, angular_radius]], axis=0)

        self.lc = np.delete(self.lc, (0), axis=0)

        self.write_lbol_plaintext(self.lc, 'direct')

    def lqbol(self):
        """Calculate the quasi-bolometric lightcurve using direct integration
           with trapezoidal integration of the fluxes"""
        self.convert_magnitudes_to_fluxes()
        self.deredden_fluxes()
        self.get_lbol_epochs()
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()
        
        self.qbol_lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for jd in self.lbol_epochs:
            wavelengths = np.array([x['wavelength'] for x in self.converted_obs
                                    if x['jd'] == jd])
            fluxes = np.array([x['flux'] for x in self.converted_obs
                               if x['jd'] == jd])
            flux_errs = np.array([x['uncertainty'] for x in self.converted_obs
                                  if x['jd'] == jd])
            names = np.array([x['name'] for x in self.converted_obs
                                  if x['jd'] == jd])

            sort_indices = np.argsort(wavelengths)
            wavelengths = wavelengths[sort_indices]
            fluxes = fluxes[sort_indices]
            flux_errs = flux_errs[sort_indices]
            names = names[sort_indices]

            fqbol, fqbol_err = fqbol_trapezoidal(wavelengths, fluxes, flux_errs)

            lqbol = fqbol * 4.0 * np.pi * self.distance_cm**2.0
            lqbol_err = np.sqrt((4.0 * np.pi * self.distance_cm**2 * fqbol_err)**2
                              +(8.0*np.pi * fqbol * self.distance_cm * self.distance_cm_err)**2)
            phase = jd - self.parameter_table.cols.explosion_JD[0]
            phase_err = self.parameter_table.cols.explosion_JD_err[0]
            # Quick and dirty fix for IR-only nights (don't want those in qbol calc)
            if min(wavelengths) < 10000.0:
                self.qbol_lc = np.append(self.qbol_lc, [[jd, phase, phase_err, lqbol, lqbol_err]], axis=0)
        self.qbol_lc = np.delete(self.qbol_lc, (0), axis=0)
        self.write_lbol_plaintext(self.qbol_lc, 'qbol')

    def lbol_bc_bh09(self, filter1, filter2):
        """Calculate the bolometric lightcurve using the bolometric corrections
           found in Bersten & Hamuy 2009 (2009ApJ...701..200B). These require 
           specifying a color, taken to be filter1 - filter2"""
        self.get_magnitudes()
        self.deredden_UBVRI_magnitudes()
        self.get_bc_epochs(filter1, filter2)
        self.distance_cm, self.distance_cm_err = self.get_distance_cm()

        self.bc_lc = np.array([[0.0, 0.0, 0.0, 0.0, 0.0]])
        
        for i in range(len(self.bc_epochs)):
            jd = self.bc_epochs[i]
            color = self.get_bc_color(jd, filter1, filter2)
            color_err = self.get_bc_color_uncertainty(jd, filter1, filter2) 
            v_mag = np.array([x['magnitude'] for x in self.photometry 
                               if x['jd'] == jd and x['name'] == 'V'])
            v_mag_err = np.array([x['uncertainty'] for x in self.photometry 
{'fraction_non_alphanumeric': 0.09380451297057003, 'fraction_numerical': 0.018119408024309266, 'mean_word_length': 1.5483223401204473, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 89, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2022721', 'n_tokens_mistral': 5451, 'n_tokens_neox': 5028, 'n_words': 1182}