Ground Based Photometry

Here I want to calculate some photometric points from spectra for comparison with published values for a bunch of known brown dwarfs.

In order to get the true magnitude for an object, I first need to calculate the instrumental magnitude and then correct for a number of effects. That is, I calculate the apparent magnitude from a particular place on the Earth and then add corrections to determine what it would be if we measured from space.

After all our corrections are made, the magnitude is given by:

Instrumental Magnitude

The first term on the right in the equation above is the magnitude measured by the instrument on the ground, given by:

Where is the energy flux density of the source in units of [erg s-1 cm-2 A-1] and is the scalar filter throughput for the band of interest.

Since I will be comparing my calculated magnitudes to photometry taken with photon counting devices, the factor of converts to a photon flux density in units [photons s-1 cm-2 A-1].

Zero Point Correction

The second term in our magnitude equation is a first order correction to compare to some standard we define as zero. I will use a flux calibrated spectrum of the A0 star Vega to calculate the zero point magnitude for the band:

Just as we obtained our instrumental magnitude above.

Extinction Correction

The third term is to correct for the extinction of the source flux due to atmospheric absorption. We can get closer to the true apparent magnitude (above the atmosphere) by adding an extinction term:

Where is the extinction coefficient for the band of interest and is the airmass.

The airmass is the optical path length of the atmosphere, which attenuates the source flux depending on its angle from the zenith . Approximating the truly spherical atmosphere as plane-parallel, the airmass goes from at to at . At zenith angles greater than that, the plane-parallel approximation falls apart and the airmass term gets complicated.

Where the airmass is the amount of atmosphere in the line of sight, the extinction coefficient is the amount by which the incident light is attenuated as it travels through the airmass. The extinction coefficient is related to the optical depth of the atmosphere as:

Where and are the magnitudes below and above the atmosphere respectively.

Example: J21512543-2441000

As an example, I'd like to calculate the 2MASS J-band magnitude of the brown dwarf at 21h51m25.43s -24d41m00s given a low resolution NIR energy flux density from the SpeX Prism instrument on the 3m NASA Infrared Telescope Facility.

J21512543-2441000

Interpolating the filter throughput to the object spectrum and then integrating as in the equation above, I get as my instrumental magnitude in the J-band.

Performing the same procedure on the flux calibrated spectrum of Vega, I get for my J-band zero point magnitude.

Checking the FITS file header, I will use for the airmass. The mean extinction coefficient for the MKO system J-band is given as in Tokunaga & Vacca (2007), making the atmospheric correction term .

The corrected magnitude is then:

Which is only 0.007 magnitudes off from the value of from the 2MASS catalog.

Remaining Problems

As shown in the example above, this works... but not for every object.

2MASS apparent J magnitudes vs. my calculated apparent j magnitudes for 67 brown dwarfs. The solid black line is for perfect agreement and the dashed line is a best fit of the data.

2MASS apparent J magnitudes vs. my calculated apparent j magnitudes for 67 brown dwarfs. The solid black line is for perfect agreement and the dashed line is a best fit of the data.

2MASS apparent H magnitudes vs. my calculated apparent h magnitudes for 67 brown dwarfs. The solid black line is for perfect agreement and the dashed line is a best fit of the data.

2MASS apparent H magnitudes vs. my calculated apparent h magnitudes for 67 brown dwarfs. The solid black line is for perfect agreement and the dashed line is a best fit of the data.

I whittled down my sample of 875 to only those objects with flux units and airmass values taken at Mauna Kea so that I could use the same extinction coefficient and make sure they are all in the same units of [erg s-1 cm-2 A-1].

Then I pulled the 2MASS catalog J and H magnitudes with uncertainties for these remaining objects and plotted them against my calculated values with uncertainties.

To the left are the plots of the 67 objects that fit the selection criteria in J-band (above) and H-band (below).

Though it's not the biggest sample, the deviation of the best fit line from unity suggests I'm off by a factor of 0.9 from the 2MASS catalog value across the board.

But more worrisome is the fact that most of the calculated magnitudes are not within the errors of the 2MASS magnitudes. This deviation ranges from very good agreement of a few thousandths of a magnitude up to the worst offenders of about 0.8 mags.

Filter Effective Wavelength(s)

The effective wavelength of a filter for narrow band photometry can easily be approximated by a constant and just looked up when needed. For broad band photometry, however, the width of the filter and the amount of flux in the band being measured actually come into play.

The effective wavelength of a filter is given by:

Where is the scalar filter throughput and is the flux density in units of [erg s-1 cm-2 A-1] or [photons s-1 cm-2 A-1] depending upon whether you are using an energy measuring or a photon counting detector, respectively.

Here are the results for 67 brown dwarfs with complete spectrum coverage of the 2MASS J-band:

Effective wavelength values for the 2MASS J-band filter. Blue and green circles indicate lambda calculated using photon flux densities (PFD) and energy flux densities (EFD) respectively. Filled circles are for 67 confirmed brown dwarfs. Open circles are for Vega.

Effective wavelength values for the 2MASS J-band filter. Blue and green circles indicate lambda calculated using photon flux densities (PFD) and energy flux densities (EFD) respectively. Filled circles are for 67 confirmed brown dwarfs. Open circles are for Vega.

The red line on the plot shows the specified value given by 2MASS. For fainter objects like brown dwarfs (filled circles), the calculated effective wavelength of the J-band filter can shift redward by as much as 150 angstroms. Vega (open circles) shifts it blueward by about 30 angstroms.

The difference is small but measurable and demonstrates the dependence of the filter width, source spectrum, and detector type on the effective wavelength while doing broad band photometry.

Brown Dwarf Synthetic Photometry

The goal here was to get the synthetic colors in the SDSS, 2MASS and WISE filters of ~2000 model objects generated by the PHOENIX stellar and planetary atmosphere software.

Since it would be silly (and incredibly slow... and much more boring) to just calculate and store every single color for all 12 filter profiles, I wrote a module to calculate colors a la carte.

The Filters

I got the J, H, and K band relative spectral response (RSR) curves in the 2MASS documentation, the u, g, r, i and z bands from the SDSS documentation, and the W1, W2, W3, and W4 bands from the WISE documentation.

I dumped all my .txt filter files into one directory and wrote a function to grab them all, pull out the wavelength and transmission values, and output the filter name in position [0], x-values in [1], and y-values in [2]:

def get_filters(filter_directory):
  import glob, os
  files = glob.glob(filter_directory+'*.txt')
 
  if len(files) == 0:
    print 'No filters in', filter_directory
  else:
    filter_names = [os.path.splitext(os.path.basename(i))[0] for i in files]
    RSR = [open(i) for i in files]
    filt_data = [filter(None,[map(float,i.split()) for i in j if not i.startswith('#')]) for j in RSR]
    for i in RSR: i.close()
 
    RSR_x = [[x[0] for x in i] for i in filt_data]
    RSR_y = [[y[1] for y in i] for i in filt_data]
    filters = {}
    for i,j,k in zip(filter_names,RSR_x,RSR_y):
      filters[i] = j, k, center(i)
 
    return filters

Calculating Apparent Magnitudes

We can't have colors without magnitudes so here's a function to grab the Teff and log g specified spectra, and calculate the apparent magnitudes in a particular band:

def mags(band, teff='', logg='', bin=1):
  from scipy.io.idl import readsav
  from collections import Counter
  from scipy import trapz, log10, interp
  
  s = readsav(path+'modelspeclowresdustywise.save')
  Fr, Wr = [i for i in s.modelspec['fsyn']], [i for i in s['wsyn']]
  Tr, Gr = [int(i) for i in s.modelspec['teff']], [round(i,1) for i in s.modelspec['logg']]
  
  # The band to compute
  RSR_x, RSR_y, lambda_eff = get_filters(path)[band]
  
  # Option to specify an effective temperature value
  if teff:
    t = [i for i, x in enumerate(s.modelspec['teff']) if x == teff]
    if len(t) == 0:
      print "No such effective temperature! Please choose from 1400K to 4500K in 50K increments or leave blank to select all."
  else:
    t = range(len(s.modelspec['teff']))
  
  # Option to specify a surfave gravity value
  if logg:
    g = [i for i, x in enumerate(s.modelspec['logg']) if x == logg]
    if len(g) == 0:
      print "No such surface gravity! Please choose from 3.0 to 6.0 in 0.1 increments or leave blank to select all."
  else:
    g = range(len(s.modelspec['logg']))
  
  # Pulls out objects that fit criteria above
  obj = list((Counter(t) & Counter(g)).elements())
  F = [Fr[i][::bin] for i in obj]
  T = [Tr[i] for i in obj]
  G = [Gr[i] for i in obj]
  W = Wr[::bin]
  
  # Interpolate to find new filter y-values
  I = interp(W,RSR_x,RSR_y,left=0,right=0)
  
  # Convolve the interpolated flux with each filter (FxR = RxF)
  FxR = [convolution(i,I) for i in F]
  
  # Integral of RSR curve over all lambda
  R0 = trapz(I,x=W)
  
  # Integrate to find the spectral flux density per unit wavelength [ergs][s-1][cm-2] then divide by R0 to get [erg][s-1][cm-2][cm-1]
  F_lambda = [trapz(y,x=W)/R0 for y in FxR]
  
  # Calculate apparent magnitude of each spectrum in each filter band
  Mags = [round(-2.5*log10(m/F_lambda_0(band)),3) for m in F_lambda]
  
  result = sorted(zip(Mags, T, G, F, I, FxR), key=itemgetter(1,2))
  result.insert(0,W)
  
  return result

Calculating Colors

Now we can calculate the colors. Next, I wrote a function to accept any two bands with options to specify a surface gravity and/or effective temperature as well as a bin size to cut down on computation. Here's the code:

def colors(first, second, teff='', logg='', bin=1):
  (Mags_a, T, G) = [[i[j] for i in get_mags(first, teff=teff, logg=logg, bin=bin)[1:]] for j in range(3)]
  Mags_b = [i[0] for i in get_mags(second, teff=teff, logg=logg, bin=bin)[1:]]
  colors = [round(a-b,3) for a,b in zip(Mags_a,Mags_b)]
  
  print_mags(first, colors, T, G, second=second)
  
  return [colors, T, G]

The PHOENIX code gives the flux as Fλ in cgs units [erg][s-1][cm-2][cm-1] but as long as both spectra are in the same units the colors will be the same.

Makin' It Handsome

Then I wrote a short function to print out the magnitudes or colors in the Terminal:

def print_mags(first, Mags, T, G, second=''):
  LAYOUT = "{!s:10} {!s:10} {!s:25}"
  
  if second:
    print LAYOUT.format("Teff", "log g", first+'-'+second)
  else:
    print LAYOUT.format("Teff", "log g", first)
  
  for i,j,k in sorted(zip(T, G, Mags)):
    print LAYOUT.format(i, j, k)

The Output

Then if I just want the J-K color for objects with log g = 4.0 over the entire range of effective temperatures, I launch ipython and just do:

In [1]: import syn_phot as s
In [2]: s.colors('J','K', logg=4)
Teff -------- log g -------- J-K
1400.0 ------ 4.0 ---------- 4.386
1450.0 ------ 4.0 ---------- 4.154
...
4450.0 ------ 4.0 ---------- 0.756
4500.0 ------ 4.0 ---------- 0.733

Similarly, I can specify just the target effective temperature and get the whole range of surface gravities. Or I can specify an effective temperature AND a specific gravity to get the color of just that one object with:

In [3]: s.colors('i','W2', teff=3050, logg=5)
Teff -------- log g -------- J-K
3050.0 ------ 5.0 ---------- 3.442

I can also reduce the number of data points in each flux array if my sample is very large. I just have to specify the number of data points to skip with the "bin" optional parameter. For example:

In [4]: s.colors('W1','W2', teff=1850, bin=3)

This will calculate the W1-W2 color for all the objects with Teff = 1850K and all gravities, but only take every third flux value.

I also wrote functions to generate color-color, color-parameter and color-magnitude plots but those will be in a different post.

Plots!

Here are a few color-parameter animated plots I made using my code. Here's how I made them. Click to animate!

And here are a few colorful-colorful color-color plots I made:

Plots with observational data

Just to be sure I'm on base, here's a color-color plot of J-H vs. H-Ks for objects with a log surface gravity of 5 dex (blue dots) plotted over some data for the Chamaeleon I Molecular Cloud (semi-transparent) from Carpenter et al. (2002).

The color scale is for main sequence stars and the black dots are probable members of the group. Cooler dwarfs move up and to the right.

And here's a plot of J-Ks vs. z-Ks as well as J-Ks vs. z-J. Again, the blue dots are from my synthetic photometry code at log(g)=5 and the semi-transparent points with errors are from Dahn et al. (2002).