# -*- coding: utf-8 -*-
from __future__ import print_function, division

"""
.. note::
         These are the spectrophotometry functions for SPLAT 
"""

# imports - internal
import copy
import os

# imports - external
import numpy
from astropy import units as u            # standard units
from astropy import constants as const        # physical constants in SI units
import matplotlib.patches as patches
import matplotlib.pyplot as plt
from scipy.integrate import trapz        # for numerical integration
from scipy.interpolate import interp1d

# splat functions and constants
from .initialize import *
from .utilities import *



#####################################################
###############   SPECTROPHOTOMETRY   ###############
#####################################################

# this function has been obseleted
def checkFilter(filt,verbose=True):
    output = False
    f = copy.deepcopy(filt)
    f = f.replace(' ','_').upper()
    for k in list(FILTERS.keys()):
        if f==k.upper() or f.lower() in FILTERS[k]['altnames']:
            output = k
    if output == False and verbose == True: 
        print('\nFilter '+filt+' not currently available for SPLAT; contact '+EMAIL+'\n')
        filterInfo()
    return output


def filterProfile(filt,**kwargs):
    '''
    :Purpose: Retrieve the filter profile for a SPLAT filter. Returns two arrays: the filter wavelength and filter transmission curve.

    :param filter: String giving the name of one of the predefined filters listed in splat.FILTERS.keys() (required)
    
    :param filterFolder: folder containing the filter transmission files (optional, default = splat.FILTER_FOLDER)

    :Example:
    >>> import splat
    >>> import splat.photometry as spphot
    >>> sp = splat.getSpectrum(shortname='1507-1627')[0]
    >>> sp.fluxCalibrate('2MASS J',14.5)
    >>> spphot.filterMag(sp,'MKO J')
        (14.345894376898123, 0.027596454828421831)
    '''
# keyword parameters
    filterFolder = kwargs.get('filterFolder',SPLAT_PATH+FILTER_FOLDER)
    if not os.path.exists(filterFolder):
        filterFolder = SPLAT_URL+FILTER_FOLDER

# check that requested filter is in list
    f0 = checkFilterName(filt, verbose=True)
    if f0 == False: raise ValueError
    filt = f0

# read in filter
    fwave,ftrans = numpy.genfromtxt(os.path.normpath(filterFolder+FILTERS[filt]['file']), comments='#', unpack=True, missing_values = ('NaN','nan'), filling_values = (numpy.nan))
#    print(type(fwave),type(ftrans),isinstance(fwave,numpy.ndarray),isinstance(ftrans,numpy.ndarray),not isinstance(fwave,numpy.ndarray) or not isinstance(ftrans,numpy.ndarray))
    if not isinstance(fwave,numpy.ndarray) or not isinstance(ftrans,numpy.ndarray):
        raise ValueError('\nProblem reading in {}'.format(filterFolder+FILTERS[filt]['file']))
    fwave = fwave[~numpy.isnan(ftrans)]*u.micron
    ftrans = ftrans[~numpy.isnan(ftrans)]
    return fwave,ftrans



def filterMag(sp,filt,*args,**kwargs):
    '''
    :Purpose: 

    Determine the photometric magnitude of a source based on its
    spectrum. Spectral fluxes are convolved with the filter profile specified by
    the ``filter`` input.  By default this filter is also
    convolved with a model of Vega to extract Vega magnitudes,
    but the user can also specify AB magnitudes, photon flux or energy flux.

    :Required Parameters:

        **sp**: Spectrum class object, which should contain wave, flux and noise array elements.
        **filter**: String giving name of filter, which can either be one of the predefined filters listed in splat.FILTERS.keys() or a custom filter name

    :Optional Parameters:
    
        **custom** = None: A 2 x N vector array specifying the wavelengths and transmissions for a custom filter
        **notch** = None: A 2 element array that specifies the lower and upper wavelengths for a notch filter (100% transmission within, 0% transmission without)
        **vega** = True: compute Vega magnitudes (may be set by filter)
        **ab** = False: compute AB magnitudes (may be set by filter)
        **energy** = False: compute energy flux
        **photon** = False: compute photon flux
        **filterFolder** = splat.FILTER_FOLDER: folder containing the filter transmission files
        **vegaFile** = 'vega_kurucz.txt': name of file containing Vega flux file, must be within ``filterFolder``
        **nsamples** = 100: number of samples to use in Monte Carlo error estimation
        **info** = False: List the predefined filter names available
        **verbose** = True: List the predefined filter names available

    :Example:
    >>> import splat
    >>> import splat.photometry as spphot
    >>> sp = splat.getSpectrum(shortname='1507-1627')[0]
    >>> sp.fluxCalibrate('2MASS J',14.5)
    >>> spphot.filterMag(sp,'MKO J')
        (14.345894376898123, 0.027596454828421831)
    '''
# keyword parameters
    filterFolder = kwargs.get('filterFolder',SPLAT_PATH+FILTER_FOLDER)
    if not os.path.exists(filterFolder):
        filterFolder = SPLAT_URL+FILTER_FOLDER
    vegaFile = kwargs.get('vegaFile',VEGAFILE)
    info = kwargs.get('info',False)
    custom = kwargs.get('custom',False)
    notch = kwargs.get('notch',False)
    vega = kwargs.get('vega',True)
    ab = kwargs.get('ab',not vega)
    rsr = kwargs.get('rsr',False)
    nsamples = kwargs.get('nsamples',100)
    verbose = kwargs.get('verbose',False)


# check that requested filter is in list
    if isinstance(custom,bool) and isinstance(notch,bool):
        f0 = checkFilterName(filt,verbose=True)
        if f0 == False: 
            return numpy.nan, numpy.nan
        filt = f0
        
# reset filter calculation methods based on filter design
        if 'ab' in FILTERS[filt]['method']: 
            ab = kwargs.get('ab',True)
            vega = not ab
        if 'vega' in FILTERS[filt]['method']: 
            vega = kwargs.get('vega',True)
            ab = not vega
        rsr = FILTERS[filt]['rsr']

# other possibilities
    photons = kwargs.get('photons',False)
    photons = kwargs.get('photon',photons)
    energy = kwargs.get('energy',False)
    energy = kwargs.get('flux',energy)
    if (photons or energy):
        vega = False
        ab = False
    if photons: energy = False
    if energy: photons = False


# Read in filter
    if isinstance(custom,bool) and isinstance(notch,bool):
        fwave,ftrans = filterProfile(filt,**kwargs)
# notch filter
    elif isinstance(custom,bool) and isinstance(notch,list):
        dn = (notch[1]-notch[0])/1000
        fwave = numpy.arange(notch[0]-5.*dn,notch[1]+5.*dn,dn)
        ftrans = numpy.zeros(len(fwave))
        ftrans[numpy.where(numpy.logical_and(fwave >= notch[0],fwave <= notch[1]))] = 1.
# custom filter
    else:
        fwave,ftrans = custom[0],custom[1]

# units
    if isinstance(fwave,u.quantity.Quantity) == True:
        fwave = fwave.to(u.micron)
    else:
        fwave = fwave*u.micron

# check that spectrum and filter cover the same wavelength ranges
    if numpy.nanmax(fwave) < numpy.nanmin(sp.wave) or numpy.nanmin(fwave) > numpy.nanmax(sp.wave):
        if verbose==True: print('\nWarning: no overlap between spectrum for {} and filter {}'.format(sp.name,filt))
        return numpy.nan, numpy.nan

    if numpy.nanmin(fwave) < numpy.nanmin(sp.wave) or numpy.nanmax(fwave) > numpy.nanmax(sp.wave):
        if verbose==True: print('\nWarning: spectrum for {} does not span full filter profile for {}'.format(sp.name,filt))

# interpolate spectrum onto filter wavelength function
    wgood = numpy.where(~numpy.isnan(sp.noise))
    if len(sp.wave[wgood]) > 0:
        d = interp1d(sp.wave[wgood].value,sp.flux[wgood].value,bounds_error=False,fill_value=0.)
        n = interp1d(sp.wave[wgood].value,sp.noise[wgood].value,bounds_error=False,fill_value=0)
# catch for models
    else:
        if verbose==True: print('\nWarning: data values in range of filter {} have no uncertainties'.format(filt))
        d = interp1d(sp.wave.value,sp.flux.value,bounds_error=False,fill_value=0.)
        n = interp1d(sp.wave.value,sp.flux.value*1.e-9,bounds_error=False,fill_value=0.)

    result = []
    if (vega):
# Read in Vega spectrum
        vwave,vflux = numpy.genfromtxt(os.path.normpath(filterFolder+vegaFile), comments='#', unpack=True, \
            missing_values = ('NaN','nan'), filling_values = (numpy.nan))
        vwave = vwave[~numpy.isnan(vflux)]*u.micron
        vflux = vflux[~numpy.isnan(vflux)]*(u.erg/(u.cm**2 * u.s * u.micron))
        vflux.to(sp.flux_unit,equivalencies=u.spectral_density(vwave))
# interpolate Vega onto filter wavelength function
        v = interp1d(vwave.value,vflux.value,bounds_error=False,fill_value=0.)
        if rsr:
            val = -2.5*numpy.log10(trapz(ftrans*fwave.value*d(fwave.value),fwave.value)/trapz(ftrans*fwave.value*v(fwave.value),fwave.value))
        else:
            val = -2.5*numpy.log10(trapz(ftrans*d(fwave.value),fwave.value)/trapz(ftrans*v(fwave.value),fwave.value))
        for i in numpy.arange(nsamples):
#            result.append(-2.5*numpy.log10(trapz(ftrans*numpy.random.normal(d(fwave),n(fwave))*sp.flux_unit,fwave)/trapz(ftrans*v(fwave)*sp.flux_unit,fwave)))
            if rsr:
                result.append(-2.5*numpy.log10(trapz(ftrans*fwave.value*(d(fwave.value)+numpy.random.normal(0,1.)*n(fwave.value)),fwave.value)/trapz(ftrans*fwave.value*v(fwave.value),fwave.value)))
            else:
                result.append(-2.5*numpy.log10(trapz(ftrans*(d(fwave.value)+numpy.random.normal(0,1.)*n(fwave.value)),fwave.value)/trapz(ftrans*v(fwave.value),fwave.value)))
        outunit = 1.

    elif (ab):
        nu = sp.wave.to('Hz',equivalencies=u.spectral())
        fnu = sp.flux.to('Jy',equivalencies=u.spectral_density(sp.wave))
        noisenu = sp.noise.to('Jy',equivalencies=u.spectral_density(sp.wave))
        filtnu = fwave.to('Hz',equivalencies=u.spectral())
        fconst = 3631*u.jansky
        d = interp1d(nu.value,fnu.value,bounds_error=False,fill_value=0.)
        n = interp1d(nu.value,noisenu.value,bounds_error=False,fill_value=0.)
        b = trapz((ftrans/filtnu.value)*fconst.value,filtnu.value)
        val = -2.5*numpy.log10(trapz(ftrans*d(filtnu.value)/filtnu.value,filtnu.value)/b)
        for i in numpy.arange(nsamples):
            a = trapz(ftrans*(d(filtnu.value)+numpy.random.normal(0,1)*n(filtnu.value))/filtnu.value,filtnu.value)
            result.append(-2.5*numpy.log10(a/b))
        outunit = 1.

    elif (energy):
        outunit = u.erg/u.s/u.cm**2
        if rsr:
            a = trapz(ftrans*fwave.value*d(fwave.value),fwave.value)*sp.wave.unit*sp.flux.unit
            b = trapz(ftrans*fwave.value,fwave.value)*sp.wave.unit
            c = trapz(ftrans*fwave.value*fwave.value,fwave.value)*sp.wave.unit*sp.wave.unit
            val = (a/b * c/b).to(outunit).value
        else:
            a = trapz(ftrans*d(fwave.value),fwave.value)*sp.wave.unit*sp.flux.unit
            b = trapz(ftrans,fwave.value)*sp.wave.unit
            c = trapz(ftrans*fwave.value,fwave.value)*sp.wave.unit*sp.wave.unit
            val = (a/b * c/b).to(outunit).value
        for i in numpy.arange(nsamples):
            if rsr:
                result.append((trapz(ftrans*fwave.value*(d(fwave.value)+numpy.random.normal(0,1.)*n(fwave.value)),fwave.value)*sp.wave.unit*sp.flux.unit).to(outunit).value)
            else:
                result.append((trapz(ftrans*(d(fwave.value)+numpy.random.normal(0,1.)*n(fwave.value)),fwave.value)*sp.wave.unit*sp.flux.unit).to(outunit).value)

    elif (photons):
        outunit = 1./u.s/u.cm**2
        convert = const.h.to('erg s')*const.c.to('micron/s')
        val = (trapz(ftrans*fwave.value*convert.value*d(fwave.value),fwave.value)*sp.wave.unit*sp.flux.unit*convert.unit).to(outunit).value
        for i in numpy.arange(nsamples):
            result.append((trapz(ftrans*fwave.value*convert.value*(d(fwave.value)+numpy.random.normal(0,1.)*n(fwave.value)),fwave.value)*sp.wave.unit*sp.flux.unit*convert.unit).to(outunit).value)
    else:
        raise NameError('\nfilterMag not given a correct physical quantity (vega, ab, energy, photons) to compute photometry\n\n')


#    val = numpy.nanmean(result)*outunit
    err = numpy.nanstd(result)
    if len(sp.wave[wgood]) == 0:
        err = 0.
    return val*outunit,err*outunit


def vegaToAB(filt,vegafile=VEGAFILE,filterfolder=SPLAT_PATH+FILTER_FOLDER,custom=False,notch=False,rsr=False,**kwargs):

# check that requested filter is in list
    if isinstance(custom,bool) and isinstance(notch,bool):
        f0 = checkFilterName(filt,verbose=True)
        if f0 == False: 
            return numpy.nan, numpy.nan
        filt = f0
        rsr = FILTERS[filt]['rsr']

# Read in filter
    if isinstance(custom,bool) and isinstance(notch,bool):
        fwave,ftrans = filterProfile(filt,**kwargs)
# notch filter
    elif isinstance(custom,bool) and isinstance(notch,list):
        dn = (notch[1]-notch[0])/1000
        fwave = numpy.arange(notch[0]-5.*dn,notch[1]+5.*dn,dn)
        ftrans = numpy.zeros(len(fwave))
        ftrans[numpy.where(numpy.logical_and(fwave >= notch[0],fwave <= notch[1]))] = 1.
# custom filter
    else:
        fwave,ftrans = custom[0],custom[1]


# Read in Vega spectrum
    vwave,vflux = numpy.genfromtxt(os.path.normpath(filterfolder+vegafile), comments='#', unpack=True, \
        missing_values = ('NaN','nan'), filling_values = (numpy.nan))
    vwave = vwave[~numpy.isnan(vflux)]*u.micron
    vflux = vflux[~numpy.isnan(vflux)]*(u.erg/(u.cm**2 * u.s * u.micron))

# trim spectrum
    vflux = vflux[vwave>=numpy.nanmin(fwave)]
    vwave = vwave[vwave>=numpy.nanmin(fwave)]
    vflux = vflux[vwave<=numpy.nanmax(fwave)]
    vwave = vwave[vwave<=numpy.nanmax(fwave)]

# convert to fnu
    nu = vwave.to('Hz',equivalencies=u.spectral())
    fnu = vflux.to('Jy',equivalencies=u.spectral_density(vwave))
    filtnu = fwave.to('Hz',equivalencies=u.spectral())
    fconst = 3631*u.jansky
    d = interp1d(nu.value,fnu.value,bounds_error=False,fill_value=0.)
    b = trapz((ftrans/filtnu.value)*fconst.value,filtnu.value)
    return -2.5*numpy.log10(trapz(ftrans*d(filtnu.value)/filtnu.value,filtnu.value)/b)



def filterInfo(*args,**kwargs):
    '''
    :Purpose: Prints out the current list of filters in the SPLAT reference library.
    '''

    verbose = kwargs.get('verbose',True)

    if len(args) > 0: 
        fname = list(args)
    elif kwargs.get('filter',False) != False: 
        fname = kwargs['filter']
    else: 
        fname = sorted(list(FILTERS.keys()))
    if isinstance(fname,list) == False: 
        fname = [fname]

    output = {}
    for k in fname:
        f = checkFilterName(k)
        if f != False:
            output[f] = {}
            output[f]['description'] = FILTERS[f]['description']
            output[f]['zeropoint'] = FILTERS[f]['zeropoint']
            fwave,ftrans = filterProfile(f,**kwargs)
            try:
                fwave = fwave.to(u.micron)
            except:
                fwave = fwave*u.micron
            fw = fwave[numpy.where(ftrans > 0.01*numpy.nanmax(ftrans))]
            ft = ftrans[numpy.where(ftrans > 0.01*numpy.nanmax(ftrans))]
            fw05 = fwave[numpy.where(ftrans > 0.5*numpy.nanmax(ftrans))]
            output[f]['lambda_mean'] = trapz(ft*fw,fw)/trapz(ft,fw)
            output[f]['lambda_pivot'] = numpy.sqrt(trapz(fw*ft,fw)/trapz(ft/fw,fw))
            output[f]['lambda_central'] = 0.5*(numpy.max(fw)+numpy.min(fw))
            output[f]['lambda_fwhm'] = numpy.max(fw05)-numpy.min(fw05)
            output[f]['lambda_min'] = numpy.min(fw)
            output[f]['lambda_max'] = numpy.max(fw)
            if verbose ==True: 
                print(f.replace('_',' ')+': '+output[f]['zeropoint'])
                print('Zeropoint = {} Jy'.format(output[f]['zeropoint']))
                print('Central wavelength: = {:.3f}'.format(output[f]['lambda_central']))
                print('Mean wavelength: = {:.3f}'.format(output[f]['lambda_mean']))
                print('Pivot point: = {:.3f}'.format(output[f]['lambda_pivot']))
                print('FWHM = {:.3f}'.format(output[f]['lambda_fwhm']))
                print('Wavelength range = {:.3f} to {:.3f}\n'.format(output[f]['lambda_min'],output[f]['lambda_max']))             
        else:
        	if verbose ==True: print('  Filter {} not in SPLAT filter list'.format(k))
    kys = list(output.keys())
    if len(kys) == 1: return output[kys[0]]
    else: return output


def filterProperties(filt,**kwargs):
    '''
    :Purpose: Returns a dictionary containing key parameters for a particular filter.

    :param filter: name of filter, must be one of the specifed filters given by splat.FILTERS.keys()
    :type filter: required
    :param verbose: print out information about filter to screen
    :type verbose: optional, default = True

    :Example:
    >>> import splat
    >>> data = splat.filterProperties('2MASS J')
    Filter 2MASS J: 2MASS J-band
    Zeropoint = 1594.0 Jy
    Pivot point: = 1.252 micron
    FWHM = 0.323 micron
    Wavelength range = 1.066 to 1.442 micron
    >>> data = splat.filterProperties('2MASS X')
    Filter 2MASS X not among the available filters:
      2MASS H: 2MASS H-band
      2MASS J: 2MASS J-band
      2MASS KS: 2MASS Ks-band
      BESSEL I: Bessel I-band
      FOURSTAR H: FOURSTAR H-band
      FOURSTAR H LONG: FOURSTAR H long
      FOURSTAR H SHORT: FOURSTAR H short
      ...
    '''
    filterFolder = kwargs.get('filterFolder',SPLAT_PATH+FILTER_FOLDER)
    if not os.path.exists(filterFolder):
        filterFolder = SPLAT_URL+FILTER_FOLDER

# check that requested filter is in list
    filt = checkFilterName(filt)
    if filt == False: return None

    report = {}
    report['name'] = filt
    report['description'] = FILTERS[filt]['description']
    report['zeropoint'] = FILTERS[filt]['zeropoint']
    report['method'] = FILTERS[filt]['method']
    report['rsr'] = FILTERS[filt]['rsr']
    fwave,ftrans = filterProfile(filt,**kwargs)
    try:
        fwave = fwave.to(u.micron)
    except:
        fwave = fwave*u.micron
    fw = fwave[numpy.where(ftrans > 0.01*numpy.nanmax(ftrans))]
    ft = ftrans[numpy.where(ftrans > 0.01*numpy.nanmax(ftrans))]
    fw05 = fwave[numpy.where(ftrans > 0.5*numpy.nanmax(ftrans))]
#        print(trapz(ft,fw))
#        print(trapz(fw*ft,fw))
    report['lambda_mean'] = trapz(ft*fw,fw)/trapz(ft,fw)
    report['lambda_pivot'] = numpy.sqrt(trapz(fw*ft,fw)/trapz(ft/fw,fw))
    report['lambda_central'] = 0.5*(numpy.max(fw)+numpy.min(fw))
    report['lambda_fwhm'] = numpy.max(fw05)-numpy.min(fw05)
    report['lambda_min'] = numpy.min(fw)
    report['lambda_max'] = numpy.max(fw)
    report['wave'] = fwave
    report['transmission'] = ftrans
# report values out
    if kwargs.get('verbose',False):
        print('\nFilter '+filt+': '+report['description'])
        print('Zeropoint = {} Jy'.format(report['zeropoint']))
        print('Pivot point: = {:.3f}'.format(report['lambda_pivot']))
        print('FWHM = {:.3f}'.format(report['lambda_fwhm']))
        print('Wavelength range = {:.3f} to {:.3f}\n'.format(report['lambda_min'],report['lambda_max']))
    return report


def magToFlux(mag,filt,**kwargs):
    '''
    :Purpose: Converts a magnitude into an energy, and vice versa.

    :param mag: magnitude on whatever system is defined for the filter or provided (required)
    :param filter: name of filter, must be one of the specifed filters given by splat.FILTERS.keys() (required)
    :param reverse: convert energy into magnitude instead (optional, default = False)
    :param ab: magnitude is on the AB system (optional, default = filter preference)
    :param vega: magnitude is on the Vega system (optional, default = filter preference)
    :param rsr: magnitude is on the Vega system (optional, default = filter preference)
    :param units: units for energy as an astropy.units variable; if this conversion does not work, the conversion is ignored (optional, default = erg/cm2/s)
    :param verbose: print out information about filter to screen (optional, default = False)

    WARNING: THIS CODE IS ONLY PARTIALLY COMPLETE
    '''

# keyword parameters
    filterFolder = kwargs.get('filterFolder',SPLAT_PATH+FILTER_FOLDER)
    if not os.path.exists(filterFolder):
        filterFolder = SPLAT_URL+FILTER_FOLDER
    vegaFile = kwargs.get('vegaFile','vega_kurucz.txt')
    vega = kwargs.get('vega',True)
    ab = kwargs.get('ab',not vega)
    rsr = kwargs.get('rsr',False)
    nsamples = kwargs.get('nsamples',100)
    custom = kwargs.get('custom',False)
    notch = kwargs.get('notch',False)
    base_unit = u.erg/(u.cm**2 * u.s)
    return_unit = kwargs.get('unit',base_unit)
    e_mag = kwargs.get('uncertainty',0.)
    e_mag = kwargs.get('unc',e_mag)
    e_mag = kwargs.get('e_mag',e_mag)
    if not isinstance(mag,u.quantity.Quantity): mag=mag*u.s/u.s
    if not isinstance(e_mag,u.quantity.Quantity): e_mag=e_mag*mag.unit

# check that requested filter is in list
    filt = checkFilterName(filt)
    if filt == False: return numpy.nan, numpy.nan

# reset filter calculation methods based on filter design
    if 'ab' in FILTERS[filt]['method']: 
        ab = kwargs.get('ab',True)
        vega = not ab
    if 'vega' in FILTERS[filt]['method']: 
        vega = kwargs.get('vega',True)
        ab = not vega
    if 'rsr' in FILTERS[filt]['method']: 
        rsr = kwargs.get('rsr',True)


# Read in filter
    if isinstance(custom,bool) and isinstance(notch,bool):
        fwave,ftrans = filterProfile(filt,**kwargs)
# notch filter
    elif isinstance(custom,bool) and isinstance(notch,list):
        dn = (notch[1]-notch[0])/1000
        fwave = numpy.arange(notch[0]-5.*dn,notch[1]+5.*dn,dn)*u.micron
        ftrans = numpy.zeros(len(fwave))
        ftrans[numpy.where(numpy.logical_and(fwave >= notch[0],fwave <= notch[1]))] = 1.
# custom filter
    else:
        fwave,ftrans = custom[0],custom[1]
    if isinstance(fwave,u.quantity.Quantity) == False: fwave=fwave*u.micron
    if isinstance(ftrans,u.quantity.Quantity) == True: ftrans=ftrans.value
    fwave = fwave[~numpy.isnan(ftrans)]
    ftrans = ftrans[~numpy.isnan(ftrans)]

    result = []
    err = 0.
# magnitude -> energy
    if kwargs.get('reverse',False) == False:
        
        if vega == True:
    # Read in Vega spectrum
            vwave,vflux = numpy.genfromtxt(os.path.normpath(filterFolder+vegaFile), comments='#', unpack=True, \
                missing_values = ('NaN','nan'), filling_values = (numpy.nan))
            vwave = vwave[~numpy.isnan(vflux)]*u.micron
            vflux = vflux[~numpy.isnan(vflux)]*(u.erg/(u.cm**2 * u.s * u.micron))
    # interpolate Vega onto filter wavelength function
            v = interp1d(vwave.value,vflux.value,bounds_error=False,fill_value=0.)
            if rsr: fact = trapz(ftrans*fwave.value*v(fwave.value),fwave.value)
            else: fact = trapz(ftrans*v(fwave.value),fwave.value)
            val = 10.**(-0.4*mag.value)*fact*u.erg/(u.cm**2 * u.s)
    # calculate uncertainty        
            if e_mag.value > 0.:
                for i in numpy.arange(nsamples): result.append(10.**(-0.4*(mag.value+numpy.random.normal(0,1.)*e_mag.value))*fact)
                err = (numpy.nanstd(result))*u.erg/(u.cm**2 * u.s)
            else: err = 0.*u.erg/(u.cm**2 * u.s)
        elif ab == True:
            fconst = 3631*u.jansky
            ftrans = (ftrans*fconst).to(u.erg/(u.cm**2 * u.s * u.micron),equivalencies=u.spectral_density(fwave))
            if rsr: fact = trapz(ftrans.value*fwave.value,fwave.value)
            else: fact = trapz(ftrans.value,fwave.value)
            val = (10.**(-0.4*mag.value)*fact)*u.erg/(u.cm**2 * u.s)
    # calculate uncertainty        
            if e_mag.value > 0.:
                for i in numpy.arange(nsamples): result.append(10.**(-0.4*(mag.value+numpy.random.normal(0,1.)*e_mag.value))*fact)
                err = (numpy.nanstd(result))*u.erg/(u.cm**2 * u.s)
            else: err = 0.*u.erg/(u.cm**2 * u.s)
        else:
            raise ValueError('\nmagToFlux needs vega or ab method specified')

# convert to desired energy units
#        try:
        val.to(return_unit)
        err.to(return_unit)
#        except:
#            print('\nWarning: unit {} is not an energy flux unit'.format(return_unit))
        try:
            val.to(base_unit)
            err.to(base_unit)
        except:
            print('\nWarning: cannot convert result to an energy flux unit'.format(base_unit))
            return numpy.nan, numpy.nan
        return val, err

# energy -> magnitude
# THIS NEEDS TO BE COMPLETED
    else:
        print('passed')
        pass
# check that input is an energy flux
#        try:
#            mag.to(base_unit)        
#            e_mag.to(base_unit)        
#        except:
#            raise ValueError('\nInput quantity unit {} is not a flux unit'.format(mag.unit))


def visualizeFilter(filters,verbose=True,xra=[],yra=[0,1.2],**kwargs):
    '''
    :Purpose: Plots a filter profile or set of filter profiles, optionally on top of a spectrum

    WARNING: THIS CODE IS CURRENTLY UNDER DEVELOPMENT, BUGS MAY BE COMMON

    '''
    filt = copy.deepcopy(filters)
    wave_unit = kwargs.get('wave_unit',DEFAULT_WAVE_UNIT)

# single filter name  
    if isinstance(filt,str):
        filt = [filt]

    if isinstance(filt,list):

# list of filter names
        if isinstance(filt[0],str):
            for f in filt:
                fc = checkFilterName(f)
                filt.remove(f)
                if fc == False: 
                    if verbose==True: print('Removed filter {}: not included in SPLAT'.format(f))
                else:
                    filt.insert(len(filt),fc)
            if len(filt) == 0:
                raise ValueError('Did not recognize any of the input filters {}'.format(filters))

# prep parameters
            fwave,ftrans = filterProfile(f,**kwargs)
            if isUnit(fwave): wave_unit = kwargs.get('wave_unit',fwave.unit)

            xl = kwargs.get('xlabel','Wavelength ({})'.format(wave_unit))
            yl = kwargs.get('ylabel','Transmission Curve')
            legend = []
            fig = plt.figure(figsize=kwargs.get('figsize',[5,4]))
            for i,f in enumerate(filt):
                fwave,ftrans = filterProfile(f,**kwargs)
                if isUnit(fwave): fwave.to(wave_unit)
                else: fwave = fwave*wave_unit
                if kwargs.get('normalize',False): ftrans = ftrans/numpy.nanmax(ftrans)
                plt.plot(fwave,ftrans)
                if len(xra) == 0: xra = [numpy.nanmin(fwave.value),numpy.nanmax(fwave.value)]
                xra = [numpy.nanmin([xra[0],numpy.nanmin(fwave.value)]),numpy.nanmax([xra[1],numpy.nanmax(fwave.value)])]
                yra = [yra[0],numpy.nanmax([yra[1],numpy.nanmax(ftrans)])]
                legend.append(FILTERS[f]['description'])
                if FILTERS[f]['rsr'] == True: yl = kwargs.get('ylabel','Transmission Curve')

# list of notch ranges
        if isinstance(filt[0],int) or isinstance(filt[0],float):
            filt = [filt]

# list of notch ranges
        if isinstance(filt[0],list):
            xl = kwargs.get('xlabel','Wavelength ({})'.format(wave_unit))
            yl = kwargs.get('ylabel','Transmission Curve')
            legend = []
            fig = plt.figure(figsize=kwargs.get('figsize',[5,4]))
            for i,f in enumerate(filt):
                fwave,ftrans = numpy.linspace(f[0],f[1],1000)*wave_unit,numpy.ones(1000)
                plt.plot(fwave,ftrans)
                if len(xra) == 0: xra = [numpy.nanmin(fwave.value),numpy.nanmax(fwave.value)]
                xra = [numpy.nanmin([xra[0],numpy.nanmin(fwave.value)]),numpy.nanmax([xra[1],numpy.nanmax(fwave.value)])]
                yra = [yra[0],numpy.nanmax([yra[1],numpy.nanmax(ftrans)])]
                legend.append('Filter {}'.format(i+1))

    else:
        raise ValueError('Could not parse input {}'.format(filt))
        

# add a comparison spectrum
    sp = kwargs.get('spectrum',None)
    sp = kwargs.get('comparison',sp)

    if isinstance(sp,splat.core.Spectrum) == True:
        print(xra)
        sp.normalize(xra)
        sp.scale(numpy.nanmax(ftrans)*kwargs.get('comparison_scale',0.8))
        plt.plot(sp.wave,sp.flux,color=kwargs.get('comparison_color','k'),alpha=kwargs.get('comparison_alpha',0.5))
        legend.append(sp.name)
        yra = [yra[0],yra[1]*1.1]

# finish up
    plt.xlim(xra)
    plt.ylim(yra)
    plt.xlabel(xl)
    plt.ylabel(yl)
    plt.legend(legend)

# save if desired
    file = kwargs.get('file','')
    file = kwargs.get('filename',file)
    file = kwargs.get('output',file)
    if file != '': plt.savefig(file)
    
    return fig


#########################################
########   SED FITTING TOOLS    #########
### WARNING: THESE ARE EXPERIMENTAL!! ###
#########################################

# plan:

def modelMagnitudes(verbose=True):
    '''
    this will be a code that calculates a set of magnitudes for a model set's SED models
    saves to file that could be uploaded
    pre-save some model magnitudes
    '''
    pass


def interpolateMagnitudes(verbose=True):
    '''
    produces an interpolated value for a grid set of model magnitudes
    '''
    pass

def compareMagnitudes(mags1,mags2,unc=None,unc2=None,ignore=[],verbose=True):
    '''
    this code compares a set of magnitudes using one of several statistics
    '''
    chi = 0.
    dm,em = [],[]
    for f in list(mags1.keys()):
        if f in list(mags2.keys()) and f in list(unc.keys()) and f not in ignore: 
            dm.append(mags1[f]-mags2[f])
            em.append(unc[f])
# find best scale factor
    dm = numpy.array(dm)
    em = numpy.array(em)
    offset = numpy.sum(dm/em**2)/numpy.sum (1./em**2)
    dmo = numpy.array([m-offset for m in dm])
    return numpy.sum((dmo/em)**2), offset

def SEDFitGrid(verbose=True):
    '''
    this code will compare a set of magnitudes to a grid of model magnitudes and choose the
    closest match based on various statistics
    '''
    pass

def SEDFitMCMC(verbose=True):
    '''
    this code will conduct a comparison of a set of magnitudes to model magnitudes using an
    MCMC wrapper, and choose best/average/distribution of parameters
    '''
    pass

def SEDFitAmoeba(verbose=True):
    '''
    this code will conduct a comparison of a set of magnitudes to model magnitudes using an
    Amoeba (Nelder-Mead) wrapper, and choose the closest match
    '''
    pass

def SEDVisualize(verbose=True):
    '''
    Visualizes magnitudes on SED scale (flux = lam x F_lam), with option of also comparing to SED spectrum
    '''
    pass



#####################################################
###############    MAGNITUDE CLASS    ###############
#####################################################

class Magnitude(object):
    '''
    :Description: 

        This is a class data structure for a magnitude value

    '''

    def __init__(self, magnitude, filt, uncertainty=0., magtype='apparent', verbose=False,**kwargs):
        self.magnitude = magnitude
        self.uncertainty = uncertainty
        self.type = magtype

# check filter and rename if necessary
        self.knownfilter = True
        fflag = checkFilterName(filt,verbose=verbose)
        if fflag == False:
            if verbose== True: print('filter {} is not a standard filter; some functions may not work'.format(filt))
            self.knownfilter = False
        else: filt = fflag
        self.filter = filt

# some things that are based on presets
        if self.knownfilter == True:
            self.wave,self.transmission = filterProfile(self.filter)
            info = filterProperties(self.filter)
            for k in info.keys(): setattr(self,k,info[k])

    def __copy__(self):
        '''
        :Purpose: Make a copy of a Magnitude object
        '''
        s = type(self)(self.magnitude,self.filter,uncertainty=self.uncertainty)
        s.__dict__.update(self.__dict__)
        return s

# backup version
    def copy(self):
        '''
        :Purpose: Make a copy of a Magnitude object
        '''
        s = type(self)(self.magnitude,self.filter,uncertainty=self.uncertainty)
        s.__dict__.update(self.__dict__)
        return s

    def __repr__(self):
        '''
        :Purpose: A simple representation of the Spectrum object
        '''
        if self.uncertainty != 0. and numpy.isfinite(self.uncertainty):
            return '{} magnitude of {}+/-{}'.format(self.filter,self.magnitude,self.uncertainty)
        else:
            return '{} magnitude of {}'.format(self.filter,self.magnitude)

    def __add__(self,other,samp=1000):
        '''
        :Purpose: 

            A representation of addition for Magnitude classes that takes into account uncertainties

        :Output: 

            A new Magnitude object equal to the sum of values
        '''

# make a copy and fill in combined magnitude
        out = copy.deepcopy(self)
        out.magnitude = self.magnitude+other.magnitude
        out.uncertainty = self.uncertainty+other.uncertainty

# combine noises
        if self.uncertainty != 0 and other.uncertainty != 0:
            m1 = numpy.random.normal(self.magnitude,self.uncertainty,samp)
            m2 = numpy.random.normal(other.magnitude,other.uncertainty,samp)
            val = m1+m2
            out.uncertainty = numpy.nanstd(val)

# check filter agreement
        if self.filter != other.filter:
            out.filter = '{}+{}'.format(self.filter,other.filter)

        return out

    def __sub__(self,other,samp=1000):
        '''
        :Purpose: 

            A representation of subtraction for Magnitude classes that takes into account uncertainties

        :Output: 

            A new Magnitude object equal to the diffence of values
        '''

# make a copy and fill in combined magnitude
        out = copy.deepcopy(self)
        out.magnitude = self.magnitude-other.magnitude
        out.uncertainty = self.uncertainty+other.uncertainty

# combine noises
        if self.uncertainty != 0 and other.uncertainty != 0:
            m1 = numpy.random.normal(self.magnitude,self.uncertainty,samp)
            m2 = numpy.random.normal(other.magnitude,other.uncertainty,samp)
            val = m1-m2
            out.uncertainty = numpy.nanstd(val)

# check filter agreement
        if self.filter != other.filter:
            out.filter = '{}-{}'.format(self.filter,other.filter)

        return out

    def flux(self,type='fnu',samp=1000):
        '''
        :Purpose: 

            Report the equivalent flux density of a magnitude

        :Output: 

            astropy quantity in flux density units (default = erg/cm2/s/micron)

        '''
        pass


    def addFlux(self,other,samp=1000):
        '''
        :Purpose: 

            A representation of addition for magnitudes (addition of fluxes)

        :Output: 

            A new magnitude object equal to the equivalent sum of fluxes

        '''
# check filter agreement
        if self.filter != other.filter:
            raise ValueError('magnitudes filters {} and {} are not the same'.format(self.filter,other.filter))

# make a copy and fill in combined magnitude
        out = copy.deepcopy(self)
        out.magnitude = self.magnitude-2.5*numpy.log10(1.+10.**(-0.4*(other.magnitude-self.magnitude)))
        out.uncertainty = self.uncertainty+other.uncertainty

# combine noises
        if self.uncertainty != 0 and other.uncertainty != 0:
            m1 = numpy.random.normal(self.magnitude,self.uncertainty,samp)
            m2 = numpy.random.normal(other.magnitude,other.uncertainty,samp)
            val = m1-2.5*numpy.log10(1.+10.**(-0.4*(m2-m1)))
            out.uncertainty = numpy.nanstd(val)

        return out



[((895, 914), 'copy.deepcopy', 'copy.deepcopy', (['filt'], {}), '(filt)\n', (908, 914), False, 'import copy\n'), ((12197, 12217), 'numpy.nanstd', 'numpy.nanstd', (['result'], {}), '(result)\n', (12209, 12217), False, 'import numpy\n'), ((13945, 14010), 'scipy.interpolate.interp1d', 'interp1d', (['nu.value', 'fnu.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(nu.value, fnu.value, bounds_error=False, fill_value=0.0)\n', (13953, 14010), False, 'from scipy.interpolate import interp1d\n'), ((14015, 14072), 'scipy.integrate.trapz', 'trapz', (['(ftrans / filtnu.value * fconst.value)', 'filtnu.value'], {}), '(ftrans / filtnu.value * fconst.value, filtnu.value)\n', (14020, 14072), False, 'from scipy.integrate import trapz\n'), ((18601, 18614), 'numpy.min', 'numpy.min', (['fw'], {}), '(fw)\n', (18610, 18614), False, 'import numpy\n'), ((18642, 18655), 'numpy.max', 'numpy.max', (['fw'], {}), '(fw)\n', (18651, 18655), False, 'import numpy\n'), ((25140, 25162), 'copy.deepcopy', 'copy.deepcopy', (['filters'], {}), '(filters)\n', (25153, 25162), False, 'import copy\n'), ((28421, 28434), 'matplotlib.pyplot.xlim', 'plt.xlim', (['xra'], {}), '(xra)\n', (28429, 28434), True, 'import matplotlib.pyplot as plt\n'), ((28439, 28452), 'matplotlib.pyplot.ylim', 'plt.ylim', (['yra'], {}), '(yra)\n', (28447, 28452), True, 'import matplotlib.pyplot as plt\n'), ((28457, 28471), 'matplotlib.pyplot.xlabel', 'plt.xlabel', (['xl'], {}), '(xl)\n', (28467, 28471), True, 'import matplotlib.pyplot as plt\n'), ((28476, 28490), 'matplotlib.pyplot.ylabel', 'plt.ylabel', (['yl'], {}), '(yl)\n', (28486, 28490), True, 'import matplotlib.pyplot as plt\n'), ((28495, 28513), 'matplotlib.pyplot.legend', 'plt.legend', (['legend'], {}), '(legend)\n', (28505, 28513), True, 'import matplotlib.pyplot as plt\n'), ((29664, 29679), 'numpy.array', 'numpy.array', (['dm'], {}), '(dm)\n', (29675, 29679), False, 'import numpy\n'), ((29689, 29704), 'numpy.array', 'numpy.array', (['em'], {}), '(em)\n', (29700, 29704), False, 'import numpy\n'), ((29769, 29808), 'numpy.array', 'numpy.array', (['[(m - offset) for m in dm]'], {}), '([(m - offset) for m in dm])\n', (29780, 29808), False, 'import numpy\n'), ((2038, 2066), 'os.path.exists', 'os.path.exists', (['filterFolder'], {}), '(filterFolder)\n', (2052, 2066), False, 'import os\n'), ((2307, 2361), 'os.path.normpath', 'os.path.normpath', (["(filterFolder + FILTERS[filt]['file'])"], {}), "(filterFolder + FILTERS[filt]['file'])\n", (2323, 2361), False, 'import os\n'), ((4966, 4994), 'os.path.exists', 'os.path.exists', (['filterFolder'], {}), '(filterFolder)\n', (4980, 4994), False, 'import os\n'), ((7605, 7697), 'scipy.interpolate.interp1d', 'interp1d', (['sp.wave[wgood].value', 'sp.flux[wgood].value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(sp.wave[wgood].value, sp.flux[wgood].value, bounds_error=False,\n    fill_value=0.0)\n', (7613, 7697), False, 'from scipy.interpolate import interp1d\n'), ((7702, 7793), 'scipy.interpolate.interp1d', 'interp1d', (['sp.wave[wgood].value', 'sp.noise[wgood].value'], {'bounds_error': '(False)', 'fill_value': '(0)'}), '(sp.wave[wgood].value, sp.noise[wgood].value, bounds_error=False,\n    fill_value=0)\n', (7710, 7793), False, 'from scipy.interpolate import interp1d\n'), ((7943, 8017), 'scipy.interpolate.interp1d', 'interp1d', (['sp.wave.value', 'sp.flux.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(sp.wave.value, sp.flux.value, bounds_error=False, fill_value=0.0)\n', (7951, 8017), False, 'from scipy.interpolate import interp1d\n'), ((8026, 8112), 'scipy.interpolate.interp1d', 'interp1d', (['sp.wave.value', '(sp.flux.value * 1e-09)'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(sp.wave.value, sp.flux.value * 1e-09, bounds_error=False,\n    fill_value=0.0)\n', (8034, 8112), False, 'from scipy.interpolate import interp1d\n'), ((8606, 8676), 'scipy.interpolate.interp1d', 'interp1d', (['vwave.value', 'vflux.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(vwave.value, vflux.value, bounds_error=False, fill_value=0.0)\n', (8614, 8676), False, 'from scipy.interpolate import interp1d\n'), ((8980, 9002), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (8992, 9002), False, 'import numpy\n'), ((13258, 13299), 'os.path.normpath', 'os.path.normpath', (['(filterfolder + vegafile)'], {}), '(filterfolder + vegafile)\n', (13274, 13299), False, 'import os\n'), ((17502, 17530), 'os.path.exists', 'os.path.exists', (['filterFolder'], {}), '(filterFolder)\n', (17516, 17530), False, 'import os\n'), ((18347, 18365), 'scipy.integrate.trapz', 'trapz', (['(ft * fw)', 'fw'], {}), '(ft * fw, fw)\n', (18352, 18365), False, 'from scipy.integrate import trapz\n'), ((18363, 18376), 'scipy.integrate.trapz', 'trapz', (['ft', 'fw'], {}), '(ft, fw)\n', (18368, 18376), False, 'from scipy.integrate import trapz\n'), ((18542, 18557), 'numpy.max', 'numpy.max', (['fw05'], {}), '(fw05)\n', (18551, 18557), False, 'import numpy\n'), ((18558, 18573), 'numpy.min', 'numpy.min', (['fw05'], {}), '(fw05)\n', (18567, 18573), False, 'import numpy\n'), ((20230, 20258), 'os.path.exists', 'os.path.exists', (['filterFolder'], {}), '(filterFolder)\n', (20244, 20258), False, 'import os\n'), ((28661, 28678), 'matplotlib.pyplot.savefig', 'plt.savefig', (['file'], {}), '(file)\n', (28672, 28678), True, 'import matplotlib.pyplot as plt\n'), ((29718, 29741), 'numpy.sum', 'numpy.sum', (['(dm / em ** 2)'], {}), '(dm / em ** 2)\n', (29727, 29741), False, 'import numpy\n'), ((29738, 29762), 'numpy.sum', 'numpy.sum', (['(1.0 / em ** 2)'], {}), '(1.0 / em ** 2)\n', (29747, 29762), False, 'import numpy\n'), ((29816, 29842), 'numpy.sum', 'numpy.sum', (['((dmo / em) ** 2)'], {}), '((dmo / em) ** 2)\n', (29825, 29842), False, 'import numpy\n'), ((32934, 32953), 'copy.deepcopy', 'copy.deepcopy', (['self'], {}), '(self)\n', (32947, 32953), False, 'import copy\n'), ((33860, 33879), 'copy.deepcopy', 'copy.deepcopy', (['self'], {}), '(self)\n', (33873, 33879), False, 'import copy\n'), ((35205, 35224), 'copy.deepcopy', 'copy.deepcopy', (['self'], {}), '(self)\n', (35218, 35224), False, 'import copy\n'), ((6530, 6588), 'numpy.arange', 'numpy.arange', (['(notch[0] - 5.0 * dn)', '(notch[1] + 5.0 * dn)', 'dn'], {}), '(notch[0] - 5.0 * dn, notch[1] + 5.0 * dn, dn)\n', (6542, 6588), False, 'import numpy\n'), ((6989, 7008), 'numpy.nanmax', 'numpy.nanmax', (['fwave'], {}), '(fwave)\n', (7001, 7008), False, 'import numpy\n'), ((7011, 7032), 'numpy.nanmin', 'numpy.nanmin', (['sp.wave'], {}), '(sp.wave)\n', (7023, 7032), False, 'import numpy\n'), ((7036, 7055), 'numpy.nanmin', 'numpy.nanmin', (['fwave'], {}), '(fwave)\n', (7048, 7055), False, 'import numpy\n'), ((7058, 7079), 'numpy.nanmax', 'numpy.nanmax', (['sp.wave'], {}), '(sp.wave)\n', (7070, 7079), False, 'import numpy\n'), ((7241, 7260), 'numpy.nanmin', 'numpy.nanmin', (['fwave'], {}), '(fwave)\n', (7253, 7260), False, 'import numpy\n'), ((7263, 7284), 'numpy.nanmin', 'numpy.nanmin', (['sp.wave'], {}), '(sp.wave)\n', (7275, 7284), False, 'import numpy\n'), ((7288, 7307), 'numpy.nanmax', 'numpy.nanmax', (['fwave'], {}), '(fwave)\n', (7300, 7307), False, 'import numpy\n'), ((7310, 7331), 'numpy.nanmax', 'numpy.nanmax', (['sp.wave'], {}), '(sp.wave)\n', (7322, 7331), False, 'import numpy\n'), ((7538, 7559), 'numpy.isnan', 'numpy.isnan', (['sp.noise'], {}), '(sp.noise)\n', (7549, 7559), False, 'import numpy\n'), ((8198, 8239), 'os.path.normpath', 'os.path.normpath', (['(filterFolder + vegaFile)'], {}), '(filterFolder + vegaFile)\n', (8214, 8239), False, 'import os\n'), ((9921, 9986), 'scipy.interpolate.interp1d', 'interp1d', (['nu.value', 'fnu.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(nu.value, fnu.value, bounds_error=False, fill_value=0.0)\n', (9929, 9986), False, 'from scipy.interpolate import interp1d\n'), ((9995, 10064), 'scipy.interpolate.interp1d', 'interp1d', (['nu.value', 'noisenu.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(nu.value, noisenu.value, bounds_error=False, fill_value=0.0)\n', (10003, 10064), False, 'from scipy.interpolate import interp1d\n'), ((10073, 10130), 'scipy.integrate.trapz', 'trapz', (['(ftrans / filtnu.value * fconst.value)', 'filtnu.value'], {}), '(ftrans / filtnu.value * fconst.value, filtnu.value)\n', (10078, 10130), False, 'from scipy.integrate import trapz\n'), ((10235, 10257), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (10247, 10257), False, 'import numpy\n'), ((12951, 13009), 'numpy.arange', 'numpy.arange', (['(notch[0] - 5.0 * dn)', '(notch[1] + 5.0 * dn)', 'dn'], {}), '(notch[0] - 5.0 * dn, notch[1] + 5.0 * dn, dn)\n', (12963, 13009), False, 'import numpy\n'), ((13776, 13788), 'astropy.units.spectral', 'u.spectral', ([], {}), '()\n', (13786, 13788), True, 'from astropy import units as u\n'), ((13828, 13853), 'astropy.units.spectral_density', 'u.spectral_density', (['vwave'], {}), '(vwave)\n', (13846, 13853), True, 'from astropy import units as u\n'), ((13896, 13908), 'astropy.units.spectral', 'u.spectral', ([], {}), '()\n', (13906, 13908), True, 'from astropy import units as u\n'), ((15558, 15571), 'numpy.min', 'numpy.min', (['fw'], {}), '(fw)\n', (15567, 15571), False, 'import numpy\n'), ((15610, 15623), 'numpy.max', 'numpy.max', (['fw'], {}), '(fw)\n', (15619, 15623), False, 'import numpy\n'), ((18416, 18434), 'scipy.integrate.trapz', 'trapz', (['(fw * ft)', 'fw'], {}), '(fw * ft, fw)\n', (18421, 18434), False, 'from scipy.integrate import trapz\n'), ((18432, 18450), 'scipy.integrate.trapz', 'trapz', (['(ft / fw)', 'fw'], {}), '(ft / fw, fw)\n', (18437, 18450), False, 'from scipy.integrate import trapz\n'), ((18485, 18498), 'numpy.max', 'numpy.max', (['fw'], {}), '(fw)\n', (18494, 18498), False, 'import numpy\n'), ((18499, 18512), 'numpy.min', 'numpy.min', (['fw'], {}), '(fw)\n', (18508, 18512), False, 'import numpy\n'), ((22694, 22764), 'scipy.interpolate.interp1d', 'interp1d', (['vwave.value', 'vflux.value'], {'bounds_error': '(False)', 'fill_value': '(0.0)'}), '(vwave.value, vflux.value, bounds_error=False, fill_value=0.0)\n', (22702, 22764), False, 'from scipy.interpolate import interp1d\n'), ((32387, 32419), 'numpy.isfinite', 'numpy.isfinite', (['self.uncertainty'], {}), '(self.uncertainty)\n', (32401, 32419), False, 'import numpy\n'), ((33166, 33225), 'numpy.random.normal', 'numpy.random.normal', (['self.magnitude', 'self.uncertainty', 'samp'], {}), '(self.magnitude, self.uncertainty, samp)\n', (33185, 33225), False, 'import numpy\n'), ((33241, 33302), 'numpy.random.normal', 'numpy.random.normal', (['other.magnitude', 'other.uncertainty', 'samp'], {}), '(other.magnitude, other.uncertainty, samp)\n', (33260, 33302), False, 'import numpy\n'), ((33355, 33372), 'numpy.nanstd', 'numpy.nanstd', (['val'], {}), '(val)\n', (33367, 33372), False, 'import numpy\n'), ((34092, 34151), 'numpy.random.normal', 'numpy.random.normal', (['self.magnitude', 'self.uncertainty', 'samp'], {}), '(self.magnitude, self.uncertainty, samp)\n', (34111, 34151), False, 'import numpy\n'), ((34167, 34228), 'numpy.random.normal', 'numpy.random.normal', (['other.magnitude', 'other.uncertainty', 'samp'], {}), '(other.magnitude, other.uncertainty, samp)\n', (34186, 34228), False, 'import numpy\n'), ((34281, 34298), 'numpy.nanstd', 'numpy.nanstd', (['val'], {}), '(val)\n', (34293, 34298), False, 'import numpy\n'), ((35486, 35545), 'numpy.random.normal', 'numpy.random.normal', (['self.magnitude', 'self.uncertainty', 'samp'], {}), '(self.magnitude, self.uncertainty, samp)\n', (35505, 35545), False, 'import numpy\n'), ((35561, 35622), 'numpy.random.normal', 'numpy.random.normal', (['other.magnitude', 'other.uncertainty', 'samp'], {}), '(other.magnitude, other.uncertainty, samp)\n', (35580, 35622), False, 'import numpy\n'), ((35712, 35729), 'numpy.nanstd', 'numpy.nanstd', (['val'], {}), '(val)\n', (35724, 35729), False, 'import numpy\n'), ((2877, 2896), 'numpy.isnan', 'numpy.isnan', (['ftrans'], {}), '(ftrans)\n', (2888, 2896), False, 'import numpy\n'), ((8516, 8541), 'astropy.units.spectral_density', 'u.spectral_density', (['vwave'], {}), '(vwave)\n', (8534, 8541), True, 'from astropy import units as u\n'), ((11115, 11137), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (11127, 11137), False, 'import numpy\n'), ((13562, 13581), 'numpy.nanmin', 'numpy.nanmin', (['fwave'], {}), '(fwave)\n', (13574, 13581), False, 'import numpy\n'), ((13608, 13627), 'numpy.nanmin', 'numpy.nanmin', (['fwave'], {}), '(fwave)\n', (13620, 13627), False, 'import numpy\n'), ((13654, 13673), 'numpy.nanmax', 'numpy.nanmax', (['fwave'], {}), '(fwave)\n', (13666, 13673), False, 'import numpy\n'), ((13700, 13719), 'numpy.nanmax', 'numpy.nanmax', (['fwave'], {}), '(fwave)\n', (13712, 13719), False, 'import numpy\n'), ((15260, 15278), 'scipy.integrate.trapz', 'trapz', (['(ft * fw)', 'fw'], {}), '(ft * fw, fw)\n', (15265, 15278), False, 'from scipy.integrate import trapz\n'), ((15276, 15289), 'scipy.integrate.trapz', 'trapz', (['ft', 'fw'], {}), '(ft, fw)\n', (15281, 15289), False, 'from scipy.integrate import trapz\n'), ((15488, 15503), 'numpy.max', 'numpy.max', (['fw05'], {}), '(fw05)\n', (15497, 15503), False, 'import numpy\n'), ((15504, 15519), 'numpy.min', 'numpy.min', (['fw05'], {}), '(fw05)\n', (15513, 15519), False, 'import numpy\n'), ((21647, 21705), 'numpy.arange', 'numpy.arange', (['(notch[0] - 5.0 * dn)', '(notch[1] + 5.0 * dn)', 'dn'], {}), '(notch[0] - 5.0 * dn, notch[1] + 5.0 * dn, dn)\n', (21659, 21705), False, 'import numpy\n'), ((22072, 22091), 'numpy.isnan', 'numpy.isnan', (['ftrans'], {}), '(ftrans)\n', (22083, 22091), False, 'import numpy\n'), ((22114, 22133), 'numpy.isnan', 'numpy.isnan', (['ftrans'], {}), '(ftrans)\n', (22125, 22133), False, 'import numpy\n'), ((22337, 22378), 'os.path.normpath', 'os.path.normpath', (['(filterFolder + vegaFile)'], {}), '(filterFolder + vegaFile)\n', (22353, 22378), False, 'import os\n'), ((23068, 23090), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (23080, 23090), False, 'import numpy\n'), ((26486, 26509), 'matplotlib.pyplot.plot', 'plt.plot', (['fwave', 'ftrans'], {}), '(fwave, ftrans)\n', (26494, 26509), True, 'import matplotlib.pyplot as plt\n'), ((27494, 27517), 'matplotlib.pyplot.plot', 'plt.plot', (['fwave', 'ftrans'], {}), '(fwave, ftrans)\n', (27502, 27517), True, 'import matplotlib.pyplot as plt\n'), ((28166, 28186), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (28178, 28186), False, 'import numpy\n'), ((35268, 35338), 'numpy.log10', 'numpy.log10', (['(1.0 + 10.0 ** (-0.4 * (other.magnitude - self.magnitude)))'], {}), '(1.0 + 10.0 ** (-0.4 * (other.magnitude - self.magnitude)))\n', (35279, 35338), False, 'import numpy\n'), ((2826, 2845), 'numpy.isnan', 'numpy.isnan', (['ftrans'], {}), '(ftrans)\n', (2837, 2845), False, 'import numpy\n'), ((9654, 9666), 'astropy.units.spectral', 'u.spectral', ([], {}), '()\n', (9664, 9666), True, 'from astropy import units as u\n'), ((9712, 9739), 'astropy.units.spectral_density', 'u.spectral_density', (['sp.wave'], {}), '(sp.wave)\n', (9730, 9739), True, 'from astropy import units as u\n'), ((9790, 9817), 'astropy.units.spectral_density', 'u.spectral_density', (['sp.wave'], {}), '(sp.wave)\n', (9808, 9817), True, 'from astropy import units as u\n'), ((9864, 9876), 'astropy.units.spectral', 'u.spectral', ([], {}), '()\n', (9874, 9876), True, 'from astropy import units as u\n'), ((11783, 11805), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (11795, 11805), False, 'import numpy\n'), ((13417, 13435), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (13428, 13435), False, 'import numpy\n'), ((13465, 13483), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (13476, 13483), False, 'import numpy\n'), ((15340, 15358), 'scipy.integrate.trapz', 'trapz', (['(fw * ft)', 'fw'], {}), '(fw * ft, fw)\n', (15345, 15358), False, 'from scipy.integrate import trapz\n'), ((15356, 15374), 'scipy.integrate.trapz', 'trapz', (['(ft / fw)', 'fw'], {}), '(ft / fw, fw)\n', (15361, 15374), False, 'from scipy.integrate import trapz\n'), ((15420, 15433), 'numpy.max', 'numpy.max', (['fw'], {}), '(fw)\n', (15429, 15433), False, 'import numpy\n'), ((15434, 15447), 'numpy.min', 'numpy.min', (['fw'], {}), '(fw)\n', (15443, 15447), False, 'import numpy\n'), ((23491, 23537), 'scipy.integrate.trapz', 'trapz', (['(ftrans.value * fwave.value)', 'fwave.value'], {}), '(ftrans.value * fwave.value, fwave.value)\n', (23496, 23537), False, 'from scipy.integrate import trapz\n'), ((23560, 23592), 'scipy.integrate.trapz', 'trapz', (['ftrans.value', 'fwave.value'], {}), '(ftrans.value, fwave.value)\n', (23565, 23592), False, 'from scipy.integrate import trapz\n'), ((23755, 23777), 'numpy.arange', 'numpy.arange', (['nsamples'], {}), '(nsamples)\n', (23767, 23777), False, 'import numpy\n'), ((27461, 27477), 'numpy.ones', 'numpy.ones', (['(1000)'], {}), '(1000)\n', (27471, 27477), False, 'import numpy\n'), ((35646, 35691), 'numpy.log10', 'numpy.log10', (['(1.0 + 10.0 ** (-0.4 * (m2 - m1)))'], {}), '(1.0 + 10.0 ** (-0.4 * (m2 - m1)))\n', (35657, 35691), False, 'import numpy\n'), ((6645, 6700), 'numpy.logical_and', 'numpy.logical_and', (['(fwave >= notch[0])', '(fwave <= notch[1])'], {}), '(fwave >= notch[0], fwave <= notch[1])\n', (6662, 6700), False, 'import numpy\n'), ((8365, 8383), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (8376, 8383), False, 'import numpy\n'), ((8417, 8435), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (8428, 8435), False, 'import numpy\n'), ((10405, 10423), 'numpy.log10', 'numpy.log10', (['(a / b)'], {}), '(a / b)\n', (10416, 10423), False, 'import numpy\n'), ((10627, 10667), 'scipy.integrate.trapz', 'trapz', (['(ftrans * fwave.value)', 'fwave.value'], {}), '(ftrans * fwave.value, fwave.value)\n', (10632, 10667), False, 'from scipy.integrate import trapz\n'), ((10931, 10957), 'scipy.integrate.trapz', 'trapz', (['ftrans', 'fwave.value'], {}), '(ftrans, fwave.value)\n', (10936, 10957), False, 'from scipy.integrate import trapz\n'), ((11583, 11602), 'astropy.constants.h.to', 'const.h.to', (['"""erg s"""'], {}), "('erg s')\n", (11593, 11602), True, 'from astropy import constants as const\n'), ((11603, 11625), 'astropy.constants.c.to', 'const.c.to', (['"""micron/s"""'], {}), "('micron/s')\n", (11613, 11625), True, 'from astropy import constants as const\n'), ((13066, 13121), 'numpy.logical_and', 'numpy.logical_and', (['(fwave >= notch[0])', '(fwave <= notch[1])'], {}), '(fwave >= notch[0], fwave <= notch[1])\n', (13083, 13121), False, 'import numpy\n'), ((18105, 18125), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (18117, 18125), False, 'import numpy\n'), ((18170, 18190), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (18182, 18190), False, 'import numpy\n'), ((18235, 18255), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (18247, 18255), False, 'import numpy\n'), ((21771, 21826), 'numpy.logical_and', 'numpy.logical_and', (['(fwave >= notch[0])', '(fwave <= notch[1])'], {}), '(fwave >= notch[0], fwave <= notch[1])\n', (21788, 21826), False, 'import numpy\n'), ((23197, 23217), 'numpy.nanstd', 'numpy.nanstd', (['result'], {}), '(result)\n', (23209, 23217), False, 'import numpy\n'), ((23437, 23462), 'astropy.units.spectral_density', 'u.spectral_density', (['fwave'], {}), '(fwave)\n', (23455, 23462), True, 'from astropy import units as u\n'), ((26449, 26469), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (26461, 26469), False, 'import numpy\n'), ((26550, 26575), 'numpy.nanmin', 'numpy.nanmin', (['fwave.value'], {}), '(fwave.value)\n', (26562, 26575), False, 'import numpy\n'), ((26576, 26601), 'numpy.nanmax', 'numpy.nanmax', (['fwave.value'], {}), '(fwave.value)\n', (26588, 26601), False, 'import numpy\n'), ((27420, 27452), 'numpy.linspace', 'numpy.linspace', (['f[0]', 'f[1]', '(1000)'], {}), '(f[0], f[1], 1000)\n', (27434, 27452), False, 'import numpy\n'), ((27558, 27583), 'numpy.nanmin', 'numpy.nanmin', (['fwave.value'], {}), '(fwave.value)\n', (27570, 27583), False, 'import numpy\n'), ((27584, 27609), 'numpy.nanmax', 'numpy.nanmax', (['fwave.value'], {}), '(fwave.value)\n', (27596, 27609), False, 'import numpy\n'), ((10694, 10748), 'scipy.integrate.trapz', 'trapz', (['(ftrans * fwave.value * fwave.value)', 'fwave.value'], {}), '(ftrans * fwave.value * fwave.value, fwave.value)\n', (10699, 10748), False, 'from scipy.integrate import trapz\n'), ((10986, 11026), 'scipy.integrate.trapz', 'trapz', (['(ftrans * fwave.value)', 'fwave.value'], {}), '(ftrans * fwave.value, fwave.value)\n', (10991, 11026), False, 'from scipy.integrate import trapz\n'), ((22512, 22530), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (22523, 22530), False, 'import numpy\n'), ((22568, 22586), 'numpy.isnan', 'numpy.isnan', (['vflux'], {}), '(vflux)\n', (22579, 22586), False, 'import numpy\n'), ((23884, 23904), 'numpy.nanstd', 'numpy.nanstd', (['result'], {}), '(result)\n', (23896, 23904), False, 'import numpy\n'), ((26647, 26672), 'numpy.nanmin', 'numpy.nanmin', (['fwave.value'], {}), '(fwave.value)\n', (26659, 26672), False, 'import numpy\n'), ((26696, 26721), 'numpy.nanmax', 'numpy.nanmax', (['fwave.value'], {}), '(fwave.value)\n', (26708, 26721), False, 'import numpy\n'), ((26776, 26796), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (26788, 26796), False, 'import numpy\n'), ((27655, 27680), 'numpy.nanmin', 'numpy.nanmin', (['fwave.value'], {}), '(fwave.value)\n', (27667, 27680), False, 'import numpy\n'), ((27704, 27729), 'numpy.nanmax', 'numpy.nanmax', (['fwave.value'], {}), '(fwave.value)\n', (27716, 27729), False, 'import numpy\n'), ((27784, 27804), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (27796, 27804), False, 'import numpy\n'), ((15052, 15072), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (15064, 15072), False, 'import numpy\n'), ((15125, 15145), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (15137, 15145), False, 'import numpy\n'), ((15198, 15218), 'numpy.nanmax', 'numpy.nanmax', (['ftrans'], {}), '(ftrans)\n', (15210, 15218), False, 'import numpy\n'), ((10305, 10330), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1)'], {}), '(0, 1)\n', (10324, 10330), False, 'import numpy\n'), ((23128, 23155), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (23147, 23155), False, 'import numpy\n'), ((9272, 9299), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (9291, 9299), False, 'import numpy\n'), ((9476, 9503), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (9495, 9503), False, 'import numpy\n'), ((23815, 23842), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (23834, 23842), False, 'import numpy\n'), ((11231, 11258), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (11250, 11258), False, 'import numpy\n'), ((11410, 11437), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (11429, 11437), False, 'import numpy\n'), ((11889, 11916), 'numpy.random.normal', 'numpy.random.normal', (['(0)', '(1.0)'], {}), '(0, 1.0)\n', (11908, 11916), False, 'import numpy\n')]