Get Normalized Throughputs from Pandeia¶

pandeia_code_path = '/grp/jwst/jwstetc/pandeia/engine'
sys.path.append(pandeia_code_path)  # required to access installed pandeia codebase

import pandeia, pandeia.engine, pandeia.engine.config

# Set default path to find the reference files
# THis must be done BEFORE importnig pandeia.engine.jwst
pandeia.engine.config.default_refdata('/grp/jwst/jwstetc/pandeia_data/')

import pandeia.engine.jwst
from pandeia.engine.instrument_factory import InstrumentFactory

def get_pce(inst, normalize=False):
    wr = inst.get_wave_range()
    wave = np.linspace(wr['wmin']*0.95, wr['wmax']*1.05, num=500)
    pce = inst.get_total_eff(wave)    
    if normalize:
        pce /= pce.max()
    return (wave,pce)


def annotate_pce_plot(title="", xlim=(0.5,5.5), legend=True):
    plt.title(title)
    plt.ylabel("Photon Conversion Efficiency")
    plt.xlabel('Wavelength [$\mu$m]')
    plt.ylim(0,1)
    plt.xlim(*xlim)
    if legend:
        plt.legend(loc='upper right')
        
def plotlabelcolorfilt(wave,pce,filtname, color='blue'):
    plt.fill_between(wave,0,pce, alpha=0.5, color=color)
    plt.plot(wave,pce, color=color, label=filtname.upper())
    plt.plot(wave,pce, color='black', alpha=0.5)
    # double print text to make it a bit darker since some of the rainbow colors are hard to read
    plt.text(wave.mean(), 0.75, filtname.upper(), alpha=1, color='black', rotation=45)
    plt.text(wave.mean(), 0.75, filtname.upper(), alpha=0.75, color=color, rotation=45)

def uppercase_si_name(name):
    # the SI name with case as typically used, and specifically as used with WebbPSF. 
    name = name.upper()
    if name=='NIRCAM': return 'NIRCam'
    elif name=='NIRSPEC': return 'NIRSpec'
    else: return name
    

def output_to_fits_table(wave,pce,inst, verbose=True, filtername=None):

    data = [wave*1e4,pce/pce.max()]  # convert wave to angstroms as expected by pysynphot

    (year, month, day, hour, minute, second, weekday, DOY, DST) =  time.gmtime()

    si_name = uppercase_si_name(inst.inst_name)
    if filtername is None:
        filtername = inst.instrument['filter']
    filtername = filtername.upper()

    t = Table(data, names=['WAVELENGTH', 'THROUGHPUT'])
    t['WAVELENGTH'].unit='angstrom'
    t.meta['WAVEUNIT']=('angstrom', 'Wavelength unit')
    t.meta['TELESCOP']=('JWST', 'James Webb Space Telescope')
    t.meta['INSTRUME']=(si_name, 'Name of instrument')
    t.meta['FILTER']=(filtername, 'Bandpass filter identifier')
    t.meta['SOURCE']=('JWST ETC / Pandeia / REDCAT', 'Derived from JWST ETC database')
    t.meta["SRC_DATE"] = ( "%4d-%02d-%02d" % (year, month, day), "Date of file generation from ETC")

    t.meta['NORMALIZ']=(True, 'Normalized to peak transmission=1')

    t.meta['HISTORY'] = ['Normalized relative spectral response file for WebbPSF', '',
                         'Based on information extracted from the ',
                         'JWST Exposure Time Calculator (ETC) system',
                         '(Pandeia software database, using reference files',
                         'managed by ReDCaT, the Reference Data for ',
                         'Calibration and Tools group at STScI)', '',
                         'File created by: webbpsf/dev_tools/',
                         '  Normalized Throughputs from Pandeia.ipynb',
                         '']

    if not os.path.isdir('outputs'):
        os.mkdir('outputs')
    
    if not os.path.isdir(os.path.join('outputs',si_name)):
        os.mkdir(os.path.join('outputs',si_name))
    
    outname = os.path.join('outputs',si_name,"{}_throughput.fits".format(filtername))
    t.write(outname, overwrite=True)
    if verbose:
        print("==>> "+outname)
    return outname

# FGS is not in Pandeia - get the spectral response from someone on FGS team directly. 

from astropy.table import Table
import time
import webbpsf
_log = webbpsf.webbpsf_core._log
import pysynphot

def read_filter_data_like_webbpsf_does(filename, filtername):
    """Copied from webbpsf_core SpaceTelescopeInstrument._getSynphotBandpass() """
    # The existing FITS files all have wavelength in ANGSTROMS since that is
    # the pysynphot convention...
    
    #filterfits = fits.open(filter_info.filename)
    filterfits=fits.open(filename)
    
    waveunit = filterfits[1].header.get('WAVEUNIT')
    if waveunit is None:
        _log.warn('The supplied file, {}, does not have a WAVEUNIT keyword. Assuming it '
                  'is Angstroms.'.format(filter_info.filename))
        waveunit = 'angstrom'

    filterdata = filterfits[1].data
    try:
        band = pysynphot.spectrum.ArraySpectralElement(
            throughput=filterdata.THROUGHPUT, wave=filterdata.WAVELENGTH,
            waveunits=waveunit, name=filtername
        )
    except AttributeError:
        raise ValueError("The supplied file, %s, does not appear to be a FITS table "
                         "with WAVELENGTH and THROUGHPUT columns." % filter_info.filename)

    return band