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

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?

• It's not 100% clear what kind of fluxes you are looking for. Where you say "Online I found the fluxes for UBVRI filters," could you edit your question to give an example of a target and the flux for one of these filters? That would help to show the units of "Flux." For example, if you have BVR fluxes for a target, then your Luminance flux could be roughly approximated by B+V+R, or by (B+V+R)/3, depending on whether units are integrated over wavelength/frequency or not. This might be good enough for estimating exposure times. Dec 30, 2021 at 7:45
• Thank you, I added an example for the V filter with the units Dec 30, 2021 at 7:55

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.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

• Thank very much for the detailed answer. I will go over the code later as now I am at work. Last night I found out that indeed synphot has the Vega spectrum included. Then instead of the blackbody of 10000K I used the included source spectrum of Vega, and I was able to compute the Flux density for the Luminance filter using the effstim() function and the result is very close to your method. I will post the code. Dec 31, 2021 at 8:25

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**