A new wavelet-based approach for the automated treatment of large sets of lunar occultation data

0.67 MB·English

by  O. Fors

Astronomy & Astrophysics manuscript no. 8987 (cid:13)c ESO 2008
February 2, 2008
A new wavelet-based approach for the automated treatment of
⋆
large sets of lunar occultation data
O. Fors1,2, A. Richichi3, X. Otazu4,5, and J. Nu´n˜ez1,2
1 Departament d’Astronomia i Meteorologia, Universitat de Barcelona, Mart´ı i Franqu´es1, E-08028 Barcelona, Spain
8 e-mail: [email protected]
0 2 Observatori Fabra, Cam´ı del’Observatori s/n, E-08035 Barcelona, Spain
0 3 European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei Mu¨nchen, Germany e-mail:
2 [email protected]
4 Computer Vision Center, Universitat Aut`onomade Barcelona, E-08193 Bellaterra, Spain
n 5 Departament de Ci`encies dela Computaci´o, Universitat Aut`onoma deBarcelona, E-08193 Bellaterra, Spain
a
J
Received November4, 2007; accepted December15, 2007

4
1 ABSTRACT

] Context.Theintroductionofinfraredarraysforlunaroccultations(LO)workandtheimprovementofpredictionsbased
h
onnewdeepIRcatalogueshaveresultedinalargeincreaseinsensitivityandinthenumberofobservableoccultations.
p Aims. Weprovidethemeans for an automated reduction of large sets of LO data. This frees theuser from thetedious
-
o task of estimating first-guessparameters for thefitof each LO lightcurve.Attheendof theprocess, ready-madeplots
r and statistics enable the user to identify sources that appear to be resolved or binary, and to initiate their detailed
t interactive analysis.
s
a Methods. The pipeline is tailored to array data, including theextraction of thelightcurves from FITS cubes. Because
[ of its robustness and efficiency, the wavelet transform has been chosen to compute the initial guess of the parameters
of the lightcurvefit.

3 Results. We illustrate and discuss our automatic reduction pipeline by analyzing a large volume of novel occultation
v data recorded at Calar Alto Observatory.The automated pipeline package is available from theauthors.
7
3 Keywords. Methods:dataanalysis–Techniques:ImageProcessing–Techniques:highangularresolution–Astrometry
5 – Occultations
0
.
1
1. Introduction Besides, the limiting sensitivity achieved in the near-IR
1
by LO at the 1.5m telescope on Calar Alto is K≈ 8mag
7
For decades, lunar occultations (LO) have occupied a spe-
0 (Richichi et al. 2006a). When extrapolated to a 4-meter
cial niche as a technique for high-angular resolution with
: classtelescopeorlarger,LOappearquitecompetitive with
v excellent performance, but relatively inefficient yield. The
even the most powerful, LBI facilities (Richichi 1997).
i diffraction fringes that are created by the lunar limb as
X Asaresult,althoughthetrendisunderstandablytode-
it occults a background source, provide a unique opportu-
r nity to achieve milliarcsecond angular resolution with sin- velop more flexible, powerful and complex interferometric
a facilities, there is some balance that makes LO still attrac-
gle telescopes also of relatively small diameter. In terms
tive at least for some applications. It should not be forgot-
of instrumentation, LO have always been simple, requiring
ten that the majorityofthe hundreds ofdirectly-measured
onlyafastphotometer.Ofcourse,theyhavethesignificant
stellar angular diameters (Richichi (2007) listed 688, and
drawback that only sources included in the apparent lu-
the numbers keep increasing)were indeed obtainedby LO,
narorbitcanbeobserved(about10%ofthesky),andthen
andthatLOarestillthemajorcontributortothediscovery
onlyatarbitraryfixedtimesandwithlimitedopportunities
of small separation binary stars.
for repeated observations. If one adds that each observa-
Tworecentdevelopments,however,haveprovidedasig-
tion only provides a one-dimensional scan of the source, it
nificant boost to the performance of the LO technique,
is clear that detailed and repeated observations are better
and have significantly enlarged its range of applications:
performed with long-baseline interferometry (LBI), when
a) the introduction of IR array detectors that can be read
available. One should, however, not forget additional im-
out at fast rates on a small subarray has made it possible
portant advantages of LO: even for complicated sources,
to provide a large gain in limiting sensitivity, and b) IR
the full, one-dimensional brightness profile can be recov-
survey catalogues that have led to an exponential increase
ered according to maximum-likelihood principles without
of the number of sources for which LO can be computed.
any assumptions on the source’s geometry (Richichi 1989).
Literally,thousands ofoccultations per nightcouldnowbe
Send offprint requests to: O. Fors potentially observed with a large telescope. We describe
⋆ Algorithm tested with observations collected at Calar Alto in this paper the details and impact of these two factors
Observatory (Spain). Calar Alto is operated by the German– for LO work. We also address the new needs imposed on
Spanish Astronomical Center (CAHA). data reduction by the potential availability of a large vol-
2 O. Fors et al.: A newwavelet-based approach for automated treatment of LO
ume oflunar occultationdataper night,by describingnew While with the previous catalogues a typical night run
approaches to an automated LO data pipeline. We illus- close to the maximum lunar phase would cover 100-150
trate both the new quality of LO data and their analysis sources over several nights , predictions with 2MASS can
bymeansofexamplesdrawnfromtheobservationoftwore- include thousands of events observable with a large tele-
centpassagesoftheMoonovercrowdedregionsinthevicin- scope over one night. Special cases, like the passage of the
ity ofthe Galactic Center,carriedoutwith array-equipped Moonovercrowded,obscuredregionsinthedirectionofthe
instruments at Calar Alto and Paranalobservatories. Galactic Center, caninclude thousands ofevents predicted
over just a few hours (Richichi et al. 2006b, Fors et al.
2006a).Fig.1illustratesthetwocases.Theincompleteness
2. Infrared arrays and new catalogues
of the catalogues without 2MASS is evident already from
A number of reasons make the near-IR domain preferable theregime5≤K≤7mag.Atevenfaintermagnitudes,but
for LO work with respect to other wavelengths. still within the limits of the technique as described here,
First, LO observations are affected by the high back- the predictions based on the 2MASS catalogue are more
ground around the Moon which, being mainly reflected numerous by several orders of magnitude.
solar light, shows an intensity maximum at visible wave- Note that the increase in the number of potential oc-
lengths. Because of the atmospheric Rayleigh scattering cultation candidates is not reflected automatically in more
(∝λ−4),thebackgroundlevelgreatlydecreasesinthenear- results. The shift to fainter magnitudes implies that the
IR. At longer wavelengths (10µm − 20µm), the thermal SNR of the recorded lightcurves is on average lower; LO
emission of Earth’s atmosphere and of the lunar surface runs based on 2MASS predictions are now likely to be less
introduces a high-backgroundlevel. efficientin detecting binarieswhen compared,for example,
Second, the spacing of diffraction fringes at the tele- to studies suchas those ofEvans etal.(1986)and Richichi
scope is proportional to λ−12. Therefore, for two LO ob- et al. (2002), especially for large brightness ratios.
servations with the same temporal sampling, one recorded
in IR will obtain a higher fringe sampling than one in the 3. Automated reduction of large sets of lunar
visible.
occultation data
Finally, at least in the field of stellar diameters, there
is an advantage to observing in the near-IR because for In general, LO data are analyzed by fitting model
a given bolometric flux redder stars will present a larger lightcurves. We take as an example the Arcetri Lunar
angular diameter. Occultation Reduction software (ALOR), a general model-
Being cheap and with a fast time response, near-IR dependent lightcurve fitting algorithm first developed by
photometers have traditionally represented the detector of one of us (Richichi 1989). Two groups of parameters
choice for LO observations. Richichi (1997) showed the are simultaneously fitted using a non-linear least squares
great increase in sensitivity possible with panoramic ar- method. First, those related to the geometry of the event:
rays,whichbyreadingonlythepixelsofinterest,permitto the occultation time (t ), the stellar intensity (F ), the in-
0 0
avoid most of the shot noise generated by the high back- tensityofthe background(B )andthe limblinearvelocity
0
ground in LO. Such arrays are now becoming a viable op- with respect to the source (V ). Second, those related to
P
tion, thanks to read-outnoises,that aredecreasingateach physical quantities of the source: for resolved sources; the
new generation of chips, and to flexible electronics allow angular diameter and; for binary (or multiple) stars, the
us to address a subarray and read it out at millisecond projected separation and the brightness ratio of the com-
rates.Richichi(1997)predictedthatan8mtelescopewould ponents.
reach between K=12 and 14mag, depending on the lunar In general, the fitting procedure is approached in two
phase and background, with an integration time of 12ms steps.First,apreliminarfitassuminganunresolvedsource
at signal–to–noiseratio (SNR)=10. Observationsonone of model is performed. To ensure convergence, ALOR needs
the8.2mVLTtelescopes,equippedwiththeISAACinstru- to be provided with reliable initial guesses. We can esti-
ment in the so-called burst mode (Richichi et al. 2006b), mate the geometrical parameters with a visual inspection
show a limiting magnitude K≈12.5at SNR=1 and 3ms in- of the data, and V is predicted. The source parameters
P
tegration time, in agreement with the decade-old predic- can be refined in a second step. This is done interactively
tion. sinceitrequiresunderstandingthenatureofeachparticular
Thesenewly-achievedsensitivitiescallforacorrespond- lightcurveandthepossiblecorrelationbetweengeometrical
ing extension in the limiting magnitudes of the catalogues and physical parameters.
used for LO predictions, and their completeness. In the As a result of that great increase in the number of po-
near-IR,untilrecentlytheonlysurvey-typecatalogueavail- tentialoccultations,wesoonrealizedthatweneededasub-
able was the Two-Micron Sky Survey (TMSS, or IRC, stantial optimization in the processes of extracting the oc-
Neugebauer & Leighton 1969) that was incomplete in dec- cultation lightcurves from the raw data and of the inter-
linationandlimitedtoK<3.Already,a1m-classtelescope active evaluation of the LO lightcurves for the estimate of
equippedwithanIRphotometerexceedsthissensitivityby theinitialparametervaluesneededforthefits.Wethende-
several magnitudes (Fors et al. 2004, Richichi et al. 1996). veloped, implemented, and tested a new automatic reduc-
The release of catalogues associated with modern all-sky tion tool,the Automatic Wavelet-basedLunar Occultation
near-infrared surveys, such as 2MASS (Cutri et al. 2003) ReductionPackage(AWLORP;(Fors2006b).Thisallowsboth
andDENIS(Epchteinetal.1997), hashelped. Ourpredic- lightcurve extraction and characterization to perform the
tion software ALOP (Richichi 1985) includes about 50 other preliminary analysis of large sets of LO events in a quick
catalogues with stellar and extragalactic sources. We have and automated fashion. In the following, we describe the
now added a subset of 2MASS with K≤11,whichincludes main parts of AWLORP, which are schematically illustrated
3.7×106 sources subject to occultations. in Fig. 2.
O. Fors et al.: A newwavelet-based approach for automated treatment of LO 3
Fig.1. Frequency of lunar occultation events as a function of K magnitude, computed on the basis of all standard
catalogues in ALOP (gray bars) and of the 2MASS catalogue only (limited to K ≤ 11, clear bars). For both cases, we
◦ ◦
haveusedthe constraintsofMoon≥25 abovehorizonandSun≤−5 belowhorizon.Left:arelativelyrich5-nightrun,
from 7 thru 11 January 2006,at Calar Alto Observatory.Right: part of the night of August 5, 2006 from Paranal,when
′
the Moon reached a minimum approach of 12 from the Galactic Center. Note the logarithmic scale.
3.1. Input data and lightcurve extraction this extraction algorithm as the default in the AWLORP de-
scription.
In the cases available to us, the LO data are stored in
Flexible Image Transport System (FITS) cubes. The num-
berofcubeframesisgivenbytheframeexposureandtotal 3.2. Lightcurve characterization
integrationtime. Additional information,such as telescope
Inaccuraciesin catalogue coordinatesand lunar limb irreg-
diameter,filter andidentificatorofthe occultedobject,are
ularities introduce an uncertainty in the predicted occul-
extracted from the FITS cube header and saved in a sepa-
tation time of about 5 to 10 seconds. To secure the effec-
rate file. In addition, the limb linear velocity and the dis-
tive registering of an occultation event, the acquisition se-
tance to the Moon as predicted by ALOP are available in a
quenceisstartedwellbeforethepredictedoccultationtime.
separate file.
This results in a very long extracted lightcurve, typically
An occultation lightcurve must be extracted from the
spanning several tens of seconds. In contrast, the fringes
recorded FITS cube file. We explored several methods for
that contain the relevant high-resolution information ex-
this purpose, among them fixed aperture integration, bor-
tend only a few tenths of a second. In addition, to accom-
der clipping, Gaussian profile and brightest-faintest pixels
plish a proper fitting of this much shorter lightcurve sub-
extraction. We found these partly unsatisfactory, among
sample, as mentioned before, we need reliable estimates of
otherthings,becauseoflackofconnectivityacrossthestel-
t , B F .
0 0 0
larimageandbecauseofsensitivitytofluxandimageshape
The problem corresponds to detecting a slope with a
variations.
known-frequency range in a noisy, equally sampled data
We addressed the problem of connectivity with the use
series. The key idea here is to note that the drop from the
of masking extraction, and two methods were considered.
first fringe intensity (close to t ) is always characterized
0
The first method, called 3D-SExtractor,consists of a cus-
by a signature of a given spatial frequency. Of course, this
tomization of the object detection package SExtractor
frequency depends on the data sampling but, once this is
(Bertin&Arnouts1996)forthecaseof3DFITSLOcubes.
fixed, the aimed algorithm should be able to detect that
The algorithm invokes SExtractor for every frame and
signatureandprovideanestimateoft ,regardlessitsSNR.
0
evaluates its output to decide if the source has been ef-
Once t is known, the other two parameters (B and F )
0 0 0
fectivelydetected.Thesegmentationmap(orsourcemask)
can be estimated.
provided by SExtractor defines the object (background)
Thisproblemcallsforatransformationofthedatathat
pixels in case of positive (negative) detection. These pixels
wouldbecapableofisolatingsignaturesinfrequencyspace,
are used to compute the source (background)intensity be-
whilesimultaneouslykeepingthetemporalinformationun-
fore and after the occultation. The second method, called
touched. Wavelet transform turns out to be convenient for
Average mask,consistsinperformingsimpleaperturepho-
this purpose.
tometry using a predefined source mask. This is obtained
by averaging a large number of frames previous to the oc-
cultation and by applying a 3σ thresholding. 3.2.1. Wavelet transform overview
We empirically compared3D-SExtractorandAverage
The wavelet transform of a distribution f(t) can be ex-
mask methods under a variety of SNR, scintillation, and
pressed as:
pixel sampling situations. Although the 3D-SExtractor
makes use of a more exactmask definition for every frame,
Average mask was found to provide less noisy lightcurves −1 +∞ t−b
W(f(t))(a,b)=|a| 2 f(t)ψ dt, (1)
with no evident fringe smoothing. Therefore, we adopted Z−∞ (cid:18) a (cid:19)
4 O. Fors et al.: A newwavelet-based approach for automated treatment of LO
Input FITS cube
Telescope
Filter
t Object id
Lightcurve extraction
6000
4000
2000
00 20000 40000 60000 80000
Long lightcurve
Lightcurve characterization
by wavelet decomposition
6000
4000 t0
2000 F0 B0
ct t0-∆t t0+∆t
e 02500 2700 2900 3100 3300 3500
bj
Trimmed lightcurve
t o
c
e xt
bj e
o N
Faint source?
xt Object excluded Wide binary?
e Y Edge cutting?
N
N
6000
Prediction parameters
Lightcurve fitting
4000 (Lunar limb rate, Moon distance)
t0F0B0
( , , )
2000
Data Model Residuals
8000
400
0
−400 Plot generation
−8002500 2700 2900 3100 3300
Fig.2. Flow-chartdescription of AWLORP.
where a and b are scaling and translational parameters The differences between two consecutive approxima-
respectively. Each base (or scaling) function ψ(t−ab) is a tions fi−1(t) and fi(t) are the wavelet (or detail) planes,
scaledandtranslatedversionofa function ψ calledmother w (t). Letting f (t)=f(t), we can reconstruct the original
i 0
wavelet, satisfying the relation ψ(t−b)=0. signal from the expression:
a
We followed the a` trouRs algorithm (Starck &
n
Murtagh 1994) to obtain the discrete wavelet decomposi-
f(t)= w (t)+f (t) , (3)
tion of f(t) into a sequence of approximations: i r
Xi=1
F (f(t))=f (t), F (f (t))=f (t),··· . (2)
1 1 2 1 2
where f (t) is a residual signal that contains the global
r
f (t) (i = 1,···,n) are computed by performing suc- energyoff(t).Notethatn=r,butweexplicitlysubstitute
i
cessive convolutions with a filter derived from the scaling n with r to clearly express the concept of residual. Each
function,whichinthis caseis aB cubic spline.The useof wavelet plane can be understood as a localized frequential
3
a B cubic spline leads to a convolution with a mask of 5 representationatagivenscaleaccordingtothewaveletbase
3
elements, scaled as (1,4,6,4,1). function used in the decomposition.
O. Fors et al.: A newwavelet-based approach for automated treatment of LO 5
In our case, we are using a multiresolution decomposi- 5. F is computed by subtracting B to I .
0 0 p
tionscheme,whichmeansthe originalsignalf(t)hastwice
the resolution of f (t). This latter has twice the resolution
1
As in the case of the 7th plane, the contribution in the
of f (t), and so on.
2
5th plane is dominated by signal features represented at
thisscale,whilenoise,eventhescintillationcomponent,has
3.2.2. Algorithm description a minor presence. Therefore, again, the estimation criteria
for B and F is likely to be well behaved and robust in
0 0
We developed a program to perform a discrete decomposi-
presence of high noise.
tion of the lightcurve into n wavelet planes. Note that
wav Although AWLORP was demonstrated on a particular
the choice of n depends exclusively on the data sam-
wav data set, its applicability is totally extensible to any sam-
pling and will be discussed later. For example, n = 7
wav pling of the lightcurve and also to reappearances. To show
was empirically found to be a suitable value for represent-
this, we repeated the previous algorithm description for 6
ing all the features in the frequency space of the lightcurve
sets of 100 simulated2 lightcurves of different samplings
when the sampling was 8.4 ms. The 2nd to 7th wavelet
(1,2,4,6,8 and 10ms). For these six samplings, n was
planes resulting from the decomposition of the lightcurve wav
found to be 8,7,6,6,5 and 5, respectively. Note these values
of the bright star SAO 190556 (SNR=43) are represented
are proportionalto a geometric sequence of ratio 2 and ar-
in Fig. 3. The 1st plane was excluded as it nearly exclu-
gument (8−n ), which is in agreement with the dyadic
sively contains noise features not relevant for this discus- wav
nature of the wavelet transform we adopted.
sion.Forthesakeofsimplicity,wewillconsiderthispartic-
ular lightcurve and sampling value in the description that
follows. 3.3. Lightcurve fitting
Wedesignedanalgorithmwhichestimatest ,B andF
0 0 0
from the previous wavelet planes. This consists of the fol- The algorithmjustdescribedhas beenintegratedinanau-
lowingtwosteps:First,itwasempiricallydetermined1 that tomated pipeline. As shown in the scheme of Fig. 2, the
the 7th plane serves as an invariant indicator of the occul- characterization of the lightcurve is used to decide if a fit
tation time (t ). In particular, t coincides approximately can be performed succesfully with ALOR. The cases of very
0 0
withthezerolocatedbetweentheabsoluteminimum(tmin) faintsources,widebinariesandthoselightcurveswithsome
0
and maximum (tmax) of that plane (see upper right panel data truncation (i.e. very short time span on either side of
0
in Fig. 3 for a zoomed display of the 7th plane). The good the diffraction fringes) are the typical exclusions, and are
localization of t in this plane is justified because the first discussedinSect.4.3.Incaseofpositiveevaluation,ALORis
0
fringemagnitudedropismostlyrepresentedatthiswavelet executedusingthedetectedvaluesoft ,F ,andB asinitial
0 0 0
scale. In addition, the presence of noise is greatly dimin- guesses.After the preliminaryfit is performed,a quicklook
ishedinthisplane.Thisisbecausenoisesources(electronics plot of lightcurve data, model, and residual files is gener-
orscintillation)contributeathigherfrequencies,andthere- ated. This process is iterated for all the observed sources.
forearebetterrepresentedatlowerwaveletscales(planes). This automatic pipeline frees us from the most tedious
Inother words,this criteriafor estimating t0 is likely to be anderror-pronepartofALORreduction.Thepipelinespends
insensitive to noise, even for the lowest SNR cases. a few seconds per occultation to complete the whole pro-
Second,onceafirstestimateoft0 wasobtained,B0 and cess described in Fig. 2. For comparison, an experienced
F0 could be derived by considering the 5th wavelet plane. user takes 10-20 minutes per event for reaching the same
We found that this plane indicates those values with fairly stageofthe reductionpipeline. Incaseswhenthe data sets
good approximation.The procedure is illustrated in Fig. 3 included hundreds of occultation events, this difference is
and is described as follows: substantial. The pipeline was coded entirely in Perl pro-
gramming language, which turns our to be a powerful and
1. Weconsidertheabscissainthe5thplane,corresponding flexible tool for concatenating the I/O streams of indepen-
to t0 found in the 7th plane. dent programs.
2. From t0, we search for the nearby zeroes in the 5th Once AWLORPhas automaticallygeneratedallthe single
plane,beforeandaftertheaboveabscissa.Wecallthem sourcefitplots,the usercanperforma quickvisualinspec-
tb and ta. tion.Theobjectiveofthisfirstevaluationistoseparatethe
3. We estimate B0 by averaging the lightcurve values unresolved, relatively uninteresting events from those that
aroundta within aspecifiedtime range.We empirically bear the typical marks of a resolved angular diameter, of
fixed this to [−8,8] samples because it provided a good an extended component or of a multiple source. These lat-
compromise between improving noise attenuation and ter will still need an interactive data reduction with ALOR,
suffering from occasional backgroundslopes. buttheywillrepresenttypicallyonlyasmallfractionofthe
4. The same window average is computed around tb. The whole data set.
obtainedvalue(I )representsameanvalueoftheinten-
p
sity at the plateau regionbefore the onset of diffraction
fringes. Note that the 5th waveletplane waschosenbe-
cause its zero at t is safely before the fringes region in 4. Performance evaluation
b
the lightcurve, where the intensity is not constant and,
We have verified the performance of AWLORP by analysing
thus, not appropriate for I calculation.
p
both simulated and real LO data sets.
1 This was realized by repeating the same analysis to many
other lightcurves of different SNR values and same time sam- 2 The procedure folowed to simulate these data sets is ex-
pling (8.4 ms). plained in Sect. 4.1.
6 O. Fors et al.: A newwavelet-based approach for automated treatment of LO
750 750 750 750
tmax
0
500 500 500 500
250 250 250 250
0
0 0 0
−250 t0
−250 −250 −250
−500
−500 −500 −500 tmin
i=7 0
i=2 i=3 i=4 −75104000 15000 16000 17000 18000 19000
−750 −750 −750
0 10000 20000 30000 0 10000 20000 30000 0 10000 20000 30000
750 750 750 5000
tb
4000
500 500 500 3000
250 250 250 2000 F0 ta
1000
B0
0 0 0 7500
500
−250 −250 −250 250
0
−500 −500 −500 −250 tb ta
−500
−750 i=5 −750 i=6 −750 i=7 −75105500i=5 16000 16500 17000 17500
0 10000 20000 30000 0 10000 20000 30000 0 10000 20000 30000
Fig.3.Schematicofthewavelet-basedalgorithmfortheestimationoft ,F ,andB tobeusedinAWLORP.Thelightcurve
0 0 0
corresponds to an occultation of SAO 190556 observed at the Calar Alto Observatory, sampled every 8.4 ms. Left: box
with 2nd to 7th waveletplanes resulting from the waveletdecomposition of the originallightcurve. Upper right: the 7th
planeisfoundtobe agoodindicatoroft .Azoomeddisplayoftheregionaroundt isshown.Lowerright:aboxdisplay
0 0
of 5th plane (bottom part of this panel) provides the abscissae t ,t to compute F and B in the original lightcurve
b a 0 0
(upper part of the same panel).
4.1. Simulated data regime. In addition, the typical width of the ∆t distribu-
0
tion is inversely proportional to the SNR value. A gaus-
Thanks to a specific module included in ALOR, a set of
sian function was fitted to every histogram, and we found
simulated LO lightcurves was generated for varying SNR
the values σ =23.0,11.7,4.6,2.3,1.1,0.5for the cases with
values. The noise model assumes three independent noises
SNR=1,2,5,10,20,50.
sources:detectorelectronics,photonshot-noise,andscintil-
lation, which are of Gaussian, Poisson, and multiplicative Second, note that the histograms in Fig. 4 are not ex-
nature, respectively (Richichi 1989). With a realistic com- actly centered at ∆t0 = 0, but systematically shifted 4ms
binationofthesethreenoisesources,wegeneratedsixseries to 2ms (only 2 to 1 sampling points). This error is about
with SNR 50,20,10,5,2 and 1, each of them consisting of the Nyquist cut-off frequency of our data sampling. It can
10000lightcurves.We chosethe samplingto be2ms,which be assumed as a limitation imposed by the data and not
is a realistic value considering what is offered by current as an intrinsic constraint of AWLORP. The difference could
detectors. becorrectedbysubtractingthissmalloffsettoallanalyzed
AWLORPwasexecutedforallthe60000simulatedevents. lightcurves, but it is in any case of no consequence for the
purpose of the subsequent interactive analysis.
For each lightcurve, we found an estimate of the triplet
(t ,F ,B ). The AWLORP only failed to characterize the
0 0 0
lightcurve in 10 cases of the noisiest series for which the
4.2. Real data
ALORfitscouldnotconverge.Fortheremaining59990cases,
wecomputedthedifference(∆t )betweenthedetectedand
0
We considered a set of six real lightcurves. These were
thesimulatedoccultationtimeandplottedthesedifferences
recorded in the course of Calar Alto Lunar Occultation
as shown in Fig. 4. Two comments can be made.
Program(CALOP)(Richichietal.2006a,Forsetal.2004).
First, the ∆t distribution is, to a good approxima-
0
They correspond to a series of SNR values similar to the
tion, Gaussian-shaped. This is in agreement with the fact
one discussed in Sect. 4.1.
that the first fringe signature is primarily dominated by
Gaussian noise at the wavelet plane (n = 7) employed The robustness of t estimation is shown in Fig. 5,
wav 0
to estimate t . This noise distribution has its origins in where even in the lightcurves at the limit of detection
0
the detector read-out for the faint end (low SNR) and (SNR=1.2,2.1)thevalueoft iscorrectlydetected.Thisis
0
in the shot-noise for the bright end (high SNR), which confirmed by visual inspection and by an comparison with
can be approximated by a Gaussian distribution in this the predicted values.
O. Fors et al.: A newwavelet-based approach for automated treatment of LO 7
500
3000
SNR=50 400 SNR=5
2000 300
200
1000
100
0 0
-100 -50 0 50 100 -100 -50 0 50 100
500
3000
SNR=20 400 SNR=2
2000 300
200
1000
100
0 0
-100 -50 0 50 100 -100 -50 0 50 100
500
3000
400
SNR=10 SNR=1
2000 300
200
1000
100
0 0
-100 -50 0 50 100 -100 -50 0 50 100
∆t (ms) ∆t (ms)
0 0
Fig.4. Application of AWLORP to six sets of 10000 simulated lightcurves at 2ms sampling and of different SNRs values.
As explained in the text, the offset between the simulated occultation time and the time detected by AWLORP (∆t ) is
0
Gaussian distributed with FWHMs inversely proportionalto the SNR value, and the histogrampeaks are sistematically
shifted within the range ∆t ∼[−4,−2]ms (only 1 to 2 sampling points).
0
To verify this concordance, we ran ALOR fits for all six of the lightcurves, since this is the size of the scaling
lightcurves with the AWLORP-detectedtriplets (t ,F ,B ) as functionatthe scaleofthe 7thplanefor the giventem-
0 0 0
initial values. Even in the faintest cases, ALOR converged poral sampling.
for all parameters of the lightcurve model. With regard as 3. Depending on the subarray size employed, the image
t , the difference between the initial and the fitted values scale, the seeing conditions or telescope tracking, part
0
never exceeded 13.6 ms (1.6 sample points) as can be seen ofthe stellar imagemightbe displacedoutside the sub-
in Fig. 5. arraysothattheextractedfluxdecreasesandtheshape
of lightcurve is affected. Under these circumstances,
AWLORP is likely to produce false t detections. Again,
4.3. Problematic cases 0
the small number of cases affected does not justify the
Thepipelinejustdescribedworkswellforabout98%ofthe substantialeffortrequiredtoimprovetheAWLORPtreat-
recordedevents.Thereare,however,afewspecialsituations ment.
where the algorithm of Fig. 2 fails. Those can be classified
in three distinctive groups: 5. Summary
1. The current version of wavelet-based lightcurve char- Theobservationoflunaroccultation(LO)eventswithmod-
acterization does not support wide binary events. In ern infrared array detectors at large telescopes, combined
other words, the pipeline cannot simultaneously deter- with the use of infrared survey catalogues for the predic-
minethevalues(tA,BA andFA)and(tB,BB andFB) tions, has shown that even a few hours of observation can
0 0 0 0 0 0
for two components A andB separatedby morethan a resultinmanytens if nothundredsofrecordedoccultation
hundred of milliarcseconds. Since these cases represent lightcurves. The work to bring these data sets to a stage
atmostafewpercentoftheoverallvolumeofLOevents whereanexperiencedobservercanconcentrateonaccurate
and they are also relatively uninteresting, this feature interactive data analysis for the most interesting events is
has not been implemented yet. long and tedious.
2. Due to observational constraints, to an unusually large We have designed, implemented, and tested an au-
prediction error or simply by mistake, sometimes the tomated data pipeline that takes care of extracting the
recording of an event is started too close to the actual lightcurves from the original array data (FITS cubes in
occultation time. Since the scaling function has a given our case); of restricting the range from the original tens
size at each wavelet scale, there is a filter ramp that of seconds to the few seconds of interest near the occul-
extends over an initial span of data depending on the tation event; of estimating the initial guesses for a model-
waveletplane.Forexample,inthecaseofdatainSect.5 dependentfit;ofperformingthefit;andfinallyofproducing
thishappensupto4000millisecondsfromthebeginning compact plots for easy visual inspection. This effectively
8 O. Fors et al.: A newwavelet-based approach for automated treatment of LO
3500 3500
3000 3000 tf
0=28501.2
2500 2500
2000 2000 td0=28505.6
1500 1500
1000 1000
500 500
0 0
0 10000 20000 30000 40000 50000 60000 27500 28000 28500 29000 29500
1500 1500
1000 1000 ttd0f0==3333118858..74
500 500
0 0
0 10000 20000 30000 40000 50000 60000 32000 32500 33000 33500 34000
700 700
560000 560000 td0=27430.9
tf
400 400 0=27441.0
300 300
200 200
100 100
0 0
0 10000 20000 30000 40000 50000 60000 26500 27000 27500 28000 28500
700 700
600 600 ttd0f=27366.1
500 500 0=27367.2
400 400
300 300
200 200
100 100
0 0
0 10000 20000 30000 40000 50000 60000 26500 27000 27500 28000 28500
500 500
400 400 td0t=f0=303502584.29.5
300 300
200 200
100 100
0 0
0 10000 20000 30000 40000 50000 60000 29500 30000 30500 31000 31500
400 400
ttd0f=27240.8
300 300 0=27241.2
200 200
100 100
0 0
0 10000 20000 30000 40000 50000 60000 26000 26500 27000 27500 28000
Fig.5.ApplicationofAWLORPto6lightcurveswithdifferentSNRs(fromtoptobottom:47.2,22.3,10.9,5.9,2.1and1.2)
observed as part of the CALOP program. The left side panels show the whole lightcurves (60 seconds). The right side
panelsshowthe trimmedlightcurves(spanningonly2 seconds)aroundthetd value detectedby AWLORP.The occultation
0
timefittedbyALORusingtd asinitialvalue,tf,arealsodisplayed.NotethateveninthefaintestSNRcase,theoccultation
0 0
time is correctly detected.
reduces the time needed for the initial preprocessing from lightcurvecharacterizationwas correctandconsistentwith
severaldaystoafewhours,andfreestheuserfromarather the simulated values. Convergence could not be reached
tedious and error-prone task. The pipeline is based on an due to poor signal-to-noise ratio in only in ten cases out
algorithm for automated extraction of the lightcurves, and of 60000. These cases would, of course, be challenging for
on a wavelet-basedalgorithm for the estimation of the ini- an interactive data analysis by an experienced observer as
tial parameter guesses. well. We also tested the pipeline ona set of realdata,with
The pipeline has been tested on a large number of sim- similar conclusions. We identified and discussed the cases
ulatedlightcurvesspanningawiderangeofrealisticsignal- that may prove problematic for our scheme of automated
to-noiseratios.Theresulthasbeencompletelysatisfactory: preprocessing.
in all cases in which the algorithm converged, the derived
O. Fors et al.: A newwavelet-based approach for automated treatment of LO 9
Acknowledgements. This work is partially supported by the ESO
Director General’s Discretionary Fund and by the MCYT-SEPCYT
Plan Nacional I+D+I AYA#2005-082604.
References
Bertin,E.,&Arnouts,S.1996,A&AS,117,393
Cutri, R. M., Skrutskie, M. F., van Dyk, S. et al. 2003, The IRSA
2MASS All-Sky Point Source Catalog, NASA/IPAC Infrared
ScienceArchive
Epchtein,N.,etal.1997, TheMessenger,87,27
Evans, D. S., McWilliam, A., Sandmann, W. H., & Frueh, M. 1986,
AJ,92,1210
Fors,O.,Richichi,A.,Nu´n˜ez,J.,Prades,A.2004, A&A,419,285
Fors O., Richichi, A., Mason, E., Stegmeier, J., Chandrasekhar, T.
2006, Highlights of Spanish Astrophysics IV, Figueras et al.
(Eds.),Springer2007
Fors, O. 2006, Ph.D. Thesis, Departament d’Astronomia i
Meteorologia, UniversitatdeBarcelona
Heckmann, O. 1975, Hamburg-Bergedorf: Hamburger Sternwarte,
1975, editedbyDieckvoss,W.
Neugebauer, G., & Leighton, R. B. 1969, NASA SP, Washington:
NASA,1969
Richichi,A.1985,Thesis,FacultyofPhysics,FlorenceUniversity(in
italian)
Richichi,A.1989,A&A,226,366
Richichi,A.,Baffa,C.,Calamai,G.,Lisi,F.1996,AJ,112,2786
Richichi, A. 1997, IAU Symp. 158: Very High Angular Resolution
Imaging,158,71
Richichi,A.,Calamai,G.,&Stecklum,B.2002,A&A,382,178
Richichi,A.,&Percheron,I.2002,A&A,386,492
Richichi,A.,Fors,O.,Merino,M.etal.2006,A&A,445,1081
Richichi,A.,Fors,O.,Mason,E.,Stegmeier,J.2006,TheMessenger,
126,24
Richichi, A. 2007, ESO Workshop The power of optical/IR interfer-
ometry:recentscientificresultsand2ndgenerationVLTIinstru-
mentationi, Richichi, A., Delplancke, F., Paresce, F., Chelli, A.
(eds.),Springer2007,p.31.
Starck,J.-L.,&Murtagh,F.1994,A&A,288,342
List of Objects
‘SAO 190556’on page 5
‘SAO 190556’on page 6