Stellar Aperture Photometry

Thomas Sherlock, Wolfram Research

ABSTRACT

Stellar CCD aperture photometry is the technique of extracting information about the brightness of stars from a series of images collected over time.  The light curve of a variable star can reveal useful information about the physics of the star, including a measure of the star’s intrinsic brightness.  Light curve analysis can yield information about eclipsing binary systems and also lead to exoplanet discoveries where a planet alters the brightness of a star by crossing its disk as viewed from Earth.  

The Wolfram Language is an ideal tool for extracting and analyzing this data because it combines image processing, data reduction, plotting, curve fitting, file import and astronomical data into one package.  This talk will address the basic techniques of stellar aperture CCD photometry using the Wolfram Language.

Introduction

In CCD photometry, we want to be able to determine a measure of the amount of radiation coming from a given star arriving on our CCD detector.  As plotted as a function of time, this can reveal important information about the star or star system.

Cepheid Variables

With Cepheid variable stars, there is a strong correlation between the star’s light curve and the star’s intrinsic brightness, or absolute magnitude.  This period-luminosity relationship, first published in 1912, is called Leavitt’s Law, after Henrietta Swan Leavitt, based on her work at the Harvard College Observatory:

Leavitt worked as a ‘human computer’ at the observatory, manually extracting observation data from glass plates, for which she was paid $0.30 per hour (eventually).

For Population I Cepheid’s, Leavitt’s Law relates the log of the star’s period (in days) to the mean absolute magnitude:

absoluteMagnitude[P_]:=-2.81*Log10[P]-1.43

absoluteMagnitude[3.0]

LogLinearPlot[absoluteMagnitude[p],{p,.1,100},PlotTheme->"Scientific",GridLines->Automatic, FrameLabel->{"Period (days)", "Absolute Magnitude"},ImageSize->800]

Since the observed apparent magnitude diminishes the further the star is from Earth, the distance (in parsecs) can be determined from the period and the apparent magnitude.

distanceInParsecs[AbsoluteMagnitude_,ApparentMagnitude_]:=Power[10.0,(ApparentMagnitude-AbsoluteMagnitude+5)/5];

distanceInParsecsFromPeriod[P_,ApparentMagnitude_]:=Power[10.0,(ApparentMagnitude-absoluteMagnitude[P]+5)/5];

Hubble’s Discovery

In 1924, Edwin Hubble identified a Cepheid variable in what was then know as the ‘Andromeda Nebula’.

By extracting its absolute magnitude from the Cepheid's light curve, and comparing it to the observed magnitude, he was able to show the distance to the star was millions of light years and prove that it was part of an entirely separate galaxy thus vastly expanding the size of the known universe.

Assumptions

Image frames must be calibrated

Dark Frames Subtracted (Removes Thermal Noise)

Flat Frames Subtracted (Removes Illumination Variation)

Image Frames Need to be Aligned (allows photometric comparison across multple images)

These allow us to compare different stars across a CCD field and get accurate magnitude comparisons.  Bias noise associated with reading the CCD array is automatically subtracted out by the main photometry computation.

Technique

Divide Area Around Star into 3 Regions

Star Aperture isolates star photons

Sky Annulus isolates sky background photons

Gap isolates the star aperture from the sky annulus

Use Mathematica to generate disk and annular masks to isolate regions on the raster

starAperture[im_,starDiameter_]:=Module[{id=ImageDimensions[im]},Rasterize[Graphics[{Black,Rectangle[{0,0},id],White,Disk[{id[[1]]/2,id[[2]]/2},starDiameter/2]},PlotRangePadding->None,ImagePadding->None],"Image",RasterSize->id]];

skyAnnulus[im_,innerDiameter_]:=Module[{id=ImageDimensions[im]},Rasterize[Graphics[{Black,Rectangle[{0,0},id],White,Annulus[{id[[1]]/2,id[[2]]/2},{innerDiameter/2,id[[1]]/2}]},PlotRangePadding->None,ImagePadding->None],"Image",RasterSize->id]];

Example using Peppers

peppers=ExampleData[{"TestImage","Peppers"}]

ImageMultiply[peppers,starAperture[peppers,200]]

ImageMultiply[peppers,skyAnnulus[peppers,300]]

Compute Total Star Brightness Minus Sky Background

Isolate the star of interest from the full CCD frame:

isolateRegion[im_Image,{row_,column_},radius_]:=ImageTake[im,{row-radius,row+radius},{column-radius,column+radius}];

Total signal coming from the star can then be computed by subtracting out the sky background

Photometry[im_Image,{row_,column_},starApertureDiameter_,skyAnnulusInnerDiameter_,skyAnnulusOuterDiameter_]:=Module[{baseImage,starImage,skyImage,skyBackgroundPerPixel,starTotal,skyTotal,skyAnnulusArea,starApertureArea},

(*Carve out area just around the star*)baseImage=isolateRegion[im,{row,column},skyAnnulusOuterDiameter/2.0];

(*Extract the star image*)
starImage=ImageMultiply[baseImage,starAperture[baseImage,starApertureDiameter]];

(*Extract the pixels in the sky annulus*)
skyImage=ImageMultiply[baseImage,skyAnnulus[baseImage,skyAnnulusInnerDiameter]];

(*Total signal from star*)
starTotal=Total@Flatten@ImageData@starImage;

(*Total signal from sky annulus*)
skyTotal=Total@Flatten@ImageData@skyImage;
skyAnnulusArea=Pi*(skyAnnulusOuterDiameter^2-skyAnnulusInnerDiameter^2)/4.0;

(*Compute the per pixel sky signal*)starApertureArea=Pi*(starApertureDiameter^2)/4.0;
skyBackgroundPerPixel=skyTotal/skyAnnulusArea;

(*Subtract out the total sky background from the star image*)starTotal-starApertureArea*skyBackgroundPerPixel
];

Determining the Time Image Was Taken

Use Julian Dates:

A Julian date is a real number representing the number of days which have transpired since noon Universal Time on January 1, 4713 BC (on the Julian calendar).  This is a good way to represent an observation time.

Bundle FITS Image and Metadata Together

LoadFITSObject[imagePath_String]:=Association[{Image->Import[imagePath,"FITS"],MetaInformation->Import[imagePath,{"FITS", "MetaInformation"}]}];

Extract Julian Date From FITS Object

Unfortunately there are several ways that the observation date is represented in FITS files.

GetJulianDateAlt[fitsObject_]:=
Module[{date, time,do, jd},
date=StringSplit[fitsObject[MetaInformation][1]["GENERAL INFORMATION"]["DATE-OBS"][[1]],{"-","/"}];
time=StringSplit[fitsObject[MetaInformation][1]["GENERAL INFORMATION"]["UT"][[1]],":"];
If[time=={"KeyAbsent"},time=StringSplit[fitsObject[MetaInformation][1]["GENERAL INFORMATION"]["TIME-OBS"][[1]],":"]];
do=DateObject[ToExpression/@{date[[3]],date[[1]],date[[2]],time[[1]],time[[2]],time[[3]]}];
JulianDate[do]
];

GetJulianDate[fitsObject_]:=
Module[{jd},
jd=fitsObject[MetaInformation][1]["GENERAL INFORMATION"]["JD"][[1]];
If[jd=="KeyAbsent",jd=GetJulianDateAlt[fitsObject]];
jd
];

Extract Brightness as a function of time

This assumes that the comparison star has a known magnitude and isn’t a variable star as well.

PhotometryObs[fitsObject_, star_List, compstar_List,compStarMag_, starDiam_,gapDiam_, skyDiam_]:=
Module[{imgobj,jd,ms, mcs, diffMag},
(* Get the image from the FITS object *)
imgobj=fitsObject[Image][1];
(* Get the Julian date of the obs time from FITS object *)
jd=GetJulianDate[fitsObject];
(* Target star brightness *)
ms=Photometry[imgobj,star,starDiam,gapDiam, skyDiam];
(* Reference star brightness *)
mcs=Photometry[imgobj,compstar,starDiam,gapDiam, skyDiam];
(* Compute the magnitude difference *)
diffMag=-2.5 Log10[ms/mcs];
(* Return the {julian date, magnitude difference pair, offset to the mag of the reference star *)
{jd,diffMag + compStarMag}
];

Example: Cepheid Variable Light Curve

Note : this is a synthetic star field.

Data Source: ftp://deimos.physics.uiowa.edu/labimage/archive/GeneralAstronomy/Cepheids/

Input FITS files

fileNames=FileNames[FileNameJoin[{NotebookDirectory[],"Data","Cepheids","*.fts"}]];

Load ‘bundled’ FITS objects

cepheidImages=LoadFITSObject/@fileNames;

Star chart provided with this set

Import[FileNameJoin[{NotebookDirectory[],"Data","Cepheids","S366Chart.jpg"}]]

Load the first image, so we can locate the cepheid on the image.

ImageAdjust[cepheidImages[[1]][Image][1]]

Once we have the {row, column} coordinates of the target star and a reference star on the image, extract the photometry data as a function of time, adjusted to the magnitude of the reference star (here taken to be 10.0)

Here the star diameter is 6 pixels, the outer gap diameter is 8 pixels and the sky annulus outer diameter is 20.

sortedData=SortBy[PhotometryObs[#,{513,513},{533,289},10.0,6,8,20]&/@cepheidImages, First];

First data point showing {Julian Date, Apparent Magnitude}

sortedData[[1]]

Plot the light curve starting from the earliest observation:

startJD=sortedData[[1]][[1]];

offsetSortedData={#[[1]]-startJD,#[[2]]}&/@sortedData;

ListLinePlot[offsetSortedData,PlotTheme->"Scientific"]

Analyzing Data

Try a simple fit to a sin() function

model=Association[FindFit[offsetSortedData,(Subscript[M, v])+A Sin[ω t + ϕ],{(Subscript[M, v]), A, ω, ϕ},t]]

Here, (Subscript[M, v]) is the apparent visual magnitude.  Unfortunately, FindFit isn’t the most accurate way to determine the frequency parameter ω when fitting to periodic data.

Fast Fourier Transforms (FFTs) require data points which are regularly spaced.  Photometric data isn’t usually of this type because observations may be spaced over days, or weeks or longer.

IrregularPeriodogram

The best way to analyze data like this is to use IrregularPeriodogram, which is available on our Function Repository.  The math behind IrregularPeriodogram grew out of the needs the astronomical community which frequently needs to analyze data from non-regularly spaced observations.

observationTimes=offsetSortedData[[All,1]];

observationValues=offsetSortedData[[All,2]];

IrregularPeriodogram needs the data to be zero offset.

observationValues2=observationValues-Mean[observationValues];

ResourceFunction["IrregularPeriodogram"]

Plot[ResourceFunction["IrregularPeriodogram"][w,observationTimes,observationValues2],{w,.1,1},PlotPoints->300,PlotRange->All,AxesLabel->{"ω (radians)","\!\(\*SubscriptBox[\(P\), \(D\)]\)(ω)"},ImageSize->800]

This shows a strong peak near .4 radians which corresponds to the fundamental angular frequency of the data.  Numerically, this frequency can be determined by NArgMax:

freq=NArgMax[{ResourceFunction["IrregularPeriodogram"][w,observationTimes,observationValues2],0≤w≤10},w]

The period which corresponds to the lowest order frequency is the period of the variable star, in days:

Period=2π/freq

Then using our distance-period relationship and the apparent visual magnitude of the star, we can compute the distance, in parsecs:

distanceInParsecsFromPeriod[Period,model[(Subscript[M, v])]]

Further Work

Real star fields tend to swim around a bit

Example (IP Pegasi)

Data courtesy of J. Bryant.

IP Pegasi. Its a U Geminorum type eclipsing cataclysmic variable star.  The light you are seeing if from a glowing accretion disk around a white dwarf star (1.16 solar masses) that gets periodically eclipsed by a 
cooler and much fainter secondary star that is a small type M star (~.56 solar masses). Mass is being siphoned off of the M star and onto the white dwarf via the accretion disk and mass transfer stream. Both stars 
are on the order of a solar mass or less.

fileNames2=FileNames[FileNameJoin[{NotebookDirectory[],"Data","IPPegasi","*.fit"}]];

ippegasiImages=LoadFITSObject/@fileNames2;

goodippegasiImages=Take[ippegasiImages,{228,317}];

Manipulate[ImageAdjust[goodippegasiImages[[i]][Image][1],{-.1,-.98}],{i,1,Length[goodippegasiImages],1}]

Aligning Star Fields

AlignStarFieldSimple[frames_List]:=
Module[{trackedStars,width,height,geomTransforms, scaledFrames, croppedFrames},
(*If I omit the conversion to Bit16 here,ImageFeatureTrack will fail to find any features.*)scaledFrames=Map[Image[#,"Bit16"]&,frames];
width=Min[(ImageDimensions/@scaledFrames)[[All,1]]];
height = Min[(ImageDimensions/@scaledFrames)[[All,2]]];
Print["cropping frames to "<>ToString@width <> " by "<>ToString@height];
croppedFrames=ImageCrop[#,{width,height}]&/@scaledFrames;
Print["aligning "<>ToString@Length[croppedFrames]<>" cropped frames"];
(*For each frame,create a list of stars that can be tracked across the set*)trackedStars=ImageFeatureTrack[croppedFrames,MaxFeatures->3];
Print["Tracking "<>ToString@Length[trackedStars[[1]]]<>" stars"];
If[Length[trackedStars[[1]]] < 2,Print["Found fewer than 2 tracked stars!"];Return[frames]];
(*Find the set of geometric transforms indicated by our set of good stars tracked across all good frames*)geomTransforms=FindGeometricTransform[trackedStars[[1]],#]&/@trackedStars;
(*Apply the transformations to our good frames to get them into alignment*)MapThread[ImagePerspectiveTransformation[#1,#2,DataRange->Full]&,{croppedFrames,geomTransforms[[All,2]]}]
];

SelectGoodFitsObjects[fitsObject_List]:=
Module[{pickList,frames,trackedStars,width,height,geomTransforms, scaledFrames, croppedFrames},
frames=#[1]&/@(#[Image]&/@fitsObject);
(*If I omit the conversion to Bit16 here,ImageFeatureTrack will fail to find any features.*)scaledFrames=Map[Image[#,"Bit16"]&,frames];
width=Min[(ImageDimensions/@scaledFrames)[[All,1]]];
height = Min[(ImageDimensions/@scaledFrames)[[All,2]]];
croppedFrames=ImageCrop[#,{width,height}]&/@scaledFrames;
(*For each frame,create a list of stars that can be tracked across the set*)trackedStars=ImageFeatureTrack[croppedFrames,MaxFeatures->3];
pickList=If[ContainsAny[#,{Missing[]}],False,True]&/@trackedStars;
Pick[fitsObject,pickList]
];

AlignFITSObjectList[fitsObject_List]:=
Module[{goodfitsObjects,unalignedFrames, alignedFrames},
goodfitsObjects = SelectGoodFitsObjects[fitsObject];
unalignedFrames=#[1]&/@(#[Image]&/@goodfitsObjects);
alignedFrames=AlignStarFieldSimple[unalignedFrames];
Table[
Association[{
Image->Association[{1->alignedFrames[[i]]}],MetaInformation->goodfitsObjects[[i]][MetaInformation]
}],
{i,Length[goodfitsObjects]}]
];

BL Camelopardalis

fileNames3=FileNames[FileNameJoin[{NotebookDirectory[],"Data","BL Cam","*.fit"}]];

varStarImagesBLC=LoadFITSObject/@fileNames3;

Manipulate[ImageAdjust[varStarImagesBLC[[i]][Image][1]],{i,1,Length[varStarImagesBLC],1},SaveDefinitions->True]

varStarImagesAlignedBLC=AlignFITSObjectList[varStarImagesBLC];

Manipulate[ImageAdjust[varStarImagesAlignedBLC[[i]][Image][1]],{i,1,Length[varStarImagesAlignedBLC],1},SaveDefinitions->True]

ImageAdjust[varStarImagesAlignedBLC[[1]][Image][1]]

ListLinePlot[SortBy[PhotometryObs[#,{221,337},{282,579},10,6,8,20]&/@varStarImagesAlignedBLC, First],PlotTheme->"Scientific"]

References

Lomb NR. 1976. Least-Squares Frequency Analysis of Unequally Spaced Data.  Astrophysics and Space Science 39: 447-462

Lafler J. and Kinman TD. 1965. An RR Lyrae Star Survey With The Lick 20-Inch Astrograph II.  The Calculation of RR Lyrae Periods By Electronic Computer.  Astrophysical Journal Supplement, 11: 216

Deeming TJ.  1975. Fourier Analysis With Unequally-Spaced Data. Astrophysics and Space Science, 36: 137-158

Scargle JD. 1982.  Studies in Astronomical Time Series Analysis.  II Statistical Aspects of Spectral Analysis of Unevenly Spaced Data.  The Astrophysical Journal, 263: 835-853

Clarke D. 2002.  String/Rope Length Methods Using the Lafler-Kinman Statistic.  Astronomy & Astrophysics, 386: 763-774