How to use pysynphot - 10 common examples

import numpy as np
    from matplotlib.ticker import FixedLocator
    import matplotlib.pyplot as plt

    import astropy.table

    try:
        from .. import fitting, prep, utils
        from . import auto_script
    except:
        from grizli import prep, utils, fitting
        from grizli.pipeline import auto_script

    if filter_bandpasses is None:
        import pysynphot as S
        filter_bandpasses = [S.ObsBandpass(bpstr) for bpstr in ['acs,wfc1,f814w', 'wfc3,ir,f105w', 'wfc3,ir,f110w', 'wfc3,ir,f125w', 'wfc3,ir,f140w', 'wfc3,ir,f160w']]

    if os.path.exists('{0}.info.fits'.format(field_root)):
        orig = utils.read_catalog('{0}.info.fits'.format(field_root))
        all_files = glob.glob('{0}*full.fits'.format(field_root))
        all_files.sort()

        print('{0}.info.fits: {1} objects.  Found {2} full.fits files, checking modify dates.'.format(field_root, len(orig), len(all_files)))
        info_mtime = os.stat('{0}.info.fits'.format(field_root)).st_mtime
        keep = np.ones(len(orig), dtype=bool)

        files = []
        for file in all_files:
            id = int(file.split('_')[1].split('.full')[0])
            if id not in orig['id']:
                files.append(file)
            else:

MLv = mass_norm / Lv_norm
        
    SFR_norm = (coeffs_norm*tab_temp['sfr']).sum()*u.solMass/u.yr
    SFRv = SFR_norm / Lv_norm
    
    mass_var = ((rms_norm*tab_temp['mass'])**2).sum()
    Lv_var = ((rms_norm*tab_temp['Lv'])**2).sum()
    SFR_var = ((rms_norm*tab_temp['sfr'])**2).sum()
    
    MLv_var = MLv**2 * (mass_var/mass_norm.value**2 + Lv_var/Lv_norm.value**2)
    MLv_rms = np.sqrt(MLv_var)

    SFRv_var = SFRv**2 * (SFR_var/SFR_norm.value**2 + Lv_var/Lv_norm.value**2)
    SFRv_rms = np.sqrt(SFRv_var)
    
    vband = S.ObsBandpass('v')
    vbandz = S.ArrayBandpass(vband.wave*(1+z), vband.throughput)
    
    best_templ = utils.SpectrumTemplate(templ['wave'], templ['full'])
    fnu = best_templ.integrate_filter(vbandz)*(u.erg/u.s/u.cm**2/u.Hz)
    
    dL = cosmology.luminosity_distance(z).to(u.cm)
    Lnu = fnu*4*np.pi*dL**2
    pivotV = vbandz.pivot()*u.Angstrom
    nuV = (const.c/pivotV).to(u.Hz) 
    Lv = (nuV*Lnu).to(u.L_sun)
            
    mass = MLv*Lv
    SFR = SFRv*Lv
    
    sps = {'Lv':Lv, 'MLv':MLv, 'MLv_rms':MLv_rms, 'SFRv':SFRv, 'SFRv_rms':SFRv_rms, 'templ':best_templ}

import astropy.coordinates as coord
        
        coo = coord.SkyCoord(ra=self.ra, dec=self.dec, unit=(u.deg, u.deg),
                             frame='icrs')
                                                          
        v = Vizier(catalog=SDSS_CATALOG, columns=['+_r','*'])
        try:
            tab = v.query_region(coo, radius="{0}s".format(radius),
                              catalog=SDSS_CATALOG)[0]
            
            ix = np.argmin(tab['rmag'])
            tab = tab[ix]
        except:
            return None
                               
        filters = [FilterDefinition(bp=S.ObsBandpass('sdss,{0}'.format(b))) for b in bands]
        pivot = {}
        for ib, b in enumerate(bands):
            pivot[b] = filters[ib].pivot()
            
        #to_flam = 10**(-0.4*(22.5+48.6))*3.e18 # / pivot(Ang)**2
        #flam = np.array([xid['modelFlux_{0}'.format(b)][0]*to_flam/pivot[b]**2 for b in bands])
        #eflam = np.array([np.sqrt(1/xid['modelFluxIvar_{0}'.format(b)][0])*to_flam/pivot[b]**2 for b in bands])
        
        to_flam = 10**(-0.4*(48.6))*3.e18 # / pivot(Ang)**2
        flam = np.array([10**(-0.4*(tab[b+'mag']))*to_flam/pivot[b]**2 for ib, b in enumerate(bands)])
        eflam = np.array([tab['e_{0}mag'.format(b)]*np.log(10)/2.5*flam[ib] for ib, b in enumerate(bands)])
        
        phot = {'flam':flam, 'eflam':eflam, 'filters':filters, 'tempfilt':None}
        
        if templ is None:
            return phot

"""
    import os
    import numpy as np
    from matplotlib.ticker import FixedLocator
    import matplotlib.pyplot as plt
    
    try:
        from .. import fitting
        from . import auto_script
    except:
        from grizli import prep, utils, fitting
        from grizli.pipeline import auto_script
    
    if filter_bandpasses is None:
        import pysynphot as S
        filter_bandpasses = [S.ObsBandpass(bpstr) for bpstr in ['acs,wfc1,f814w','wfc3,ir,f105w', 'wfc3,ir,f110w', 'wfc3,ir,f125w', 'wfc3,ir,f140w', 'wfc3,ir,f160w']]
        
    ### SUmmary catalog
    fit = fitting.make_summary_catalog(target=field_root, sextractor=None,
                                       filter_bandpasses=filter_bandpasses)
                                       
    fit.meta['root'] = field_root
    
    ## Add photometric catalog
    sex = utils.GTable.gread('{0}-ir.cat.fits'.format(field_root))
    # try:
    # except:
    #     sex = utils.GTable.gread('../Prep/{0}-ir.cat.fits'.format(field_root), sextractor=True)
        
    idx = np.arange(len(sex))
    sex_idx = np.array([idx[sex['NUMBER'] == id][0] for id in fit['id']])

'F775W': 7693.297933335407,
            'F814W': 8058.784799323767,
            'VISTAH':1.6433e+04,
            'GRISM': 1.6e4}

no_newline = '\x1b[1A\x1b[1M' # character to skip clearing line on STDOUT printing

### Demo for computing photflam and photplam with pysynphot
if False:
    import pysynphot as S
    n = 1.e-20
    spec = S.FlatSpectrum(n, fluxunits='flam')
    photflam = {}
    photplam = {}
    for filter in ['F098M', 'F105W', 'F110W', 'F125W', 'F140W', 'F160W']:
        bp = S.ObsBandpass('wfc3,ir,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    for filter in ['F435W', 'F606W', 'F775W', 'F814W']:
        bp = S.ObsBandpass('acs,wfc1,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    
class GrismFLT(object):
    """
    Scripts for simple modeling of individual grism FLT images
    
    tbd: 

>>> import pysynphot as S
           >>> bandpasses = [S.ObsBandpass('wfc3,ir,f140w')]
    
    
    Returns
    -------
    ab_mags : dict
        Dictionary with keys from `bandpasses` and the integrated magnitudes
        
    """
    import pysynphot as S
    flat = S.FlatSpectrum(0, fluxunits='ABMag')
    ab_mags = OrderedDict()
    
    for bp in bandpasses:
        flat_obs = S.Observation(flat, bp)
        spec_obs = S.Observation(spectrum, bp)
        ab_mags[bp.name] = -2.5*np.log10(spec_obs.countrate()/flat_obs.countrate())
    
    return ab_mags

def calc_fluxes():
    return np.array([S.Observation(sp,i).effstim("flam") for i in pfilts])

if False:
    import pysynphot as S
    n = 1.e-20
    spec = S.FlatSpectrum(n, fluxunits='flam')
    photflam = {}
    photplam = {}
    for filter in ['F098M', 'F105W', 'F110W', 'F125W', 'F140W', 'F160W']:
        bp = S.ObsBandpass('wfc3,ir,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    for filter in ['F435W', 'F606W', 'F775W', 'F814W']:
        bp = S.ObsBandpass('acs,wfc1,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    
class GrismFLT(object):
    """
    Scripts for simple modeling of individual grism FLT images
    
    tbd: 
        o helper functions for extracting 2D spectra
        o lots of book-keeping for handling SExtractor objects & catalogs
        ...
        
    """
    def __init__(self, flt_file='ico205lwq_flt.fits', sci_ext=('SCI',1),
                 direct_image=None, refimage=None, segimage=None, refext=0,
                 verbose=True, pad=100, shrink_segimage=True,

>>> bandpasses = [S.ObsBandpass('wfc3,ir,f140w')]
    
    
    Returns
    -------
    ab_mags : dict
        Dictionary with keys from `bandpasses` and the integrated magnitudes
        
    """
    import pysynphot as S
    flat = S.FlatSpectrum(0, fluxunits='ABMag')
    ab_mags = OrderedDict()
    
    for bp in bandpasses:
        flat_obs = S.Observation(flat, bp)
        spec_obs = S.Observation(spectrum, bp)
        ab_mags[bp.name] = -2.5*np.log10(spec_obs.countrate()/flat_obs.countrate())
    
    return ab_mags

'VISTAH':1.6433e+04,
            'GRISM': 1.6e4}

no_newline = '\x1b[1A\x1b[1M' # character to skip clearing line on STDOUT printing

### Demo for computing photflam and photplam with pysynphot
if False:
    import pysynphot as S
    n = 1.e-20
    spec = S.FlatSpectrum(n, fluxunits='flam')
    photflam = {}
    photplam = {}
    for filter in ['F098M', 'F105W', 'F110W', 'F125W', 'F140W', 'F160W']:
        bp = S.ObsBandpass('wfc3,ir,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    for filter in ['F435W', 'F606W', 'F775W', 'F814W']:
        bp = S.ObsBandpass('acs,wfc1,%s' %(filter.lower()))
        photplam[filter] = bp.pivot()
        obs = S.Observation(spec, bp)
        photflam[filter] = n/obs.countrate()
        
    
class GrismFLT(object):
    """
    Scripts for simple modeling of individual grism FLT images
    
    tbd: 
        o helper functions for extracting 2D spectra
        o lots of book-keeping for handling SExtractor objects & catalogs