Aperture Photometry with Mathematica

by Jeff Bryant

This notebook shows the steps involved in performing aperture photometry on a sequence of star images taken with a CCD camera.  Techniques required to process a full data set include centroid tracking of a moving object and isolation of specified regions.  Mathematica's support for the FITS image format makes this task much easier to accomplish.
   All of the images used here were collected in Flagstaff, Arizona by the author using a 31" telescope at the National Undergraduate Research Observatory (NURO) on Anderson Mesa.  This observatory is part of the Lowell Observatory complex.  The camera was a liquid nitrogen cooled CCD camera with a 498x511 pixel usable data area.  Exposures were all 60 seconds in length.  The images have all been reduced which only involved flat-fielding and bias frames.  No dark frame was used since the operating temperature kept the thermal noise to a minimum.  For a demonstration of how a CCD image can be reduced using Mathematica, visit the Wolfram Research website using the following URL.

http://library.wolfram.com/howtos/reduce/

Installation and Initialization

Once the CCDPhotometry package has been downloaded, place it into the AddOns\Applications folder in the Mathematica layout.  Once this is done, rebuild the help browser.  Now the documentation should be found under the AddOns section of the Help Browser.

First, we load the package

Needs["CCDPhotometry`"]

We can see all of the symbols defined in the package as follows.

Names["CCDPhotometry`*`*"]

Import the Data

data = 
      ImageLoad[
        ToFileName["CCDPhotometry", "11080154.fit"]]

Parsing the Header

The header of a FITS file can be obtained as follows.  Also, specific fields can be extracted.

Header[data]

HeaderField[data, "OBJECT"]

HeaderField[data, "FILTER"]

The following function extracts the time of the start of the exposure from the header.  The information is processed and converted to return the heliocentric Julian date of the observation.

ImageTime[data]

The right ascension and declination of the FITS object can be found as well.  This information is required for proper photometry calculations.

ImageRA[data]

ImageDEC[data]

Heliocentric Julian Date

The time recorded in the header is the time for the light from the object to reach the earth and be recorded.  Using the Gregorian calendar, it is difficult to plot data versus time, especially if the data spans over more than one day and definitely over a month.  The length of the time-axis would vary depending on the month.  The Julian calendar is more appropriate.  It expresses the date as the number of days since a fixed calibration day.  Time units less than a day are expressed as fractional days.

AccountingForm[JulianDate[ImageTime[data]], 16]

Although the Julian calendar is more linear, it still expresses the time required for the light to travel from the source to the earth.  Since the earth moves around the sun, its distance from the source object varies.  This varying distance can result in as much as (8^m) (13^s)time difference, depending on where the earth is in its orbit.  This can be a problem if high precision timings are needed.  To correct this time, the timings are calculated for a heliocentric coordinate system.  The Sun's distance doesn't change as much so this provides a more constant reference point.

The following function calculates the necessary correction time

HeliocentricCorrectionTime[ImageTime[data], ImageRA[data], 
      ImageDEC[data]]

When added to the Julian date, we have a linear timescale with a more stable reference point for our timing calculations

?HeliocentricJulianDate

AccountingForm[
      HeliocentricJulianDate[ImageTime[data], ImageRA[data], ImageDEC[data]], 
      16]

Visualize the Data

We need to know where the star of interest is located in a FITS image.  The following image can be used to locate the coordinates of the star in question. Once generated, click once on the image and hold the Control (or the Command key on a Mac or Shift-Mod on UNIX) key down as you move your mouse over the image.  You should see the mouse coordinates in the lower left corner of the notebook window.  These coordinates can be use to locate the coordinates of the star in question, a comparison star, and a tracking star.

DisplayImage[data, PlotRange -> {4, 75}];

Locate the Star in Question and the Comparison Star and Highlight Them

Using the image above, define the coordinates of the tracking star, the star in question, and a comparison star.

trackStar = {244, 183};

starOfInterest := {35, 70} + trackStar;

comparisonStar := {210, -29} + trackStar;

Performing the above operations on a single image is trivial.  When working with multiple images, a problem arises.  Telescopes rarely track perfectly.  As a result, the stars will typically move slightly in the image.  We need to be able to track the stars as they drift around.  The easiest methos is to provide an inital position, and then use a centroid tracking algorithm to update the stars coordinates.  If the stars position changes only a little, the following simple algorithm should work.  A function needs to be defined that allows the star to be centered within an aperture.  Given a best guess, or a previous location that was near the current position, this function should re-center the star in the aperture.

updateCoords[imgobj_, oldLoc_List] := 
      First[Position[imgobj[[1, 1, 2]], 
          Max[Flatten[IsolateRegion[imgobj, oldLoc, 20]]]]]

To correctly perform differential aperture photometry, we need a comparison star that does not vary in brightness.  Both stars will need to be referenced with respect to the tracked star.

DisplayImage[data, PlotRange -> {4, 75}, 
        Epilog -> {RGBColor[1, 1, 1], 
            Circle[Reverse@starOfInterest, 10], 
            Circle[Reverse@comparisonStar, 10], 
            Circle[Reverse@trackStar, 10], 
            Text["IP Peg", 
              Reverse@starOfInterest + {10, 0}, {-1, 0}], 
            Text["Tracking Star", 
              Reverse@trackStar + {10, 0}, {-1, 0}], 
            Text["comp. star", 
              Reverse@comparisonStar + {10, 0}, {-1, 0}]}];

Differential Aperture Photometry

Aperture photometry requires that all of the pixels within a specified radius be summed to get a measure of the brightness of the star in question.  The following function allows a circular region of the data around and containing the star to be extracted.

?Photometry

mystar = Photometry[data, starOfInterest, 10]

mycompstar = Photometry[data, comparisonStar, 10]

differentialMagnitude = 2.5 Log[mystar/mycompstar]

Phase[data, 2445922.96565, 0.158206424]

To construct a light curve of the star in question, a sequence of images of the star need to be collected an analyzed.  The code in this package can be used to build an algorithm to load a series of images, one at a time, and perform photometry on the image of the star in each frame.  This exercise requires a full set of images taken over a specific time span.  There are too many data files required to be able to realistically provide such an example here.  This task is left to the user.