def model_IIP_nugent(self):
        """
        """
        sncosmo.get_source('nugent-sn2p', version='1.2')
        p, w, f = sncosmo.read_griddata_ascii(
            os.path.join(sncosmo.builtins.get_cache_dir(),
                         'sncosmo/models/nugent/sn2p_flux.v1.2.dat')
        )
        mask = (p < 130)

        return sncosmo.Model(source=sncosmo.TimeSeriesSource(p[mask], w, f[mask]),
                             effects=[sncosmo.CCM89Dust()],
                             effect_names=['host'],
                             effect_frames=['rest'])

    def model_Ia_hsiao(self):
        """
        """
        sncosmo.get_source('hsiao', version='3.0')
        p, w, f = sncosmo.read_griddata_fits(
            os.path.join(sncosmo.builtins.get_cache_dir(),
                         'sncosmo/models/hsiao/Hsiao_SED_V3.fits')
        )

        return sncosmo.Model(
            source=sncosmo.StretchSource(p, w, f, name='hsiao-stretch'),
            effects=[sncosmo.CCM89Dust()],
            effect_names=['host'],
            effect_frames=['rest']
        )

def showcurve(sn, source_name):
    try:
        hml=np.load('/home/fcm1g13/Documents/Supernova/CorecollapseLC/ccdata/' + sn)    
    
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
    # print len(hml[:,1]), 'data points'
    
        # Adding zpsys and filter columns. zp system used is ab and filter used is ptf48r
        ab = np.zeros(len(hml[:, 1]), dtype='|S2')
        band = np.zeros(len(hml[:, 1]), dtype='|S6')
        for i in range(len(ab)):
            ab[i] = 'ab'
            band[i] = 'ptf48r'
        hml = np.column_stack((hml, ab, band))

        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        
        ###fitting model
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(0.005,0.35)}, nburn=10000, nsamples=50000)
            #print 'peak', float(res.parameters[1])
        print res                              
        sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name), xfigsize=10)
        #plt.clfig()
        plt.show()

        #print #'### Parameters ###'
        #print 'SN',str(hml[:,0][0]), 'z:',float(res.parameters[0])# float(res.errors['z']), float(res.parameters[1]), float(res.errors['t0']),float(res.parameters[2]), float(res.errors['x0']),  float(res.parameters[3]), float(res.errors['x1']), float(res.parameters[4]), float(res.errors['c']), float(hml[:,8][0]), float(hml[:,9][0])
        print 'Done:', hml[:,0][0], 'z:',float(res.parameters[0]), 'Reduced chi^2:', res.chisq/res.ndof,#'Data points:', len(hml[:,1]),' Type'+types[hml[:,0][0]]+'Model name: '+ source.name,'\n' #'Dof:', res.ndof

    except ValueError:
        sn, source_name, 'cannot be plotted'

def lightcurve_Ibc(filter, z, hostebv_Ibc):
    """Given a filter, redshift z at given phase, generates the observed
    magnitude for SNe Type Ibc"""
    zp = zero_point[filter]
    zpsys = zero_point['zpsys']
    model_Ibc = ['s11-2005hl', 's11-2005hm', 's11-2006fo', 'nugent-sn1bc',
                 'nugent-hyper', 's11-2006jo', 'snana-2004fe',
                 'snana-2004gq', 'snana-sdss004012', 'snana-2006fo',
                 'snana-sdss014475', 'snana-2006lc', 'snana-04d1la',
                 'snana-04d4jv', 'snana-2004gv', 'snana-2006ep',
                 'snana-2007y', 'snana-2004ib', 'snana-2005hm',
                 'snana-2006jo', 'snana-2007nc']
    obsflux_Ibc = []
    phase_arrays = []
    for i in model_Ibc:
        model_i = sncosmo.Model(source=sncosmo.get_source(i), effects=[dust],
                                effect_names=['host'], effect_frames=['rest'])
        mabs = -17.56
        model_i.set(z=z)
        phase_array_i = np.linspace(model_i.mintime(), model_i.maxtime(), 100)
        model_i.set_source_peakabsmag(mabs, 'bessellb', 'ab')
        p_core_collapse = {'z': z, 't0': t0, 'hostebv': hostebv_Ibc,
                           'hostr_v': hostr_v}
        model_i.set(**p_core_collapse)
        phase_arrays.append(phase_array_i)
        obsflux_i = model_i.bandflux(filter, phase_array_i, zp, zpsys)
        obsflux_Ibc.append(obsflux_i)
    keys = model_Ibc
    values = []
    for i, item in enumerate(model_Ibc):
        values.append([obsflux_Ibc[i], phase_arrays[i]])
    dict_Ibc = dict(zip(keys, values))
    return (dict_Ibc)

def fitcurve(x, source_name, hml):
    try:
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
        
        #adding zpsys and filter columns
        ab=np.zeros(len(hml[:,1]), dtype='|S2')
        for i in range(len(ab)):
            ab[i]='ab'
        hml=np.column_stack((hml,ab))   
        band=np.zeros(len(hml[:,1]), dtype='|S6')
        for i in range(len(band)):
            band[i]='ptf48r'
        hml=np.column_stack((hml,band))
        #fit to model
        z0 = float(spec_z(x[:-7]))
        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(z0,z0 + 0.0001)}, nburn=10000, nsamples=50000)
        
        #The following excludes data points not in the range of the model and data sets with fewer than 4 data points
        limit = modellimit(source_name, x[:-7], res.parameters[1])
        hml2 = []
        for j in range(len(hml[:,1])):
            datapoint = hml [:,1][j]
            if (res.parameters[1]- limit[0])< float(datapoint) < (res.parameters[1]+limit[1]):
                hml2.append(hml[j]) 
        hml2 = np.array(hml2)
        if len(hml2)>3:
            
            return finalfitcurve(x, source_name, hml2) 
               
        
    except ValueError:
        print 'error'  

def select_model(x):
    # Load file containing data set
    hml = np.load(os.path.abspath('') + '/ccdata/' + x)
    # print len(hml[:,1]), 'data points'             #number of rows and therefore data points

    # Adding zpsys and filter columns. zp system used is ab and filter used is ptf48r
    ab = np.zeros(len(hml[:, 1]), dtype='|S2')
    band = np.zeros(len(hml[:, 1]), dtype='|S6')
    for i in range(len(ab)):
        ab[i] = 'ab'
        band[i] = 'ptf48r'
    hml = np.column_stack((hml, ab, band))

    # Converting into a table using astropy with titles: ptfname, time,
    # magnitude, mag_err, flux, flux_err, zp_new, zp, ra, dec, zpsys and filter
    hml_dat = astropy.table.Table(data=hml, names=(
            'ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys',
            'filter'),
                                      dtype=(
                                          'str', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
                                            'float', 'str', 'str'))
    #sn types are Ib, Ib/c, Ic, Ic-BL, IIP
    sn_type = types[hml[:, 0][0]]

    # Selecting only Type II supernovae
    if sn_type[1] == 'I':
        # Selecting model
        source = sncosmo.get_source('s11-2004hx',version='1.0')
        fit_curve(hml_dat, source, hml)

    # Selecting only Type Ic supernovae
    elif sn_type[1]== 'c':
        # Selecting model
        source = sncosmo.get_source('snana-2004fe',version='1.0')
        fit_curve(hml_dat, source, hml)

    # Selecting TypeIb sn only
    elif len(sn_type)== 2 and (sn_type[1]!= 'b'):
        # Selecting model
        source = sncosmo.get_source('snana-2006ep', version='1.0')
        fit_curve(hml_dat, source, hml)

    # Selecting TypeIbc sn only
    elif len(sn_type)> 2 and (sn_type[1]!= 'b'):
        # Selecting model
        source = sncosmo.get_source('nugent-sn1bc', version='1.1')
        fit_curve(hml_dat, source, hml)

def Gen_SN():
    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, 'bessellb', 'ab',
                                cosmo=FlatLambdaCDM(H0=70, Om0=0.3))  # Fixing my peak absolute magnitude
    # model.set(x1=x_1, c=colour)
    '''
    absmagb =model.source_peakabsmag('bessellb','ab', cosmo=FlatLambdaCDM(H0=70,Om0=0.3))
    absmag_r =model.source_peakabsmag('bessellr','ab', cosmo=FlatLambdaCDM(H0=70,Om0=0.3))


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

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

    data = np.loadtxt('/home/fcm1g13/Documents/Supernova/PTF48R filter/PTF48R.dat')
    wavelength = np.array([row[0] for row in data])
    transmission = np.array([row[1] for row in data])

    band = sncosmo.Bandpass(wavelength, transmission, name='PTF48R')

    wave = np.arange(5580., 7500., 20)
    spec = model.flux([20., 30., 35., 40., 50.], wave)

    adjspec = transmission * spec[1]
    '''
        for point in transmission:
            for flux in spec[1]:
                adjspec.append(point*flux)
        '''
    print len(transmission)
    print len(spec[1])
    print adjspec

    plt.plot(wave, spec[1], color='#27ae60', label='Flux')
    model.bandflux('PTF48R', 30)
    plt.plot(wave, np.array(adjspec), color='#2980b9', label='Observed flux')
    plt.plot(wavelength, transmission * max(spec[1]))
    # plt.plot(wave, spec[2], color = '#27ae60',label = '5 days after peak')
    # plt.plot(wave, spec[3], color = '#c0392b',label = '10 days after peak')
    # plt.plot(wave, spec[4], color = '#8e44ad',label = '20 days after peak')

    plt.title('Model spectrum of Type Ia supernova at peak')
    plt.ylabel('Flux in ergs / s / cm^2 / $\AA$')
    plt.xlabel('Wavelength in $\AA$')
    plt.minorticks_on()
    plt.ylim([0, max(spec[1]) * 1.2])
    plt.legend()
    plt.draw()

    return timearray

    def luminosity(**kwargs):
        source = sncosmo.get_source("nugent-sn1bc")
        source.set_peakmag(-17.5 + kwargs["x0"], "desg", "ab")  # magnitude at 10pc
        dust = sncosmo.CCM89Dust()
        model = sncosmo.Model(source=source, effects=[dust], effect_names=["host"], effect_frames=["rest"])
        kwargs2 = dict(kwargs)
        del kwargs2["x0"]

        model.set(**kwargs2)
        return model

	def luminosity(**kwargs):
		source = sncosmo.get_source("nugent-sn1bc")
		source.set_peakmag(-17.5+kwargs['x0'], 'desg', 'ab') # magnitude at 10pc	
		dust = sncosmo.CCM89Dust()
		model = sncosmo.Model(source=source,effects=[dust], effect_names=['host'], effect_frames=['rest'])
		kwargs2 = dict(kwargs)
		del kwargs2['x0']

		model.set(**kwargs2)
		return model

def get_snana_filenames(sntype):
    """
    """
    if not sntype.startswith('SN '):
        sntype = 'SN %s'%sntype

    reg = sncosmo.registry._get_registry(sncosmo.Source)
    source_tuples = [(v['name'], v['version'])
                     for v in reg.get_loaders_metadata()
                     if v['name'].startswith('snana')
                     and v['type'] == sntype]

    filenames = []
    for name, version in source_tuples:
        filepath = os.path.join(sncosmo.builtins.get_cache_dir(),
                                'sncosmo',
                                reg._loaders[(name, version)][1])
        if not os.exists(filepath):
            sncosmo.get_source(name, version=version)
        filenames.append(filepath)

    return filenames

def model_limits(sou):
    source = sncosmo.get_source(sou)
    model = sncosmo.Model(source=source)

    timearray = range(-20, 300)
    fluxlc = model.bandflux('ptf48r',timearray)     # Creating a flux array of the light curve
    modelbreak = []
    for i in range(len(fluxlc)-1):
        if fluxlc[i] == fluxlc[i+1]:
            modelbreak.append(timearray[i+1])
    for i in range(len(modelbreak)-1):
        if (modelbreak[i+1]-modelbreak[i]) > 20:
            return sou, modelbreak[i], modelbreak[i+1]

def fitcurve(x, source_name):
    # Get data set needed
    hml=np.load('/home/fcm1g13/Documents/Supernova/CorecollapseLC/ccdata/' + x)    
    try:
        # Import model
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
       # print len(hml[:,1]), 'data points'
       
        # Adding zpsys and filter columns. zp system used is ab and filter used is ptf48r
        ab = np.zeros(len(hml[:, 1]), dtype='|S2')
        band = np.zeros(len(hml[:, 1]), dtype='|S6')
        for i in range(len(ab)):
            ab[i] = 'ab'
            band[i] = 'ptf48r'
        hml = np.column_stack((hml, ab, band))

        # Converting into a table using astropy with titles: ptfname, time,
        # magnitude, mag_err, flux, flux_err, zp_new, zp, ra, dec, zpsys and filter
        hml_dat = astropy.table.Table(data=hml, names=(
             'ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys',
               'filter'),
                                      dtype=(
                                          'str', 'float', 'float', 'float', 'float', 'float', 'float', 'float', 'float',
                                            'float', 'str', 'str'))

        # Fitting model: model parameter z bound.
        res, fitted_model = sncosmo.fit_lc(hml_dat, model, ['z', 't0', 'amplitude'], bounds={'z': (0.005, 0.35)})

        #The following excludes data with fewer than 4 data points and not enough distribution..
        j=0; k=0
        for i in hml[:,1]:
            if float(i)>float(res.parameters[1]):
                j+=1
            else:
                k+=1
        # print hml[:,0][0],k,j
        if j>=2 and k>=2:            
            if len(hml[:,1])>3:                                
                #sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name), xfigsize=10)
                #plt.close()                
                
                return np.array([float(res.chisq/res.ndof), hml[:,0][0], source_name]) #returns reduced chi squared, sn name and model name. 
            else:
                pass
        else:
            pass
                
    except ValueError:
        pass

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 showcurve(sn, source_name, hml):
    try:
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 

        #fit to model
        z0 = float(spec_z(sn[:-7]))
        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(z0, z0+0.001)}, nburn=10000, nsamples=50000)

        #plot model
        sncosmo.plot_lc(hml_dat, model=fitted_model, errors=res.errors, color='blue', figtext=str(hml[:,0][0]+' Type'+types[hml[:,0][0]]+'\n'+'Model name: '+ source.name + '\n'+'Reduced Chi Squared: ')+ str(res.chisq/res.ndof), xfigsize=10)
        plt.show()

        #print 'Parameters:''z:',float(res.parameters[0]), float(res.errors['z']), float(res.parameters[1]), float(res.errors['t0']),float(res.parameters[2]), float(res.errors['x0']),  float(res.parameters[3]), float(res.errors['x1']), float(res.parameters[4]), float(res.errors['c']), float(hml[:,8][0]), float(hml[:,9][0])'
        print 'Done:', hml[:,0][0], 'z:',float(res.parameters[0]), 'Reduced chi^2:', res.chisq/res.ndof,#'Data points:', len(hml[:,1]),' Type'+types[hml[:,0][0]]+'Model name: '+ source.name,'\n' #'Dof:', res.ndof

    except ValueError:
        print sn, source_name, 'cannot be plotted'

def get_model_cc(z, sn_type, all_model):
    """Given a filter and redshift z, generates the observed
    magnitude for SNe Type Ibc or Type II"""
    my_model = choice(all_model)
    if my_model in ['s11-2004hx', 'nugent-sn2l']:
        sn_type = 'IIl'
    elif my_model in ['nugent-sn2n', 'snana-2006ez', 'snana-20069ix']:
        sn_type = 'IIn'
    model = sncosmo.Model(source=sncosmo.get_source(my_model),
                          effects=[dust], effect_names=['host'],
                          effect_frames=['rest'])
    mabs = all_mabs(sn_type)
    model.set(z=z)
    model.set_source_peakabsmag(mabs, 'bessellb', 'vega')
    p = {'z': z, 't0': t0, 'hostebv': uniform(-0.1, 0.65), 'hostr_v': hostr_v}
    model.set(**p)
    p['model_name'] = my_model
    max_phase = max_cc_phase(my_model)
    my_phase = uniform(-2, 3) + max_phase
    return model, p, my_phase

def lightcurve_II(filter, z, hostebv_II):
    """Given a filter and redshift z, generates the observed magnitude for
    SNe Type II"""
    zp = zero_point[filter]
    zpsys = zero_point['zpsys']
    model_II = ['s11-2005lc', 's11-2005gi', 's11-2006jl', 'nugent-sn2p',
                'snana-2004hx', 'snana-2005gi', 'snana-2006gq',
                'snana-2006kn', 'snana-2006jl', 'snana-2006iw',
                'snana-2006kv', 'snana-2006ns', 'snana-2007iz',
                'snana-2007nr', 'snana-2007nr', 'snana-2007kw',
                'snana-2007ky', 'snana-2007lj', 'snana-2007lb',
                'snana-2007ll', 'snana-2007nw', 'snana-2007ld',
                'snana-2007md', 'snana-2007lz', 'snana-2007lx',
                'snana-2007og', 'snana-2007ny', 'snana-2007nv',
                'snana-2007pg', 's11-2004hx', 'nugent-sn2l', 'nugent-sn2n',
                'snana-2006ez', 'snana-2006ix']
    obsflux_II = []
    phase_arrays = []
    for i in model_II:
        model_i = sncosmo.Model(source=sncosmo.get_source(i), effects=[dust],
                                effect_names=['host'],
                                effect_frames=['rest'])
        if i == 's11-2004hx' == 'nugent-sn2l':
            mabs = -17.98
        else:
            mabs = -16.75
        model_i.set(z=z)
        phase_array_i = np.linspace(model_i.mintime(), model_i.maxtime(), 100)
        model_i.set_source_peakabsmag(mabs, 'bessellb', 'ab')
        p_core_collapse = {'z': z, 't0': t0, 'hostebv': hostebv_II,
                           'hostr_v': hostr_v}
        model_i.set(**p_core_collapse)
        phase_arrays.append(phase_array_i)
        obsflux_i = model_i.bandflux(filter, phase_array_i, zp, zpsys)
        obsflux_II.append(obsflux_i)
    keys = model_II
    values = []
    for i, item in enumerate(model_II):
        values.append([obsflux_II[i], phase_arrays[i]])
    dict_II = dict(zip(keys, values))
    return (dict_II)

def get_model_Ia(z, min_phase, max_phase):
    """Given a filter, redshift z at given phase, generates the observed
    magnitude for SNe Type Ia"""

    alpha = 0.12
    beta = 3.
    x1 = normal(0., 1.)
    c = normal(0., 0.1)
    my_phase = uniform(min_phase, max_phase)
    mabs = normal(-19.1 - alpha*x1 + beta*c, scale=0.15)
    salt_name, salt_version = which_salt(z)
    model_Ia = sncosmo.Model(source=sncosmo.get_source(salt_name,
                                                       version=salt_version))
    model_Ia.set(z=z)
    model_Ia.set_source_peakabsmag(mabs, 'bessellb', 'vega')
    p = {'z': z, 't0': t0, 'x1': x1, 'c': c}
    p['x0'] = model_Ia.get('x0')
    model_Ia.set(**p)
    p['salt_name'] = salt_name
    p['salt_version'] = salt_version
    return model_Ia, p, my_phase

def modellimit(source_name, x, t0):
    source=sncosmo.get_source(source_name) 
    model=sncosmo.Model(source=source)
    
    #Creating a flux array of the model lightcurve
    timearray= np.arange(t0-120,t0+500, 0.01)
    model.set(z=spec_z(x),t0=t0)
    fluxlc=model.bandflux('ptf48r',timearray) 
    modelbreak = [t0-120.0]
    for i in range(len(fluxlc)-1):
        if fluxlc[i] == fluxlc[i+1]:
            modelbreak.append(timearray[i+1])

    #find the limits of the model by defining where the flux does not change. 
    for i in range(len(modelbreak)-1):
        if abs(modelbreak[i+1]-modelbreak[i]) > 20:
            minimum = t0 - modelbreak[i]
            maximum = modelbreak[i+1] - t0
            
            #return the difference of t0 and the edges of the model and the time of peak flux
            ##this is done as some models use t0 as peak and others use it for start of model
            return minimum, maximum, timearray[list(fluxlc).index(max(fluxlc))]

def finalfitcurve(x, source_name, hml):
    try:
        source=sncosmo.get_source(source_name) 
        model=sncosmo.Model(source=source) 
        
        #fit to model
        z0 = float(spec_z(x[:-7]))
        hml_dat=astropy.table.Table(data=hml, names=('ptfname', 'time', 'magnitude', 'mag_err', 'flux', 'flux_err', 'zp_new', 'zp', 'ra', 'dec', 'zpsys', 'filter'), dtype=('str','float','float','float','float','float','float','float','float','float','str', 'str'))
        res, fitted_model=sncosmo.fit_lc(hml_dat, model, ['z','t0','amplitude'], bounds={'z':(z0, z0+0.0001)}, nburn=10000, nsamples=50000)
        #The following excludes data with not enough distribution..
        ##Ensuring at least 2 data points on either side of peak.
        limit = modellimit(source_name, x[:-7], res.parameters[1])
        j=0; k=0
        for i in hml[:,1]:
            if float(i)>float(limit[2]):
                j+=1
            else:
                k+=1
        if j>=2 and k>=2:          
            return np.array([float(res.chisq/res.ndof), hml[:,0][0], source_name]), hml #returns reduced chi squared, sn name and model name. 
           
    except ValueError:
        print 'error'