How to compute a Flux through a filter using SYNPHOT in Python

Question

I am in the process of building an exposure time calculator to help nail down the exposure times I need for different targets in my observatory.

In amateur astronomy imaging, we use a filter called Luminance filter, that lets all the light from a source pass through to the CCD sensor, from 400 nm to 700 nm. It is basically a very wide broadband filter. Online I found the fluxes for UBVRI filters but not for such a broadband filter, probably because it doesn't have any use for scientific work. For example with a V filter, mag=0 star (Vega), the Flux = 3.75e-9 ergs/sec/cm2/A. I would like to find the Flux for the Luminance filter for mag=0 star

Recently I came across the SYNPHOT module for python, so I tried to simulate a filter, but I don't know how to actually compute the Flux for such a filter. Here is the python code:

import numpy as np
from synphot.models import BlackBody1D,BlackBodyNorm1D
from synphot import units, SourceSpectrum, SpectralElement
import matplotlib.pyplot as plt

# create a source spectrum type like Vega (mag = 0)
black_body = SourceSpectrum(BlackBodyNorm1D, temperature=10000)

#create the filter with a central waveL of 5500 A and a width of 2500 A on either side
L_filter = SpectralElement(Box1D, amplitude=1, x_0=5500, width=2500)

#plot the filter waveL vs transmission, convert A to nm
L_filter.plot(wavelengths=bp.waveset.to(u.nm), top=1.1, title='Luminance filter')
plt.ylabel("Transmission")

#Normalize the source spectrum to 1 Jy in a given box bandpass and integrate it:
sp_rn = black_body.normalize(1 * u.Jy, band=L_filter)
print(sp_rn.integrate())

#Create an observation by passing the normalized source spectrum through the box bandpass:
obs = Observation(sp_rn, L_filter)

Any idea if this is the way to do it? And how should I continue to compute the Flux?

Answer 1

For the objective of setting exposure times for different targets, it may be easiest to just use the V-filter fluxes you have for the targets, or take the average of the B-, V-, and R-filter target fluxes as a rough guess at the Luminance-filter target flux.

To support this recommendation, I'll try to answer the title question using SYNPHOT.

But first it's important to define what "flux" means. The example in the question has units of erg s^-1 cm^-2 Å^-1. These are units for spectral irradiance (or spectral flux density). If this is integrated over wavelength, it becomes irradiance (or flux density), with units of erg s^-1 cm^-2. The jansky is a unit of spectral flux density, but it is per unit frequency (Hz^-1) instead of per unit wavelength (Å^-1). Here is a table of some of these units that are available in SYNPHOT:

First compare the example source to the filter passbands mentioned in the question. We use a 10000-K blackbody as the source spectrum for Vega, but SYNPHOT contains a Vega spectrum that could be used instead. SYNPHOT also contains predefined filter transmission functions for B, V, and R. We build the Luminance filter as a boxcar function:

The units in the plot are the same spectral flux density units as in the example from the original question, with filter transmissions arbitrarily scaled to show up on the graph. The spectral flux density of the source actually changes by over a factor of 2 across the Luminance filter (or from B to R). So it's worth understanding what the 3.75e-9 erg s^-1 cm^-2 Å^-1 in the original question refers to. In the V filter, the source has a spectral flux density of 3.32e-9 erg s^-1 cm^-2 Å^-1 at a wavelength of 5479 Å, which is close enough to the example that it could be explained by differences between the actual Vega spectrum and the 10000-K blackbody. The wavelength of 5479 Å is the pivot wavelength, which is one way of defining the central wavelength or effective wavelength of the filter. I chose pivot wavelength because it's simple to get it from SYNPHOT.

SYNPHOT can be used to integrate the spectral flux density over filter bandpasses to get flux density, to determine average spectral flux density over filter bandpasses, or to sample the spectral flux density at the effective wavelengths of the filter bandpasses. Here is a table comparing the results:

          Integrated          Averaged              Sampled 
        (erg cm-2 s-1)   (erg cm-2 s-1 Å-1)   (erg cm-2 s-1 Å-1)
--------------------------------------------------------------------
B:         5.26E-06           2.63E-09              5.23E-09
R:         4.61E-06           1.05E-09              2.14E-09
V:         3.24E-06           1.41E-09              3.32E-09
B+V+R:     1.31E-05           5.09E-09              3.56E-09
Luminance: 9.69E-06           3.88E-09              3.80E-09

Actually calculating all this requires the source spectrum. If your input data only consists of broadband flux values (in B, V, R), then you don't have the necessary spectral details. But for this example, we can compare different methods of combining the broadband filter data.

The integrated flux density may be best for an exposure time calculator (but you may want to play with the units to get photons s^-1 cm^-2 instead). For this example, just adding up the B+V+R flux densities comes pretty close to the calculated Luminance flux density. But, without knowing more about the spectra of your targets, you won't really be able to integrate the spectral flux densities.

The averaged spectral flux density column again cannot be calculated without more info on the target spectra, and it's included here mainly to show that it's not useful. These values are basically the area under the curves of the spectra where the source and filter have been convolved. Here, adding B+V+R comes somewhat close to the Luminance average spectral flux density, but it's not very good.

The best result is to use the "Sampled" column, and take the average of B, V, and R spectral flux densities (which is what you have in the example for Vega in V-band). For our example blackbody, the average only has a 7% error. Second best for this example would be to just use the V-band spectral flux density, with a 13% error. The actual errors will vary for each target with a different spectrum.

---

Below is the code to make the plot and the table. For the table, I used a SYNPHOT object called an Observation: "Observation is a special type of Source Spectrum, where the source is convolved with a Bandpass." I imported Astropy units as "u" and SYNPHOT units as "units" to keep them separate, but I am not fully sure how this works. It could be redundant if the SYNPHOT units already contain all the Astropy units. I should disclose that I am definitely not a SYNPHOT (or Python) expert, so others should feel free to edit the code and text if I have made any errors. But the numbers check out pretty well for Vega, so I am confident of the conclusion that it may be easiest to just use the V-filter fluxes you have for the targets, or take the average of the B-, V-, and R-filter target fluxes as a rough guess at the Luminance-filter target flux.

import numpy as np
from astropy import units as u
from synphot.models import BlackBody1D,BlackBodyNorm1D,Box1D
from synphot import units, SourceSpectrum, SpectralElement, Observation
import matplotlib.pyplot as plt

# define our units for later (use astropy version because SYNPHOT's "FLAM" is not descriptive)
flxden = units.FLAM 
flxden = (u.erg / u.cm**2 / u.s / u.AA)

# create Luminance filter with a central waveL of 5500 A and a width of 2500 A on either side
L_filter = SpectralElement(Box1D, amplitude=1, x_0=5500, width=2500)
# also use SYNPHOT predefined filters
b = SpectralElement.from_filter('johnson_b')
v = SpectralElement.from_filter('johnson_v')
print("V-band pivot wavelength: ", v.pivot())
r = SpectralElement.from_filter('johnson_r')

# create a source spectrum type like Vega (mag = 0)
black_body = SourceSpectrum(BlackBodyNorm1D, temperature=10000)
# normalize just scales the spectrum to the given value in the given bandpass
bbspec_v = black_body.normalize(3.75e-9*flxden, band=v)

plt.figure(figsize=(8, 5.5))
plt.plot(black_body.waveset, bbspec_v(wavelengths=black_body.waveset, flux_unit=flxden), 'black')
plt.plot(b.waveset, b(b.waveset)*9e-9*flxden, 'b', v.waveset, v(v.waveset)*9e-9*flxden, 'g', r.waveset, r(r.waveset)*9e-9*flxden, 'r', L_filter.waveset, L_filter(L_filter.waveset)*9e-9*flxden, 'gold')  
plt.ylim(0, 1e-8)  
plt.xlim(3000, 9000)
plt.legend(['Vega (10000 K blackbody)', b.meta['header']['descrip'], v.meta['header']['descrip'], r.meta['header']['descrip'], 'Luminance'])
plt.title('Filter transmission curves')
plt.ylabel('Spectral flux density (erg cm-2 s-1 Å-1)') 
plt.xlabel('Wavelength (Å)')

# make table output

# blackbody flux densities
bb_v_int = Observation(bbspec_v, v).integrate(flux_unit=flxden)
bb_b_int = Observation(bbspec_v, b).integrate(flux_unit=flxden)
bb_r_int = Observation(bbspec_v, r).integrate(flux_unit=flxden)
bb_luminance_int = Observation(bbspec_v, L_filter).integrate(flux_unit=flxden)

# blackbody bandpass-average spectral flux densities
bb_v_avg = Observation(bbspec_v, v).integrate(flux_unit=flxden)/(v.waverange[1] - v.waverange[0])
bb_b_avg = Observation(bbspec_v, b).integrate(flux_unit=flxden)/(b.waverange[1] - b.waverange[0])
bb_r_avg = Observation(bbspec_v, r).integrate(flux_unit=flxden)/(r.waverange[1] - r.waverange[0])
bb_luminance_avg = Observation(bbspec_v, L_filter).integrate(flux_unit=flxden)/(L_filter.waverange[1] - L_filter.waverange[0])

# blackbody bandpass center spectral flux densities
bb_v_pivot = Observation(bbspec_v, v)(v.pivot(),flux_unit=flxden)
bb_b_pivot = Observation(bbspec_v, b)(b.pivot(),flux_unit=flxden)
bb_r_pivot = Observation(bbspec_v, r)(r.pivot(),flux_unit=flxden)
bb_luminance_pivot = Observation(bbspec_v, L_filter)(L_filter.pivot(),flux_unit=flxden)

# calculate B+V+R sums and average
sum_int = bb_b_int.value+bb_v_int.value+bb_r_int.value
sum_avg = bb_b_avg.value+bb_v_avg.value+bb_r_avg.value
avg_pivot = (bb_b_pivot.value+bb_v_pivot.value+bb_r_pivot.value)/3


print("bb_v_int.unit: ", bb_v_int.unit)
print("bb_v_avg.unit: ", bb_v_avg.unit)
print("bb_v_pivot.unit: ", bb_v_pivot.unit)
print(" ")


print("          Integrated          Averaged              Sampled ")
print("        (erg cm-2 s-1)   (erg cm-2 s-1 Å-1)   (erg cm-2 s-1 Å-1)")
print("--------------------------------------------------------------------")
print("B:         %7.2E           %7.2E              %7.2E" % (bb_b_int.value, bb_b_avg.value, bb_b_pivot.value))
print("R:         %7.2E           %7.2E              %7.2E" % (bb_r_int.value, bb_r_avg.value, bb_r_pivot.value))
print("V:         %7.2E           %7.2E              %7.2E" % (bb_v_int.value, bb_v_avg.value, bb_v_pivot.value))
print("B+V+R:     %7.2E           %7.2E              %7.2E" % (sum_int, sum_avg, avg_pivot))
print("Luminance: %7.2E           %7.2E              %7.2E" % (bb_luminance_int.value, bb_luminance_avg.value, bb_luminance_pivot.value))

Package version info:

Python 3.7.6
astropy                   4.0.2
matplotlib                3.3.2
numpy                     1.19.1
synphot                   1.1.1

Answer 2

Last night after reading more Synphot documentation I came up with this method:

bb = SourceSpectrum.from_vega()

#create the filter with a central waveL of 5500 A and a width of 3000 A on either side

L_filter = SpectralElement(Box1D, amplitude=1, x_0=5500, width=3000)

#plot the filter waveL vs transmission, convert A to nm

L_filter.plot(wavelengths=bb.waveset.to(u.nm), top=1.1, title='Luminance filter',left=300, right=1000)
plt.ylabel("Transmission")

#Create an observation by passing the source spectrum through the box bandpass:
#FLAM = erg/sec/cm^2/A

obs = Observation(bb, L_filter)
y=obs.effstim(flux_unit='flam') 
print(y)

# The result is: **3.716e-09 FLAM**