def load_models(self):
     """Load the BaSeL basis models, initialize the grid to hold resampled models,
     and load the filters"""
     self.basel = starmodel.BaSeL3()  # BaSeL 3.1 Basis
     self.basel.read_all_Z()
     self.stargrid = starmodel.SpecLibrary() #object to hold the model grid
     self.model_filterlist = observate.load_filters(self.rp['model_fnamelist']) #filter objects
     self.fit_filterlist = observate.load_filters(self.rp['fit_fnamelist'])

def selftest():
    from sedpy.observate import load_filters
    sps = sps_basis.StellarPopBasis(debug=True)
    params = {}
    params['tage'] = np.array([1, 2, 3, 4.])
    params['zmet'] = np.array([-0.5, 0.0])
    ntot = len(params['tage']) * len(params['zmet'])
    params['mass'] = np.random.uniform(0, 1, ntot)
    params['sigma_smooth'] = 100.
    outwave = sps.ssp.wavelengths
    flist = ['sdss_u0', 'sdss_r0']
    filters = load_filters(flist)

    # get a spectrum
    s, p, e = sps.get_spectrum(params, outwave, filters)
    # change parameters that affect neither the basis nor the ssp, and
    # get spectrum again
    params['mass'] = np.random.uniform(0, 1, ntot)
    s, p, e = sps.get_spectrum(params, outwave, filters)
    # lets get the basis components while we're at it
    bs, bp, be = sps.get_components(params, outwave, filters)
    # change something that affects the basis
    params['tage'] += 1.0
    bs, bp, be = sps.get_components(params, outwave, filters)
    # try a single age pop at arbitrary metallicity
    params['tage'] = 1.0
    params['zmet'] = -0.2
    bs, bp, be = sps.get_components(params, outwave, filters)

def load_obs(objid=0, phottable='demo_photometry.dat', **kwargs):
    """Load photometry from an ascii file.  Assumes the following columns:
    `objid`, `filterset`, [`mag0`,....,`magN`] where N >= 11.  The User should
    modify this function (including adding keyword arguments) to read in their
    particular data format and put it in the required dictionary.

    :param objid:
        The object id for the row of the photomotery file to use.  Integer.
        Requires that there be an `objid` column in the ascii file.

    :param phottable:
        Name (and path) of the ascii file containing the photometry.

    :returns obs:
        Dictionary of observational data.
    """
    # Writes your code here to read data.  Can use FITS, h5py, astropy.table,
    # sqlite, whatever.
    # e.g.:
    # import astropy.io.fits as pyfits
    # catalog = pyfits.getdata(phottable)

    # Here we will read in an ascii catalog of magnitudes as a numpy structured
    # array
    with open(phottable, 'r') as f:
        # drop the comment hash
        header = f.readline().split()[1:]
    catalog = np.genfromtxt(phottable, comments='#',
                            dtype=np.dtype([(n, np.float) for n in header]))

    # Find the right row
    ind = catalog['objid'] == float(objid)
    # Here we are dynamically choosing which filters to use based on the object
    # and a flag in the catalog.  Feel free to make this logic more (or less)
    # complicated.
    filternames = filtersets[ int(catalog[ind]['filterset']) ]
    # And here we loop over the magnitude columns
    mags = [catalog[ind]['mag{}'.format(i)] for i in range(len(filternames))]
    mags = np.array(mags)

    # Build output dictionary. 
    obs = {}
    # This is a list of sedpy filter objects.    See the
    # sedpy.observate.load_filters command for more details on its syntax.
    obs['filters'] = load_filters(filternames)
    # This is a list of maggies, converted from mags.  It should have the same
    # order as `filters` above.
    obs['maggies'] = np.squeeze(10**(-mags/2.5))
    # Hack.  you should use real flux uncertainties
    obs['maggies_unc'] = obs['maggies'] * 0.07
    # Here we mask out any NaNs or infs
    obs['phot_mask'] = np.isfinite(np.squeeze(mags))
    # We have no spectrum.
    obs['wavelength'] = None

    # Add unessential bonus info.  This will be sored in output
    #obs['dmod'] = catalog[ind]['dmod']
    obs['objid'] = objid

    return obs

 def load_data(self):
     """Read the catalogs, apply distance modulus,
     and determine 'good' pixels"""
     self.data_filterlist = observate.load_filters(self.rp['data_fnamelist'])
     self.data_mag, self.data_magerr, self.data_header = catio.load_image_cube(self.rp)
     self.distance_modulus = 5.0*np.log10(self.rp['dist'])+25
     self.nobj = self.data_mag.shape[0]

def load_obs(**extras):
    """the key word will be "uvj_key"
    generate UVJ fluxes based off of this
    """

    # first load the filters
    fnames = ['bessell_U','bessell_V','twomass_J']

    # generate uvj color
    uv, vj = return_uvj(int(extras['uvj_key']))

    # translate to fluxes
    # do this by fixing V-flux in maggies
    # this corresponds roughly to M~1e10
    vflux = 5e7
    maggies = np.array([10**(-uv/2.5), 1, 10**(vj/2.5)]) * vflux
    maggies_unc = maggies / 40.

    ### build output dictionary
    obs = {}
    obs['filters'] = observate.load_filters(fnames)
    obs['wave_effective'] = np.array([filt.wave_effective for filt in obs['filters']])
    obs['phot_mask'] = np.ones_like(maggies,dtype=bool)
    obs['maggies'] = maggies
    obs['maggies_unc'] =  maggies_unc
    obs['wavelength'] = None
    obs['spectrum'] = None
    obs['logify_spectrum'] = False
    obs['uv'] = uv
    obs['vj'] = vj

    return obs

def load_obs(photname='', objname='', **extras):
    """
    Custom-built because the photometric files are actually generated by the model
    """
    obs ={}

    # if the photometric files exist
    with open(photname, 'r') as f:
        hdr = f.readline().split()
    dat = np.loadtxt(photname, comments = '#',
                     dtype = np.dtype([(n, np.float) for n in hdr[1:]]))
    obj_ind = np.where(dat['id'] == int(objname))[0][0]
    
    # extract fluxes+uncertainties for all objects and all filters
    flux_fields = [f for f in dat.dtype.names if f[0:2] == 'f_']
    unc_fields = [f for f in dat.dtype.names if f[0:2] == 'e_']
    filters = [f[2:] for f in flux_fields]

    # extract fluxes for particular object, converting from record array to numpy array
    flux = dat[flux_fields].view(float).reshape(len(dat),-1)[obj_ind]
    unc  = dat[unc_fields].view(float).reshape(len(dat),-1)[obj_ind]

    # build output dictionary
    obs['filters'] = observate.load_filters(filters)
    obs['wave_effective'] = np.array([filt.wave_effective for filt in obs['filters']])
    obs['phot_mask'] = np.ones_like(flux,dtype=bool)
    obs['maggies'] = flux
    obs['maggies_unc'] =  unc
    obs['wavelength'] = None
    obs['spectrum'] = None
    obs['logify_spectrum'] = False

    return obs

def load_obs(objname=None, datloc=None, err_floor=0.05, **extras):

    ''' 
    objname: number of object in the 3D-HST COSMOS photometric catalog
    err_floor: the fractional error floor (0.05 = 5% floor)
    zp_err: inflate the errors by the zeropoint offsets from Skelton+14
    '''

    ### open file, load data
    with open(datloc, 'r') as f:
        hdr = f.readline().split()
    dtype = np.dtype([(hdr[1],'S20')] + [(n, np.float) for n in hdr[2:]])
    dat = np.loadtxt(datloc, comments = '#', dtype = dtype)

    ### extract filters, fluxes, errors for object
    obj_idx = (dat['ID'] == objname)
    filters = np.array([f[1:] for f in dat.dtype.names if f[0] == 'e'])
    flux = np.squeeze([dat[obj_idx][f] for f in filters])
    unc = np.squeeze([dat[obj_idx]['e'+f] for f in filters])

    ### define photometric mask, convert from Jy to maggies
    phot_mask = (flux != unc) & (flux != -99.0) & (unc > 0)
    maggies = flux / 3631
    maggies_unc = unc / 3631

    ### implement error floor
    maggies_unc = np.clip(maggies_unc, maggies*err_floor, np.inf)

    ### mask anything touching or bluewards of Ly-a
    zred = dat['zz'][obj_idx]
    fnames = [f+'_candels' for f in filters]
    translate = {
                 'f24_candels':'spitzer_mips_24',
                 'f70_candels':'spitzer_mips_70',
                 'f100_candels':'herschel_pacs_100',
                 'f160_candels':'herschel_pacs_160',
                 'f250_candels':'herschel_spire_250'
                }
    fnames = [translate[f] if f in translate.keys() else f for f in fnames]
    ofilters = observate.load_filters(fnames)

    wavemax = np.array([f.wavelength[f.transmission > (f.transmission.max()*0.1)].max() for f in ofilters]) / (1+zred)
    wavemin = np.array([f.wavelength[f.transmission > (f.transmission.max()*0.1)].min() for f in ofilters]) / (1+zred)
    filtered = [1230]
    for f in filtered: phot_mask[(wavemax > f) & (wavemin < f)] = False
    phot_mask[wavemin < 1200] = False

    ### build output dictionary
    obs = {}
    obs['filters'] = ofilters
    obs['wave_effective'] = np.array([filt.wave_effective for filt in obs['filters']])
    obs['phot_mask'] = phot_mask
    obs['maggies'] = maggies
    obs['maggies_unc'] =  maggies_unc
    obs['wavelength'] = None
    obs['spectrum'] = None
    obs['logify_spectrum'] = False

    return obs

def read_hdf5(filename, **extras):
    """Read an HDF5 file (with a specific format) into a dictionary of results.

    This HDF5 file is assumed to have the groups ``sampling`` and ``obs`` which
    respectively contain the sampling chain and the observational data used in
    the inference.

    All attributes of these groups as well as top-level attributes are loaded
    into the top-level of the dictionary using ``json.loads``, and therefore
    must have been written with ``json.dumps``.  This should probably use
    JSONDecoders, but who has time to learn that.

    :param filename:
        Name of the HDF5 file.
    """
    groups = {'sampling': {}, 'obs': {}}
    res = {}
    with h5py.File(filename, "r") as hf:
        # loop over the groups
        for group, d in groups.items():
            # read the arrays in that group into the dictionary for that group
            for k, v in hf[group].items():
                d[k] = np.array(v)
            # unserialize the attributes and put them in the dictionary
            for k, v in hf[group].attrs.items():
                try:
                    d[k] = json.loads(v)
                except:
                    try:
                        d[k] = unpick(v)
                    except:
                        d[k] = v
        # do top-level attributes.
        for k, v in hf.attrs.items():
            try:
                res[k] = json.loads(v)
            except:
                try:
                    res[k] = unpick(v)
                except:
                    res[k] = v
        res.update(groups['sampling'])
        res['obs'] = groups['obs']
        try:
            res['obs']['filters'] = load_filters([str(f) for f in res['obs']['filters']])
        except:
            pass
        try:
            res['rstate'] = unpick(res['rstate'])
        except:
            pass
        try:
            mp = [names_to_functions(p.copy()) for p in res['model_params']]
            res['model_params'] = mp
        except:
            pass

    return res

def load_obs_mmt(filename=None, objname=None, #dist = 1e-5, vel = 0.0,
                  wlo=3750., whi=7200., verbose=False,
                  phottable='data/f2_apcanfinal_6phot_v2.fits',
                  **kwargs):
    """
    Read one of Caldwell's MMT spectra and find the matching PHAT
    photometry, return a dictionary containing the observations.
    """
    from sedpy import observate

    obs ={}

    ####### SPECTRUM #######
    if verbose:
        print('Loading data from {0}'.format(filename))

    scale = 1e0 #testing
    #fluxconv = np.pi * 4. * (dist * 1e6 * pc)**2 / lsun #erg/s/AA/cm^2 to L_sun/AA
    fluxconv =  1.0#5.0e-20 * scale #approximate counts to erg/s/AA/cm^2
    #redshift = 0.0 #vel / 2.998e8
    dat = np.squeeze(pyfits.getdata(filename))
    hdr = pyfits.getheader(filename)
    
    crpix = (hdr['CRPIX1'] -1) #convert from FITS to numpy indexing
    try:
        cd = hdr['CDELT1']
    except (KeyError):
        cd = hdr['CD1_1']

    obs['wavelength'] = (np.arange(dat.shape[1]) - crpix) * cd + hdr['CRVAL1']
    obs['spectrum'] = dat[0,:] * fluxconv
    obs['unc'] = np.sqrt(dat[1,:]) * fluxconv
    
    #Masking.  should move to a function that reads a mask definition file
    #one should really never mask in the rest frame - that should be modeled!
    obs['mask'] =  ((obs['wavelength'] >= wlo ) & (obs['wavelength'] <= whi))
    obs['mask'] = obs['mask'] & ((obs['wavelength'] <= 5570) |
                                 (obs['wavelength'] >= 5590)) #mask OI sky line
    obs['mask'] = obs['mask'] & ((obs['wavelength'] <= 6170) |
                                 (obs['wavelength'] >= 6180)) #mask...something.

    #obs['wavelength'] /= (1.0 + redshift)

    ######## PHOTOMETRY ########
    if verbose:
        print('Loading mags from {0} for {1}'.format(phottable, objname))
    mags, mags_unc, flag = query_phatcat(objname, phottable = phottable, **kwargs)
    
    obs['filters'] = observate.load_filters(['wfc3_uvis_'+b.lower() for b in
                                             ["F275W", "F336W", "F475W", "F814W"]] +
                                             ['wfc3_ir_'+b.lower() for b in
                                              ["F110W", "F160W"]])
    obs['maggies'] = 10**(-0.4 * (mags -
                                  np.array([f.ab_to_vega for f in obs['filters']]) -
                                  2.5*np.log10(scale) ))
    obs['maggies_unc'] = mags_unc * obs['maggies'] / 1.086

    return obs

def zsfh_to_obs(sfhlist, zlist, lfbandnames=None, select_function=None,
                bandnames=None, sps=None, isocs=None, **kwargs):
    """
    Go from a list of SFHs (one for each metallicity) to a broadband SED and set of
    luminosity functions for a stellar population.
    """
    sed_values, lf_values = {}, []
    #basti = np.any(sps.zlegend[0:2] == 0.0006) #Hack to check for Basti Isochrones
    nsfh = len(sfhlist)
    assert len(zlist) == nsfh

    ###########
    # SED
    ############
    if bandnames is not None:
        filterlist = observate.load_filters(bandnames)
        spec, wave, mass = rsed.one_region_sed(copy.deepcopy(sfhlist), total_zmet, sps)
        mags = observate.getSED(wave, spec*rsed.to_cgs, filterlist=filterlist)
        maggies = 10**(-0.4 * np.atleast_1d(mags))
        sed_values['sed_ab_maggies'] = maggies
        sed_values['sed_filters'] = bandnames

    #############
    # LFs
    ############
    #create the SSP CLFs, using nearest neighbor interpolation for the metallicity
    all_lf_base = []
    bins = rsed.lfbins
    for i,zmet in enumerate(zlist):
        if isocs is not None:
            isoc = isocs[i]
        else:
            sps.params['zmet'] = np.abs(sps.zlegend - zmet).argmin() + 1
            isoc = sps.isochrones()
            print("Using Zmet={0} in place of requested "
            "Zmet={1}".format(sps.zlegend[sps.params['zmet']-1],zmet))

        ldat = isochrone_to_clfs(copy.deepcopy(isoc), lfbandnames,
                                 select_function=select_function, **kwargs)
        all_lf_base += [ldat]
    #use the SSP CLFs to generate a total LF (for each band)
    for i, band in enumerate(lfbandnames):
        lf_oneband = {}
        lf_base = [zbase[i] for zbase in all_lf_base]
        lfs_zt, lf, logages = rsed.one_region_lfs(copy.deepcopy(sfhlist), lf_base)
        lf_oneband['bandname'] = band
        lf_oneband['clf'] = lf
        lf_oneband['clf_mags'] = bins
        lf_oneband['logages'] = logages
        lf_oneband['clfs_zt'] = lfs_zt
        
        lf_values += [lf_oneband]
        
    #############
    # Write output
    ############
    return sed_values, lf_values

 def load_models(self):
     """Load the Draine & Li basis models, initialize the grid to
     hold resampled models, and load the filters
     """
     # Draine and Li Basis
     self.dl07 = dustmodel.DraineLi()
     # object to hold the model grid
     self.dustgrid = dustmodel.SpecLibrary()
     # filter objects
     self.filterlist = observate.load_filters(self.rp['fnamelist'])

def load_obs(snr=10.0, filterset=["sdss_g0", "sdss_r0"],
             add_noise=True, **kwargs):
    """Make a mock dataset.  Feel free to add more complicated kwargs, and put
    other things in the run_params dictionary to control how the mock is
    generated.

    :param snr:
        The S/N of the phock photometry.  This can also be a vector of same
        lngth as the number of filters.

    :param filterset:
        A list of `sedpy` filter names.  Mock photometry will be generated
        for these filters.

    :param add_noise: (optional, boolean, default: True)
        If True, add a realization of the noise to the mock spectrum
    """
    # We'll put the mock data in this dictionary, just as we would for real
    # data.  But we need to know which bands (and wavelengths if doing
    # spectroscopy) in which to generate mock data.
    mock = {}
    mock['wavelength'] = None # No spectrum
    mock['filters'] = load_filters(filterset)
    # 

    # We need the models to make a mock
    sps = load_sps(**kwargs)
    mod = load_model(**kwargs)

    # Now we get the mock params from the kwargs dict
    params = {}
    for p in mod.params.keys():
        if p in kwargs:
            params[p] = np.atleast_1d(kwargs[p])

    # And build the mock
    mod.params.update(params)
    spec, phot, _ = mod.mean_model(mod.theta, mock, sps=sps)
    # Now store some output
    mock['true_spectrum'] = spec.copy()
    mock['true_maggies'] = phot.copy()
    mock['mock_params'] = deepcopy(mod.params)
    # And add noise
    pnoise_sigma = phot / snr
    if add_noise:
        pnoise = np.random.normal(0, 1, len(phot)) * pnoise_sigma
        mock['maggies'] = phot + pnoise
    else:
        mock['maggies'] = phot.copy()
    mock['maggies_unc'] = pnoise_sigma
    mock['mock_snr'] = snr
    mock['phot_mask'] = np.ones(len(phot), dtype=bool)

    return mock

def load_obs(filter_key=1, **extras):

    # what are our mock parameters?
    params = {
              'logmass': 10,
              'logzsol': -0.75,
              'logsfr_ratios': np.array([0.15,0.15,0.15,0.15,0.15,0.15]), # hacked, gently rising
              'dust2': 0.5,
              'dust1_fraction': 0.9,
              'fagn': 0.05,
              'agn_tau': 20,
              'duste_qpah': 0.5,
              'duste_umin': 2,
              'duste_gamma': 0.2,
              'dust_index': -0.3
              }
    for p in params.keys(): params[p] = np.atleast_1d(params[p])

    # what are our filters?
    filters = find_filters(filter_key)

    # load model, sps
    mod = load_model(**run_params)
    sps = load_sps(**run_params)

    # we will also need an obs dictionary
    obs = {}
    obs['filters'] = observate.load_filters(filters)
    obs['wavelength'] = None

    # signal to noise ratio
    snr = 20
    if filter_key == 10:
        snr = np.array([np.repeat(10,10).tolist() + [4]])

    # Generate the photometry, add noise
    mod.params.update(params)
    spec, maggies, _ = mod.mean_model(mod.theta, obs, sps=sps)
    maggies_unc = np.atleast_1d((maggies / snr).squeeze())
    phot_mask = np.ones_like(maggies,dtype=bool)

    ### build output dictionary
    obs['wave_effective'] = np.array([filt.wave_effective for filt in obs['filters']])
    obs['phot_mask'] = phot_mask
    obs['maggies'] = maggies
    obs['maggies_unc'] =  maggies_unc
    obs['spectrum'] = None
    obs['logify_spectrum'] = False
    obs['true_params'] = params
    obs['true_spec'] = spec

    return obs

def median_by_band(x,y,avg=False,log=False):

    ##### get filter effective wavelengths for sorting
    delz = 0.15
    from brownseds_np_params import translate_filters
    from sedpy import observate

    filtnames = np.array(translate_filters(0,full_list=True))
    filts = observate.load_filters(filtnames[filtnames != 'nan'])
    wave_effective = np.array([filt.wave_effective for filt in filts])/1e4
    wave_effective.sort()

    ### remove overlapping filters
    ### and ones outside the plot
    delta_lam = 0.06
    for i in range(wave_effective.shape[0]):
        if i >= wave_effective.shape[0]:
            continue
        if ((np.abs(1-wave_effective/wave_effective[i]) < delta_lam).sum() > 1):
            wave_effective = np.delete(wave_effective,i)

    wave_effective = wave_effective[wave_effective < 25]


    avglam, outval, outval_up, outval_down, outval_upup, outval_downdown = [np.array([]) for i in range(6)]
    for lam in wave_effective:
        #in_bounds = (x <= lam) & (x > lam/(1+delz))
        in_bounds = (np.abs(1-lam/x) < delta_lam)
        avglam = np.append(avglam, np.mean(x[in_bounds]))
        
        mid, top, bot, toptop, botbot = np.percentile(y[in_bounds], [50, 84, 16, 95, 5])
        outval_up = np.append(outval_up, top)
        outval_down = np.append(outval_down, bot)
        outval_upup = np.append(outval_upup, toptop)
        outval_downdown = np.append(outval_downdown, botbot)
        
        if avg:
            mean, sig = y[in_bounds].mean(), np.std(y[in_bounds])
            outval = np.append(outval, mean)
            #outval_up = np.append(outval_up, mean+sig)
            #outval_down = np.append(outval_down, mean-sig)
            #outval_upup = np.append(outval_upup, mean+2*sig)
            #outval_downdown = np.append(outval_downdown, mean-2*sig)
        else:
            outval = np.append(outval, mid)

    if log:
        return avglam, -np.log10(1-outval), -np.log10(1-outval_up), -np.log10(1-outval_down), -np.log10(1-outval_upup), -np.log10(1-outval_downdown)
    else:
        return avglam, outval, outval_up, outval_down, outval_upup, outval_downdown

def load_obs(objname=None, datdir=None, runname=None, err_floor=0.05, zperr=True, no_zp_corrs=False, **extras):

    ''' 
    objname: number of object in the 3D-HST COSMOS photometric catalog
    err_floor: the fractional error floor (0.05 = 5% floor)
    zp_err: inflate the errors by the zeropoint offsets from Skelton+14
    '''

    filters = ['f1000w_miri','f1280w_miri','f1500w_miri','f1800w_miri','f2100w_miri','mips_24um_aegis']
    obs = {'maggies': None, 
           'phot_mask': np.ones_like(filters),
           'wavelength': None, 
           'filters': observate.load_filters(filters)}
    return obs

def do_all(runname='vis', outfolder=None,**opts):

    if outfolder is None:
        outfolder = os.getenv('APPS') + '/prospector_alpha/plots/'+runname+'/vis_plots/'
        if not os.path.isdir(outfolder):
            os.makedirs(outfolder)
            os.makedirs(outfolder+'data/')

    data = collate_data(runname,filename=outfolder+'data/dat.h5',**opts)
    for i in range(len(data['filters'])): data['filters'][i] = load_filters(data['filters'][i])

    if runname == 'vis':
        data['mod'] = vis_params.load_model(**vis_params.run_params)
    else:
        data['mod'] = vis_expsfh_params.load_model(**vis_expsfh_params.run_params)
    plot_addfilts(data,outfolder,runname)

def load_obs_lris(filename=None, objname=None, #dist = 1e-5, vel = 0.0,
                  wlo=3550., whi=5500., verbose=False,
                  phottable='data/f2_apcanfinal_6phot_v2.fits',
                  **kwargs):
    """
    Read one of the Keck LRIS spectra and find the matching PHAT
    photometry, return a dictionary containing the observations.
    """
    from sedpy import observate

    obs ={}
    
    ####### SPECTRUM #######
    if verbose:
        print('Loading data from {0}'.format())

    #fluxconv = np.pi * 4. * (dist * 1e6 * pc)**2 / lsun #erg/s/AA/cm^2 to L_sun/AA
    fluxconv = 1.0
    scale = 1e0 #testing
    #redshift = vel / 2.998e8
    dat = pyfits.getdata(filename)
    hdr = pyfits.getheader(filename)
    
    obs['wavelength'] = dat[0]['wave_opt']
    obs['spectrum'] = dat[0]['spec']
    obs['unc'] = 1./np.sqrt(dat[0]['ivar'])
    #masking
    obs['mask'] =  ((obs['wavelength'] >= wlo ) & (obs['wavelength'] <= whi))
    #obs['wavelength'] /= (1.0 + redshift)
    

    ######## PHOTOMETRY ######
    if verbose:
        print('Loading mags from {0} for {1}'.format(phottable, objname))
    mags, mags_unc, flag = query_phatcat(objname, phottable = phottable, **kwargs)
     
    obs['filters'] = observate.load_filters(['wfc3_uvis_'+b.lower() for b in
                                             ["F275W", "F336W", "F475W", "F814W"]] +
                                             ['wfc3_ir_'+b.lower() for b in
                                              ["F110W", "F160W"]])
    obs['maggies'] = 10**(-0.4 * (mags -
                                  np.array([f.ab_to_vega for f in obs['filters']]) -
                                  2.5*np.log10(scale) ))
    obs['maggies_unc'] = mags_unc * obs['maggies'] / 1.086

    return obs

def load_obs(photname='', objname='', **extras):
    """
    Custom-built because the photometric files are actually generated by the model
    """
    obs ={}

    filters = translate_filters()

    # build output dictionary
    obs['filters'] = observate.load_filters(filters)
    obs['wave_effective'] = np.array([filt.wave_effective for filt in obs['filters']])
    obs['phot_mask'] = np.ones_like(filters,dtype=bool)
    obs['maggies'] = np.ones_like(filters,dtype=float)
    obs['maggies_unc'] = np.ones_like(filters,dtype=float)
    obs['wavelength'] = None
    obs['spectrum'] = None
    obs['logify_spectrum'] = False
    return obs

def rectify_basis(wave, spectra, wlow=0, whigh=np.inf,
                  exclude=[], outwave=None, filters=None, **extras):
    """Mask a spectral basis using lists of include and exclude ranges
    """
    if filters is not None:
        flist = observate.load_filters(filters)
        sed = observate.getSED(wave, spectra, filterlist=flist)
        return np.array([f.wave_effective for f in flist]), 10**(-0.4 * sed)

    if outwave is not None:
        onedinterp = interpolate(wave, spectra, axis=-1)
        spectra = onedinterp(outwave)
        wave = outwave

    g = (wave >= wlow) & (wave <= whigh)
    for (lo, hi) in exclude:
        g = g & ((wave < lo) | (wave > hi))
    return wave[g], spectra[:, g]

def observe(fnames):

    #  units: lambda (A), flux: fnu normalized to unity
    dat = load_data()
    filters = load_filters(fnames)

    out = {}
    for name in dat.dtype.names:
        if name == 'blank' or name == 'lambda':
            continue

        # sourcewave: Spectrum wavelength (in AA), ndarray of shape (nwave).
        # sourceflux: Associated flux (assumed to be in erg/s/cm^2/AA), ndarray of shape (nsource,nwave).
        # filterlist: List of filter objects, of length nfilt.
        # array of AB broadband magnitudes, of shape (nsource, nfilter).
        out[name] = getSED(dat['lambda'], (lightspeed/dat['lambda']**2)*dat[name], filters)

    return out