spacetelescope / roman_tools

Instructions and tutorials for Roman software tools distributed by STScI for the science community
BSD 3-Clause "New" or "Revised" License
14 stars 6 forks source link

Getting the WFIRST filters using Pandeia #38

Open rmjarvis opened 4 years ago

rmjarvis commented 4 years ago

I was trying to reproduce the WFIRST filter throughputs that were posted on the NASA website. And I think I've managed to get close, but there seems to be something missing in my calculation, so I was hoping someone could let me know what I'm doing wrong.

from urllib.request import urlopen
import os
import numpy as np
import xlrd  # To read the xlsm file

from pandeia.engine.calc_utils import get_instrument_config
from pandeia.engine.instrument_factory import InstrumentFactory

# Get the filter spread sheet from nasa.gov.
file_name = 'WFIRST_WIMWSM_throughput_data_190531.xlsm'
if not os.path.isfile(file_name):
    # Linked from here: https://roman.gsfc.nasa.gov/science/Roman_Reference_Information.html
    url = 'https://roman.gsfc.nasa.gov/science/201907/' + file_name
    with urlopen(url) as u:
        with open(file_name, 'wb') as f:
            f.write(u.read())

# Read the spreadsheet
workbook = xlrd.open_workbook(file_name)
sheet = workbook.sheet_by_name('EffectiveArea')
wave = sheet.col_values(0,19)
filters = sheet.row_values(18,1,8)
print('filters = ',filters)

# Set up the Pandeia configuration for the wfirstimager
config = get_instrument_config('wfirst','wfirstimager')
config = config['defaults']['imaging']

# Get each filter and compare to value in spreadsheet
for i, filt in enumerate(filters):
    print('filter name = ',filt)

    # Get the throughputs with Pandeia:
    config['instrument']['filter'] = filt.lower()
    inst = InstrumentFactory(config)
    thru = inst.get_total_eff(wave)
    # The spreadsheet has this in terms of effective area, so multiply by coll_area
    # * 1.e-4 since coll_area is apparently in cm^2, not m^2.
    area1 = thru * inst.telescope.coll_area * 1.e-4

    # Compare this to the values in the NASA spreadsheet
    area2 = np.array(sheet.col_values(i+1,19))

    # Should be close to equal (?)
    nz = np.nonzero(area2)
    print('mean difference = ',np.mean(area1[nz] - area2[nz]))
    print('median ratio = ',np.median(area1[nz] / area2[nz]))
    big = np.where(area2 > 0.5)  # ignore where throughput is small.
    print('max ratio = ',np.max(area1[big] / area2[big]))
    # In fact most of the non-trivial values have a very similar ratio. 
    # (Mostly 0.83333 except F184)
    #print(area1[big] / area2[big])

The output when I run this is

$ python getRomanFilters.py 
filters =  ['F062', 'F087', 'F106', 'F129', 'F146', 'F158', 'F184']
filter name =  F062
mean difference =  -0.03379910524097619
median ratio =  0.7649984352572698
max ratio =  0.8333311192673077
filter name =  F087
mean difference =  -0.0381553731571408
median ratio =  0.7583318484603302
max ratio =  0.833330978298646
filter name =  F106
mean difference =  -0.046309478066203065
median ratio =  0.7545554915913855
max ratio =  0.8333310166506992
filter name =  F129
mean difference =  -0.05516289690213783
median ratio =  0.7708317019566648
max ratio =  0.8333309982462173
filter name =  F146
mean difference =  -0.18500388996636585
median ratio =  0.8041645475199117
max ratio =  0.8333310166506988
filter name =  F158
mean difference =  -0.06898710478670267
median ratio =  0.7666706348918275
max ratio =  0.8333309964445996
filter name =  F184
mean difference =  -0.01954905812563244
median ratio =  0.8502927112595088
max ratio =  0.9395410391489472

So I'm clearly close, but it doesn't quite match. Am I doing something wrong here? Or should these not match for some reason? Thanks much in advance!

ariedel commented 4 years ago

I haven't looked at your numbers in detail, but the values in Pandeia are actually from spreadsheets David Hughes at Goddard sent to me. The values on the site were in effective area (which wasn't immediately useful to me), and I later found out that the filter throughputs posted on that site were filter throughputs PLUS detector quantum efficiency and optical absorption baked in. The ones in Pandeia (from his spreadsheets) are JUST the throughputs.

Therefore, there's a reason your filters wouldn't match the ones in Pandeia.

rmjarvis commented 4 years ago

I thought that was the Instrument.get_filter_eff function. I used get_total_eff which multiplies these things together, so I thought that should match the values in the spread sheet. Certainly, just using get_filter_eff returns something that doesn't look anything much like the values in the spread sheet. Only when multiplying by QE (get_detector_qe) and optical (get_internal_eff) does it come close.

rmjarvis commented 4 years ago

In case anyone has some time to look at this, the main thing I want to confirm from this process is what are the canonical Roman filter curves, and what exactly do they include? I need this to update the GalSim filter curve files, which are currently out of date.

I had thought that Pandeia would have the bespoke filter curves, since they seem to have things specified at a higher level of detail (separating QE from filter from optical throughput), and my script was mostly to confirm that I understood how all that worked, expecting to get the curves on the website, which combine all of these. But since they don't agree, I mostly want to know, am I doing something wrong here? Or are they not expected to agree, since either Pandeia or the website is not actually up to date, and I should just trust the other one.

Thanks for any advice you can offer.

rmjarvis commented 4 years ago

OK, trimming my expectations a second time, since I'm not getting much of a response so far.

Can someone just tell me which set of filter information is expected to be more up to date? Pandeia or the gsfc web site?