AutoCal : A software application for calibrating photometric data

HOW TO CITE: Wium DJ, Van Soelen B. AutoCal: A software application for calibrating photometric data. S Afr J Sci. 2016;112(3/4), Art. #2015-0034, 8 pages. http://dx.doi.org/10.17159/ sajs.2016/20150034 We present a software application for the calibration of stellar magnitudes in the absence of standard stars. It uses an existing algorithm to match stars in the target’s field of view to catalogue entries and computes the average offset between the two sets of magnitudes using a weighted least-squares approach. This offset is used to calibrate the target’s instrumental magnitude. The software application was used to calibrate magnitudes for six Be/X-ray binary sources in the Small Magellanic Cloud and the results were compared with published results for these sources. Where comparisons were possible, our results agreed with those results within the uncertainties specified. Infrared variability was found for all six of the sources tested. The interactive outlier removal that was made possible by our software allowed for smaller uncertainties to be reported for our results.


Introduction
The measured (instrumental) magnitude of a star, a measure of its brightness, is influenced by loss of light in the atmosphere. The apparent magnitude is calculated by removing this influence through observations of standard stars with well-defined apparent magnitudes. However, this is difficult if atmospheric conditions change on short timescales (shorter than the time required to undertake the observations). Alternatively, a good estimate of the apparent magnitude of a target can be calculated by comparing the instrumental magnitudes of other stars, in the field of view (FoV), to their known apparent magnitudes, and compensating for the difference.
We developed a software application that facilitates this process by matching observed stars in the FoV of the target to a catalogue and calculating the relationship between the two sets of magnitudes using a weighted least-squares approach. An existing algorithm, proposed by Groth 1 and Stetson 2 , was used for the matching. The software decreases the statistical uncertainty in the calibrated magnitude relative to what can be achieved through a manual selection of a few stars by hand, by using a large number of stars.
The software was tested on infrared observations of a number of Be/X-ray binary systems in the Small Magellanic Cloud (SMC). A good agreement was found between our results and previous results for these systems. Additionally, using our software, we obtained smaller error margins in the calibrated magnitudes, compared to previous values calculated with the same data using similar techniques.

Be/X-ray binary systems
Be/X-ray binaries consist (in general) of a neutron star which is accreting material from the Be star it is orbiting (see Reig 3 for a review, and Sturm et al. 4 and Li et al. 5 for a possible exception). The Be star is surrounded by a circumstellar disc, as is evidenced by the presence of emission lines in the spectrum of these stars. 6 In addition to these emission lines, Be stars show a larger flux in the infrared than would be expected from a blackbody (stellar) emitter. This infrared excess is also believed to originate from the circumstellar disc and the Be star's spectrum can be fitted by including a free-free component in addition to the blackbody spectrum of the underlying star. [7][8][9] The size of the circumstellar disc is also variable as is reflected in the emission lines and in the infrared excess. For example, Telting et al. 10 showed the long term variability in X Persei between a disc and discless phase. Interferometry is now allowing for more detailed studies of these systems. See Rivinius et al. 11 for a recent review.
The X-ray emission observed from Be/X-ray binaries is the result of matter accreting onto the surface of the neutron star. Two types of outbursts are observed in these systems: Type I outbursts result when the neutron star moves through the circumstellar disc of the Be star during its orbit, while Type II outbursts are much brighter and longer lasting outbursts that arise when the circumstellar disc of the Be star is larger than normal and continuous accretion occurs. Because the Type II outbursts occur during periods of increased circumstellar disc size, this corresponds to periods of increased infrared flux from the system. Observations after a Type II flare show a decrease in the infrared flux. 12,13 Long-term studies of optical and infrared variability of Be/X-ray binaries (e.g. Rajoelimanana et al. 14 ) are therefore an important tool in studying these systems.
In this paper, we used data obtained as part of long-term monitoring of Be/X-ray binaries with the IRSF telescope, to test the robustness of the software application developed. We report on our results in the sections that follow.

Data analysis
The combination of the dithers and the photometric data reduction was done with the sirius09 16 pipeline package that is integrated into the Image Reduction and Analysis Facility (IRAF) software package. 17 Point Spread Function (PSF) photometry was performed on the science frames using the standard IRAF/NOAO packages. 18 The resulting photometric output of instrumental magnitudes was converted into a list of stars recording the x and y frame coordinates of the star, and the measured instrumental magnitudes with the associated (statistical) uncertainties. This list was used with the AutoCal program and the published IRSF catalogue 19 to calibrate the target's magnitude.

The AutoCal software application
A typical calibration workflow using AutoCal can be loosely divided into the three sections of preparation: calibration, post calibration analysis and fit correction. The steps undertaken during these three phases are briefly discussed below.

The preparation phase of AutoCal
In addition to the list of magnitudes determined by IRAF, the software takes as input the FITS (Flexible Image Transport System) file that was analysed. This allows the FITS image to be displayed and offers the user visual feedback of the matching process. Information is also read from the FITS header, including the Right Ascension (RA) and Declination (Dec) at the centre of the image and its width and height in pixels. The user is prompted for any additional information not included in the header that is required, such as the telescope's plate scale. The target star is selected by clicking on the displayed image or specifying its x,y coordinates.
In order to optimise the search for catalogue magnitudes and to ensure consistent formatting, catalogue data is imported into a database that is linked into AutoCal. This import is only performed the first time a catalogue is used. Once the catalogue has been added to the database, it will remain available during future uses of the program. All database tables contain the fields RA, Dec and magnitude, followed by a corresponding error-in-magnitude field. The interface prohibits the addition or omission of fields, ensuring that the required format is adhered to.
After the instrumental magnitudes are loaded into AutoCal, a database query is initiated to locate all stars in the catalogue that fall within the FoV of the image. The FoV is determined by the RA and Dec position from the FITS header, together with the image size (px) and plate scale (arcsec/px). The user may, alternatively, enter the ranges manually or specify only a percentage of the FoV to be used.

Calibration procedure
The calibration procedure comprises the following two main steps. Firstly, stars from the frame are matched to the corresponding catalogue stars, and secondly, a linear equation describing the relationship between the magnitudes of the pairs of matched stars is calculated and applied to the target star's instrumental magnitude. These two steps are described in the sections that follow.

Star matching
The pattern-matching algorithm implemented in the program was concur rently developed by Groth 1 and Stetson 2 , with the former publishing it and the latter describing his version in lecture notes. A more detailed description of the algorithm can be found in the discussions by these authors. The algorithm starts by creating two lists of triangles, one containing all possible triangles (with stars as vertices) constructed from the stars in the frame star list and another containing all possible triangles constructed from those in the catalogue star list, and then searching for similarly shaped triangles from these two lists. The rationale behind the use of triangles is that the shape of a triangle is independent of translation, rotation, magnification and flips, and that corresponding star fields observed with different instruments and at different angles can therefore be matched easily. The algorithm then proceeds to pair vertices of matching triangles, or, in other words, stars.
In AutoCal, an initial matching is done on a user defined number of the brightest stars on the frame (default 25). From this initial matching, a coordinate transformation is calculated and applied to convert the coordinate system of the catalogue to that of the frame.
With the two sets of stellar coordinates now in the same coordinate system, entries in the two lists that are in close proximity are matched. At this stage, most of the stars in the frame list should be matched with a catalogue entry. This larger set of matches is then used to calculate and apply a new coordinate transformation between the two lists, effectively serving as a refinement to the initial transformation. Finally, all stars with positions that are in close proximity to the final transformed catalogue positions are matched to those entries. These pairs are used to perform the magnitude calibration, which is discussed in the next section.
Calculation of the calibrated magnitude of the target Following the star matching, the best linear relationship between the magnitudes of the frame stars and their matches from the catalogue has to be calculated. We implemented a weighted least-squares approach to determine this relationship, which is applied to the instrumental magnitude of the target in order to obtain its calibrated apparent magnitude. This section motivates the rationale behind using a linear relationship and explains the weighted least-squares approach in more detail.
To a first-order approximation (see Bailer-Jones et al. 20 for an explanation of why a first-order approximation is used) for the extinction correction, the catalogue and instrumental magnitudes are related by: where k  is the extinction correction coefficient (at wavelength ), X is the airmass and  is a zero-point offset between m cat and m ins . Therefore, as all observations on the same FoV occur at the same airmass (within less than 0.3%), we have calibrated the instrumental magnitudes against the catalogue magnitudes by assuming that the relation: m cat =gm ins +c Equation 2 may be used, where g and c are two constants that can be determined by a least squares fit. The value of the gradient should, ideally, be g=1 in all cases. However, we have included it as a free parameter and have used a weighted least squares (WLS) method to determine the values of g and c. Because the crowded star fields, obtained from the SMC observations generally contained a large portion of faint stars (with small signal to noise ratios) and only a small portion of bright stars (with large signal to noise ratios), this causes the fainter stars to skew the straight line fit if weights are not applied. The weight applied by AutoCal is calculated as: where Em cat and Em ins are the reported statistical errors in the magnitudes of the frame stars and catalogue stars (m cat and m ins ) respectively. The matrix form of the weighted least squares system that must be solved is: The 95% confidence interval for the apparent magnitude of the target was calculated using a two-sided Student's t value with n -2 degrees of freedom, where n denotes the number of matched star pairs. It can be shown that the apparent magnitude lies within: http://www.sajs.co.za Here,

Analysis of calibration results
Following completion of the calibration, the magnitude and uncertainty in the magnitude is determined for the target. The program displays a plot of m cat vs. m ins with the data points, the fitted line and other information relevant to the fit. The user is able to revise this calibration graph by choosing whether points should be excluded from the WLS fit. This is an interactive process and the user is able to select and/or remove stars interactively from the plot (a change in colour indicating which stars are included) by either clicking on individual stars, or by dragging a region around the stars that should be effected. This automatically updates the WLS fit, displaying the new best fit line on the data the equation of the fit, and re-calculates the target's apparent magnitude. The interactive nature of the star exclusion/inclusion allows for a rapid first calibration of the data.

Results
In total, AutoCal was tested on 15 observations of Be/X-ray binaries, each including J, H and K s bands, and successfully calculated apparent magnitudes for all 43 FITS frames for which the target could be detected (for two K s frames the target star was too faint to detect). The procedure we followed to obtain the apparent magnitudes was similar for all the data frames, so only one frame will be discussed in detail. The results of the calibrations performed by AutoCal are then summarised and finally, our results for a single observation are compared with previously published results from the same observation.   Magnitude calibration Figure 1 shows the calibration graph for the J-band observation of SXP2.37 after it was revised by the user. The solid line represents the WLS solution through the data points. The apparent magnitude calculated by AutoCal was 14.70±0.02 mag with initial values r 2 =0.85 and a slope of g=1.02 (all uncertainties in the magnitude are recorded at the 95% confidence level).
The scatter in the data is clearly larger for fainter stars. Because of this, and because some stars may still be false matches, have incorrect instrumental or catalogue magnitudes measurements, or have presented intrinsic variability, we tested the robustness of the fit by restricting the fit to magnitude values near the range of the target. We excluded data points from the calibration graph by first removing all faint stars with m ins greater than a user selected cut-off value, then removing any obvious outliers and finally removing all points more than 3σ from the straight line fit.
For the example discussed here, all stars fainter than m ins ≥20.9 mag were removed, after which the one obvious outlier (m ins =18.9 mag, m cat =19.9 mag) and two stars that were more than 3σ from the straight line fit were discarded. This removed 441 of the original 493 matches (the blue dots in Figure 1). The selection of the 20.9 mag cut-off point was based on a visual inspection and we found that for this frame, applying a cut-off of any value between 20 mag and 22 mag yielded very similar results. The apparent magnitude that was calculated using the smaller sample was 14.70±0.01 mag with r 2 =1.00 and a slope of g=1.01. For this J-band observation of SXP2.37 there was no significant change in the calibrated apparent magnitude when using the smaller sample of matches, demonstrating that the WLS method is robust despite outliers. However, for some of the other frames analysed with this method, we found larger discrepancies between the WLS fits undertaken before and after the removal of outlier data points, especially for certain frames that contained a large number of outliers. For those frames, removing the obvious outliers and faint stars with scattered magnitudes and lower signal to noise ratios, resulted in better linear fits near the magnitudes of the target stars and noticeable changes in the calibrated magnitudes. Table 1 lists the calibration results for the six sources analysed in this study. For every frame but one, matches were found for at least 70% of the available stars from the frame star list. There was one exception where only 45 stars were available for matching and only slightly more than half of the stars could be matched by the algorithm. A possible explanation for the differences became apparent when we repeated our calibration for the J-band without performing outlier removal, the results of which can be seen in the left half of Figure 2.

Results of calibration for 43 frames
The new result was a J-band magnitude of 14.24±0.04 mag. The uncertainty now matches that of the previously published results exactly and the magnitude is also closer to theirs. A plausible explanation for our uncertainties being smaller is our use of a weighted least-squares fitting that assigns smaller weights to the majority of fainter stars on the frames, coupled with the interactive outlier removal based on the appearance of the WLS plot.

Discussion
A number of the sources considered here have been previously discussed by other authors. For example, Rajoelimanana et al. 14 presented the longterm I-band variability of 38 Be/X-ray binaries in the SMC, some of which include the time period covered by the observations presented here. In addition, some of the IRSF data used in this study have been previously analysed and presented by Townsend 13 and Townsend et al. 213,22 , as was discussed earlier in the paper. Furthermore, the magnitudes calculated in this study are compared with those from the 2 Micron All Sky Survey (2MASS). 23 The times at which the observations were performed are reported in Modified Julian Date (MJD). Our results for each source are discussed individually below, and brought into the context of the existing literature.

SXP0.72
SXP0.72 was observed on 18 December 2007 (MJD 54 452.77) and 5 January 2010 (MJD 55 201.79). As can be seen from Figure 3, a comparison of the magnitudes calculated from the two observations showed no significant variability in any of the bands. Furthermore, only a very slight change was found between the 2MASS magnitudes and those calculated here (see Table 1), corresponding to a decrease in flux of between 10% and 15% for each of the three bands. There is a significant, albeit small, increase in magnitude of ~0.3 mag observed between the first and second observing runs. This slight variation is indicative of the variability associated with the circumstellar disc. This variability is clearly shown in Figure 4a of Rajoelimanana et al. 14 , that demonstrates the long-term variation.

SXP2.76
SXP2.76 is another source that showed very little infrared variation during all three observing runs, and the magnitudes are comparable to the 2MASS observations (see Table 1 and Figure 3). A very slight brightening was observed in all bands between the first and second runs (MJD 54 451.78 and 55 201.84) corresponding to an increase in flux of between 11% and 17%. The magnitudes during the third run (MJD 55 551.88) were similar to those for the second run. Rajoelimanana et al. 14 showed that SXP 2.76 shows magnitude variation, in I band, with a period of roughly 2800 days and an amplitude of ~0.2 mag (see their Figure 5a). The increase in brightness we observed in the infrared is consistent with these results.  catalogue, but is much fainter during the third run in all three bands (see Table 1 and Figure 3).
This clearly indicates the variable nature of this system, and the decrease in infrared flux during the third observation run indicates a decrease in

SXP18.3
Of the six sources covered by this study, SXP18.3 showed by far the highest level of variability in the near infrared. It was observed on MJD 55 195.85 and 3 times consecutively around MJD 55 549.80. We used the magnitude calculated from the second of these observations (listed as observation B in Table 1), because it was the only one for which the target could be identified in the K s -band. A change in magnitude of 1.02 mag and 0.9 mag was observed from the time of the 2MASS campaign to the second observing run for the J and K s -band respectively. The J and K s -band brightening is equivalent to flux increases of 156% and 129% respectively. An even greater variability was found between the second and third observation runs, equating to a decrease in flux of 68%, 69% and 76% in J, H and K s respectively. The magnitudes calculated from the data captured during the second and third observation runs are indicated in Figure 3.
From Figure 14a in Rajoelimanana et al. 14 , the optical data clearly show an extended period of brightening, starting around MJD 52 800 and lasting until around MJD 54 600, a period of about five years. The drastic decrease in the brightness of the system from December 2009 to December 2010 seems consistent with the aftermath of another Type II outburst. This idea is supported by the same figure, that shows a sharp I-band brightening for the period around MJD 54 700 to MJD 55 000, which leads up to the IRSF observation that was taken on MJD 55 549.

Conclusion
Of the 45 FITS frames that were created from the data captured during the three IRSF observations, the target was visible on 43. AutoCal was successful in calculating calibrated magnitudes for all of these frames. The correctness of these magnitudes is demonstrated by the good agreement with previously published results for the same data.
That section further argued that the interactive outlier removal through inspection of the WLS plot, together with the use of WLS methods to calculate uncertainties, allows for smaller uncertainties in the calibrated magnitudes.
AutoCal proved robust in the sense that it managed to successfully calibrate all the frames that we tested, irrespective of their stellar density. Despite the fact that all the frames can be said to be reasonably densely populated, the number of stars they contained ranged from 41 to 907. All the calibrations were performed with the user-specified parameters at their default values. The fact that AutoCal's functioning was not sensitive to these values ensured faster calibration and a good user experience. AutoCal has the potential to be expanded further to allow blind searches for variability to be undertaken.
We were able to detect statistically significant infrared variability for each of the six Be/X-ray binaries discussed here. The results for SXP11.5 agree well with those reported in Townsend 13 , while those of SXP2.37 and SXP18.3 point towards the same events that caused the variability seen in the OGLE and MACHO data presented in Rajoelimanana et al. 14