def load_filters(FILTER_PREFIX='tophat_'):
    '''
    Load UBVRI tophat filters defined in Pereira (2013) into the sncosmo
    registry.  Also returns a dictionary of zero-point values in Vega system
    flux for each filter.  Flux is given in [photons/s/cm^2].

    Example usage:

        ZP_CACHE = loader.load_filters()
        U_band_zp = ZP_CACHE['U']
    '''
    # dictionary for zero point fluxes of filters
    ZP_CACHE = {'prefix': FILTER_PREFIX}

    for f in 'UBVRI':
        filter_name = FILTER_PREFIX+f
        file_name = dirname+'/data/filters/'+filter_name+'.dat'
        
        try:
            snc.get_bandpass(filter_name)
        except:
            data = np.genfromtxt(file_name)
            bandpass = snc.Bandpass( data[:,0], data[:,1] )
            snc.registry.register(bandpass, filter_name)

        zpsys = snc.get_magsystem('ab')
        zp_phot = zpsys.zpbandflux(filter_name)

        ZP_CACHE[f] = zp_phot

    return ZP_CACHE

def test_register():
    disp = np.array([4000., 4200., 4400., 4600., 4800., 5000.])
    trans = np.array([0., 1., 1., 1., 1., 0.])

    # create a band, register it, make sure we can get it back.
    band = sncosmo.Bandpass(disp, trans, name='tophatg')
    sncosmo.register(band)
    assert sncosmo.get_bandpass('tophatg') is band

    # test deprecated path to registry
    band = sncosmo.Bandpass(disp, trans, name='tophatg2')
    sncosmo.registry.register(band)
    assert sncosmo.get_bandpass('tophatg2') is band

def plot_bandpass_set(setname):
    """Plot the given set of bandpasses."""

    rc("font", family="serif")

    bandpass_meta = registry.get_loaders_metadata(Bandpass)

    fig = plt.figure(figsize=(9, 3))
    ax = plt.axes()

    nbands = 0
    for m in bandpass_meta:
        if m['filterset'] != setname:
            continue
        b = get_bandpass(m['name'])
        ax.plot(b.wave, b.trans, label=m['name'])
        nbands += 1

    ax.set_xlabel("Wavelength ($\\AA$)")
    ax.set_ylabel("Transmission")

    ncol = 1 + (nbands-1) // 9  # 9 labels per column
    ax.legend(loc='upper right', frameon=False, fontsize='small',
              ncol=ncol)

    # Looks like each legend column takes up about 0.125 of the figure.
    # Make room for the legend.
    xmin, xmax = ax.get_xlim()
    xmax += ncol * 0.125 * (xmax - xmin)
    ax.set_xlim(xmin, xmax)
    plt.tight_layout()
    plt.show()

    def filtcheck(self, bandpass, z, frac=0.75, survey="Euclid", f_index=-1):
        """
        check if the redshifted effective wavelength is redder than the effective
        wavelength of the reddest filter (yes, its a complicated sentence)
        Input is a bandpass (as a string) and redshift
        
        """
        bp_rest = sncosmo.get_bandpass(bandpass)
        if survey == "Euclid":
            effwave = self.effwave_arr
            filtarr = self.filtarr
            
        elif survey == "LSST":
            effwave = self.lsst_effwave_arr
            filtarr = self.lsst_filtarr
            
        if bp_rest.wave_eff*(1+z)  > effwave[f_index]:
            filtfun=filtarr[effwave==effwave[f_index]][0]

            #condition: check what fraction of the redshifted filter is 
            cond = bp_rest.wave*(1+z) < max(filtfun.wave[filtfun.trans > 1e-4])
            
            simp_prod = simps(bp_rest.trans, bp_rest.wave)
            if len(bp_rest.wave[cond])>10:
                simp_prod_cond = simps(bp_rest.trans[cond],bp_rest.wave[cond])
            else:
                simp_prod_cond=0
            if simp_prod_cond/simp_prod > frac:
                return 1
            else:
                print "rest-frame filter is too red for observation"
                return 0
        else:        
            return 1

def main():
        filters_bucket, zp_bucket = l.generate_buckets(3300, 9700, N_BUCKETS, inverse_microns=True)

        
        filter_eff_waves = np.array([snc.get_bandpass(zp_bucket['prefix']+f).wave_eff
                                     for f in filters_bucket]
                                    )
        
        sn12cu_excess, phases = load_12cu_excess(filters_bucket, zp_bucket)
        
        
        fig = plt.figure()
        
        for i, phase in enumerate(phases):
                print "Plotting phase {} ...".format(phase)
                
                ax = plt.subplot(2,6,i+1)
                
                plot_contour(i, phase, redden_fm, sn12cu_excess[i], filter_eff_waves,
                             EBV_GUESS, EBV_PAD, RV_GUESS, RV_PAD, STEPS, ax)
                             
        
        fig.subplots_adjust(left=0.04, bottom=0.08, right=0.95, top=0.92, hspace=.06, wspace=.1)
        fig.suptitle('SN2012CU: $E(B-V)$ vs. $R_V$ Contour Plot per Phase', fontsize=TITLE_FONTSIZE)
        plt.show()

def getlsstbandpassobjs(loadsncosmo=True, loadcatsim=True, plot=True):
    """
    General utility to return a list of the baseline LSST bandpasses loaded
    as catsim bandpass objects, and register them as SNCosmo bandpasses
    accessible through strings like 'LSSTu'.
    args:
        loadsncosmo: Bool, optional, defaults to True
            variable to decide whether to register the LSST bandpasses as
            SNCosmo registered bandpass objects accessible through strings
            like 'LSSTu'
        loadcatsim : Bool, optional, defaults to True
            variable to decide whether to set up catsim bandpass objects
            are return the list of u,g,r,i,z,y bandpasses
    returns:
        if loadcatsim is true, list of catsim bandpass objects corresponding
        to LSST baseline u, g, r, i, z, y filters.
        if loadcatsim is False, return is None

    Examples:

    """
    bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
    banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
    lsstbands = []
    lsstbp = {}

    for band in bandPassList:
        # setup sncosmo bandpasses
        bandfname = banddir + "/total_" + band + '.dat'
        if loadsncosmo:
            # register the LSST bands to the SNCosmo registry
            # Not needed for LSST, but useful to compare independent codes
            # Usually the next two lines can be merged,
            # but there is an astropy bug currently which affects only OSX.
            numpyband = np.loadtxt(bandfname)
            sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
                                           trans=numpyband[:, 1],
                                           wave_unit=Unit('nm'),
                                           name='LSST' + band)
            sncosmo.registry.register(sncosmoband, force=True)
        if loadcatsim:
            # Now load LSST bandpasses for catsim
            lsstbp[band] = Bandpass()
            lsstbp[band].readThroughput(bandfname, wavelen_step=wavelenstep)
            lsstbands.append(lsstbp[band])
    fifilterfigs, filterax = plt.subplots()
    if plot:
        for band in bandPassList:
            b = sncosmo.get_bandpass('LSST' + band)
            filterax.plot(b.wave, b.trans, '-k', lw=2.0)
            filterax.set_xlabel(r'$\lambda$, ($\AA$)')
            filterax.set_ylabel(r'transmission')
            filterax.set_label(r'LSST Filter Transmission Functions')
        plt.show()

    if loadcatsim:
        return lsstbands
    else:
        return None

    def test_bandflux(self):
        self.source.set(amplitude=1.0)
        f = self.source.bandflux("bessellb", 0.0)

        # Correct answer
        b = sncosmo.get_bandpass("bessellb")
        ans = np.sum(b.trans * b.wave * b.dwave) / sncosmo.models.HC_ERG_AA
        assert_approx_equal(ans, f)

 def set_mag_keys(self,key_mag,key_magerr):
     #
     # add lbda def
     #
     super(SDSSCatalogue,self).set_mag_keys(key_mag,key_magerr)
     if key_mag is not None:
         bandpass = get_bandpass("sdss%s"%key_mag[0])
         self.lbda = bandpass.wave_eff

def test_register():
    disp = np.array([4000., 4200., 4400., 4600., 4800., 5000.])
    trans = np.array([0., 1., 1., 1., 1., 0.])
    band = sncosmo.Bandpass(disp, trans, name='tophatg')
    sncosmo.registry.register(band)

    band2 = sncosmo.get_bandpass('tophatg')

    # make sure we can get back the band we put it.
    assert band2 is band

def plotexcess(phases, name, loader, EBV, RV, filters, zp, ax, AV=0.0, P=0.0, plotpl=False):
    
    print "Plotting",name,"..."
    
    ref = loader(phases, filters, zp)
    prefix = zp['prefix']
    filter_eff_waves = [snc.get_bandpass(prefix+f).wave_eff for f in filters]

    # get 11fe synthetic photometry at BMAX, get ref sn color excesses at BMAX
    sn11fe = l.interpolate_spectra(phases, l.get_11fe())

    if type(sn11fe) == type(()):  # convert from tuple to list if just one phase
        sn11fe = [sn11fe]
    
    for i, phase, sn11fe_phase in izip(xrange(len(phases)), phases, sn11fe):
        
        # calculate sn11fe band magnitudes
        sn11fe_mags = {f : -2.5*np.log10(sn11fe_phase[1].bandflux(prefix+f)/zp[f])
                       for f in filters}
        
        # calculate V-X colors for sn11fe
        sn11fe_colors = [sn11fe_mags['V']-sn11fe_mags[f] for f in filters]

        # make list of colors of reference supernova for given phase i
        ref_colors = [ref[i][f] for f in filters]

        # get colors excess of reference supernova compared for sn11fe
        phase_excesses = np.array(ref_colors)-np.array(sn11fe_colors)

        # convert effective wavelengths to inverse microns then plot
        eff_waves_inv = (10000./np.array(filter_eff_waves))
        mfc_color = plt.cm.gist_rainbow(abs(phase/24.))        
        plt.plot(eff_waves_inv, phase_excesses, 's', color=mfc_color,
                 ms=8, mec='none', mfc=mfc_color, alpha=0.8)


    x = np.arange(3000,10000,10)
    xinv = 10000./x
    ftz_curve = redden_fm(x, np.zeros(x.shape), EBV, RV, return_excess=True)
    plt.plot(xinv, ftz_curve, 'k--')

    if plotpl:
        # plot PL curve
        gpl_curve = redden_pl(x, np.zeros(x.shape), AV, P, return_excess=True)
        plt.plot(xinv, gpl_curve, 'k-')
    
 
    ax.set_title(name+': Color Excess at B-maximum (with '+prefix[:-1]+' filters)')
    plt.ylabel('$E(V-X)$')
    plt.xlabel('Wavelength ($1 / \mu m$)')
    plt.xlim(1.0, 3.0)

def mksedplot( z=1.8, color='k' ):
    """ make a set of plots showing the SN Ia SED (Hsiao template)
    at various redshifts, with bandpasses overlaid.
    """
    snIa = sncosmo.Model( source='hsiao' )
    snIa.set( z=z, t0=0 )
    snwave = np.arange( 6000., 20000., 10. )
    snflux = snIa.flux( 0, snwave )
    snwave = snwave / 10000.
    snflux = 0.5 * snflux / snflux.max()

    pl.plot( snwave, snflux, color=color, ls='-')

    f105w = sncosmo.get_bandpass( 'f105w' )
    f098m = sncosmo.get_bandpass( 'f098m' )
    f127m = sncosmo.get_bandpass( 'f127m' )
    f139m = sncosmo.get_bandpass( 'f139m' )
    f153m = sncosmo.get_bandpass( 'f153m' )

    wf127m = f127m.wave / 10000.
    wf139m = f139m.wave / 10000.
    wf153m = f153m.wave / 10000.

    pl.plot( wf127m, f127m.trans, color='darkmagenta', ls='-')
    pl.plot( wf139m, f139m.trans, color='teal', ls='-')
    pl.plot( wf153m, f153m.trans, color='darkorange', ls='-')

    intf127m = scint.interp1d( wf127m, f127m.trans, bounds_error=False, fill_value=0 )
    overlap = np.min( [snflux, intf127m(snwave)], axis=0 )
    pl.fill_between( snwave, np.zeros(len(snwave)), overlap, color='darkmagenta' )

    intf139m = scint.interp1d( wf139m, f139m.trans, bounds_error=False, fill_value=0 )
    overlap = np.min( [snflux, intf139m(snwave)], axis=0 )
    pl.fill_between( snwave, np.zeros(len(snwave)), overlap, color='teal' )

    intf153m = scint.interp1d( wf153m, f153m.trans, bounds_error=False, fill_value=0 )
    overlap = np.min( [snflux, intf153m(snwave)], axis=0 )
    pl.fill_between( snwave, np.zeros(len(snwave)), overlap, color='darkorange' )

    def matchSNANAbandnamesinregistry():
        """
        Will have to build this along as we go, as I don't know the variety
        of naming conventions

        """

        bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
        for bandpass in bandPassList:
            band = sncosmo.get_bandpass('LSST_' + bandpass)
            band.name = bandpass
            if bandpass == 'y':
                band.name = 'Y'
            sncosmo.registry.register(band)

    def is_discover(self, band, z, sys, ep, peakmag=-18.4, sig_thresh=0.3, deep='No'):

            """
            INPUTS: Filter (rest frame), Redshift, Magnitude System, Epochs of observation 
            OPTIONS: Absolute Peak magnitude

            Outputs: array of observed magnitudes that are above the detection limit

            """
            
            input_filter = filtcov(z).obs_filt(band, z)[0]
            
            try:
                fcosm = simlc().create_bandpass(input_filter[0])
            except:
                fcosm = sncosmo.get_bandpass(band)
            mod = self.set_params(band, z, peakmag=peakmag)

            
            mag_arr=mod.bandmag(fcosm, sys, ep)

            filt_arr = np.array(self.filters)
            #if the deep fields limit is set use, else use dedicated survey limits from Astier+2014
            if deep == 'Yes':
                limarr = np.array(self.limits)
            elif deep == 'No':
                limarr = np.array(self.deep_limits)
            elif deep == 'Mod':
                limarr = np.array(self.mod_limits)
            elif deep == 'Ast':
                limarr = np.array(self.ast_limits)
            #extract the limiting magnitude for the appropriate filter    
            limmag = limarr[filt_arr == input_filter[0]]
            print limmag, mag_arr

            sig_eval = self.sigma(mag_arr)
            #strict threshold on the estimated error
            ##(Do something more sophisticated??)
            disc_arr = sig_eval[sig_eval <= sig_thresh]#mag_arr[mag_arr < float(limmag[0])]
            
            
            disc_arr = list(disc_arr)
            if not disc_arr:
                print "No Observation above the threshold"
            
                return []
            else:
                print "SN is discovered by Euclid"
                return list(disc_arr)

    def get_jwst_filt(self, indexx, pl='No'):
        """
        From the filters described above, select a JWST and optionally plot it
        Input: Numerical index for the array

        """
        filt_val=self.flist[indexx]
        band = sncosmo.get_bandpass(filt_val)
        
        
        if pl=='Yes':
            plt.plot(band.wave, band.trans)
            plt.show()
        else:
            return band

def test_bandpass_bessell():
    """Check that Bessell bandpass definitions are scaled by inverse
    wavelength."""

    band = sncosmo.get_bandpass('bessellb')
    trans = band.trans[[4, 9, 14]]  # transmission at 4000, 4500, 5000

    # copied from file
    orig_wave = np.array([4000., 4500., 5000.])
    orig_trans = np.array([0.920, 0.853, 0.325])

    scaled_trans = orig_trans / orig_wave

    # scaled_trans should be proportional to trans
    factor = scaled_trans[0] / trans[0]
    assert_allclose(scaled_trans, factor * trans)

def lsq_excess_fit(ref_excess_dict, red_law, EBV, rv_guess, filters, zp):
    prefix = zp['prefix']
    filter_eff_waves = np.array([snc.get_bandpass(prefix+f).wave_eff for f in filters])
    ref_excess = np.array([ref_excess_dict[f] for f in filters])
    
    def lsq_func(Y):
        RV = Y[0]
        ftz_curve = red_law(filter_eff_waves,
                              np.zeros(filter_eff_waves.shape),
                              -EBV, RV,
                              return_excess=True)
        return ftz_curve-ref_excess
    
    Y = np.array([rv_guess])
    valid_phases = {}
    return lsq(lsq_func, Y)

def which_salt(z):
    salt_name = 'salt2'
    salt_version = '2.4'

    rest_salt_max_wav = 9200
    rest_salt_min_wav = 2000

    salt_max_wav = (1 + z) * rest_salt_max_wav
    salt_min_wav = (1 + z) * rest_salt_min_wav

    for filter in filters:
        band = sncosmo.get_bandpass(filter)
        if (band.wave[0] < salt_min_wav or band.wave[-1] > salt_max_wav):
            salt_name = 'salt2-extended'
            salt_version = '1.0'
            break
    return salt_name, salt_version

def Gen_SN():
    # Use below if on Iridis
    # source = sncosmo.SALT2Source(modeldir="/scratch/cf5g09/Monte_Carlos/salt2-4")

    ##Use below if not on iridis
    source = sncosmo.get_source('salt2', version='2.4')
    model = sncosmo.Model(source=source)
    timearray = range(100)
    mabs = -19.05  # Setting the absolute Magnitude of the Supernova
    model.set(z=0.01, t0=30)  # Setting redshift

    model.set_source_peakabsmag(mabs, 'ptf48r', 'ab',
                                cosmo=FlatLambdaCDM(H0=70, Om0=0.3))  # Fixing my peak absolute magnitude
    # model.set(x1=x_1, c=colour)

    absmag_r = model.source_peakabsmag('ptf48r', 'ab', cosmo=FlatLambdaCDM(H0=70, Om0=0.3))
    print absmag_r

    band = sncosmo.get_bandpass('ptf48r')  # Retrieving the ptf48r bandpass

    maglc = model.bandmag('ptf48r', 'ab', timearray)  # Creating a magnitude array of the lightcurve
    fluxlc = model.bandflux('ptf48r', timearray)  # Creating a flux array of the lightcurve

    model2 = model
    model2.set_source_peakabsmag(mabs, 'bessellr', 'ab',
                                 cosmo=FlatLambdaCDM(H0=70, Om0=0.3))  # Fixing my peak absolute magnitude
    # model.set(x1=x_1, c=colour)

    absmag_r = model2.source_peakabsmag('bessellr', 'ab', cosmo=FlatLambdaCDM(H0=70, Om0=0.3))
    maglc2 = model2.bandmag('bessellr', 'ab', timearray)  # Creating a magnitude array of the lightcurve

    fluxlc2 = model2.bandflux('bessellr', timearray)  # Creating a magnitude array of the lightcurve

    plt.scatter(timearray, fluxlc, color='blue', label='ptf48r')
    plt.scatter(timearray, fluxlc2, color='red', label='bessellr')

    model.bandflux('PTF48R', 30)

    # plt.gca().invert_yaxis()
    plt.title('SNTypeIa peak 30 days')
    plt.legend()
    plt.show()

    # print sn_par

    return maglc, fluxlc, timearray

def lightcurve_Ia(filter, z, x1, c, x0=None):
    """Given a filter and redshift z, generates the observed
    flux for SNe Type Ia"""
    alpha = 0.12
    beta = 3.
    mabs = -19.1 - alpha*x1 + beta*c
    zp = zero_point[filter]
    zpsys = zero_point['zpsys']

    # Checking if bandpass is outside spectral range for SALT2. If yes,
    # use salt2-extended.
    salt_name = 'salt2'
    salt_version = '2.4'

    rest_salt_max_wav = 9200
    rest_salt_min_wav = 2000

    salt_max_wav = (1 + z) * rest_salt_max_wav
    salt_min_wav = (1 + z) * rest_salt_min_wav

    band = sncosmo.get_bandpass(filter)
    if (band.wave[0] < salt_min_wav or band.wave[-1] > salt_max_wav):
        salt_name = 'salt2-extended'
        salt_version = '1.0'

    # Type Ia model
    model_Ia = sncosmo.Model(source=sncosmo.get_source(salt_name,
                                                       version=salt_version))
    if x0 is not None:
        p = {'z': z, 't0': t0, 'x0': x0, 'x1': x1,
             'c': c}
    else:
        p = {'z': z, 't0': t0, 'x1': x1, 'c': c}
        model_Ia.set(z=z)
        model_Ia.set_source_peakabsmag(mabs, 'bessellb', 'vega')
    model_Ia.set(**p)
    phase_array = np.linspace(model_Ia.mintime(), model_Ia.maxtime(), 100)
    obsflux_Ia = model_Ia.bandflux(filter, phase_array, zp=zp, zpsys=zpsys)
    keys = ['phase_array', 'obsflux']
    values = [phase_array, obsflux_Ia]
    dict_Ia = dict(zip(keys, values))
    np.savetxt('test.dat', np.c_[dict_Ia['phase_array'], dict_Ia['obsflux']])
    x0 = model_Ia.get('x0')
    return (dict_Ia, x0, salt_name, salt_version)

def getlsstbandpassobjs(plot=False):
    """
    General utility to return a list of the baseline LSST bandpasses
    and register them as SNCosmo bandpasses accessible through strings 
    like 'LSSTu'.
    
    Parameters
    ----------
    plot : bool, optional, defaults to False
        plot filter functions obtained
    Returns
    -------
    None
    
    Examples
    --------
    >>> getlsstbandpassobjs() 

    """
    bandPassList = ["u", "g", "r", "i", "z", "y"]
    banddir = os.path.join(os.getenv("THROUGHPUTS_DIR"), "baseline")
    lsstbands = []
    lsstbp = {}

    for band in bandPassList:
        # setup sncosmo bandpasses
        bandfname = banddir + "/total_" + band + ".dat"
        loadsncosmo = True
        if loadsncosmo:
            # register the LSST bands to the SNCosmo registry
            # Not needed for LSST, but useful to compare independent codes
            # Usually the next two lines can be merged,
            # but there is an astropy bug currently which affects only OSX.
            numpyband = np.loadtxt(bandfname)
            sncosmoband = sncosmo.Bandpass(
                wave=numpyband[:, 0], trans=numpyband[:, 1], wave_unit=Unit("nm"), name="LSST_" + band
            )
            sncosmo.registry.register(sncosmoband, force=True)
    if plot:
        filterfigs, filterax = plt.subplots()
        for band in bandPassList:
            b = sncosmo.get_bandpass("LSST_" + band)
            filterax.plot(b.wave, b.trans, "-k", lw=2.0)
    return None