Spidast: a new modular software to process spectrointerferometric measurements^{†}^{†}thanks: Applied on observations made with ESO telescopes at the Paranal Observatory under Belgian VISA Guaranteed Time programme ID 083.D029(A/B), 084.D0131(A/B), 086.D0067(A/B/C)
Abstract
Extracting stellar fundamental parameters from SPectroInterferometric (SPI) data requires reliable estimates of observables and with robust uncertainties (visibility, triple product, phase closure). A number of fine calibration procedures is necessary throughout the reduction process. Testing departures from centrosymmetry of brightness distributions is a useful complement. Developing a set of automatic routines, called SPIDAST (made available to the community) to reduce, calibrate and interpret raw data sets of instantaneous spectrointerferograms at the spectral channel level, we complement (and in some respects improve) the ones contained in the amdlib Data Reduction Software. Our new software SPIDAST is designed to work in an automatic mode, free from subjective choices, while being versatile enough to suit various processing strategies. SPIDAST performs the following automated operations: weighting of nonaberrant SPI data (visibility, triple product), fine spectral calibration (subpixel level), accurate and robust determinations of stellar diameters for calibrator sources (and their uncertainties as well), correction for the degradations of the interferometer response in visibility and triple product, calculation of the CentroSymmetry Parameter (CSP) from the calibrated triple product, fit of parametric chromatic models on SPI observables, to extract model parameters. SPIDAST is currently applied to the scientific study of 18 cool giant and supergiant stars, observed with the VLTI/AMBER facility at medium resolution in the K band. Because part of their calibrators have no diameter in the current catalogs, SPIDAST provides new determinations of the angular diameters of all calibrators. Comparison of SPIDAST final calibrated observables with amdlib determinations shows good agreement, under good and poor seeing conditions.
keywords:
methods: data analysis – techniques: interferometric – stars: latetype1 Introduction
The power of opticalinfrared interferometry to obtain information about the astronomical source morphology (including the angular size) is now wellestablished. To properly determine the source properties, the quality of the measurements is an issue, which is still the subject of active research. As with any other measuring apparatus, the absolute calibration of the instrument (including atmosphere) requires careful attention. Derived from the measurement of the mutual degree of coherence of the incident radiation field, on spatial frequencies sampled by the aperturearray configuration, the final interferometric observables are nonlinear mixes of noisy quantities, and of parameter estimates with their own uncertainties.
In this paper, we propose to revisit and extend the existing data processing and calibration methods, in the aim to obtain reliable estimates and robust uncertainties for calibrated measurements of visibility and complex triple product. The careful reduction process, described in the present paper, has been elaborated for the scientific study of a sample of 18 bright cool giant and supergiant stars (see Table 1). Measurements were obtained with the VLTI/AMBER facility at medium spectral resolution (=1500) in the band, using triplets of 1.80m auxiliary telescopes. Observations have been conducted during 15 observing nights between May 2009 and December 2010 (under Belgian VISA Guaranteed Time). The aim of the present work is to test the data reduction and calibration software, on data with various qualities, that we started to develop in 2006.
Name  Spec. Type  Calibrator(s)  

Car  F0II  10.6(6)  1.3(3)  Col/ Eri/HR 3282 
Cet  K0III  33.9(2)  0.3(4)  Cet 
TrA  K2II  8.4(2)  1.2(1)  TrA/ Lib/ Sgr 
Hya  K3IIIII  18.1(2)  1.1(2)  Hya 
Ara  K3III  6.7(2)  0.6(2)  TrA/ Sgr 
CMa  K2.5Iab  0.2(4)  0.3(3)  HR 2411 
Oph  M0.5III  19.1(2)  1.2(2)  Lib/ TrA 
Hyi  M2III  15.2(1)  1.0(4)  Ret 
Ori  M3III  5.0(23)  0.7(2)  HR 2411 
Lib  M3.5III  11.3(3)  1.4(2)  51 Hya 
Ret  M4III  7.0(1)  0.5(3)  Ret 
L Pup  M5IIIe  15.6(10)  1.8(1)  HR 3282/ Col 
CE Tau  M2Iabb  1.8(3)  0.9(2)  40 Ori 
T Cet  M5.5Ib/II  3.7(5)  0.8(3)  Scl/ Eri 
TX Psc  C7,2(N0)(Tc)  3.6(4)  0.5(3)  Psc 
R Scl  C6,5ea(Np)  2.1(15)  0.1(1)  Eri 
W Ori  C5,4(N5)  2.6(10)  0.5(4)  HR 2113/40 Ori 
TW Oph  C5,5(Nb)  3.7(12)  0.5(4)  Sgr/ Lib 
2 Defining the basic interferometric observables
The coherent flux is provided by the measurement of the instantaneous observables delivered by the AMBER Data Reduction Software amdlib^{1}^{1}1www.mariotti.fr/data_processing_amber.htm. This quantity traces the sinelike molulated component superimposed on the continuum component, in the observed intensity distribution. It is computed using a linear fit of the individual interferograms in the detector plane, for each spectral channel (see Tatulli et al. 2007, for details).
At this point, we define the other interferometric observables:

the squared flux, which is the product of the photometric fluxes and , associated with the baseline
(1) 
the squared visibility, which is the ratio of the squared modulus of the coherent flux to the squared flux
(2) 
the complex bispectrum (Weigelt 1977), which is the product of the 3 coherent fluxes contributing to the baseline triplet
(3) 
the triple product, which is the ratio of the bispectrum to the product of the 3 fluxes contributing to the baseline triplet
(4) 
the closure phase, which is the argument of the complex bispectrum (or triple product)
(5)
3 Reducing the raw data
To process the temporal series of frames (composing the “exposures”) produced by the AMBER instrument, we use the data reduction software amdlib, which consists of a core library of C functions plus a highlevel interface in the form of a Yorick^{2}^{2}2sourceforge.net/projects/yorick/ plugin (a highlevel language environment, similar to IDL or Matlab, which is in the public domain).
The standard procedure to extract the basic interferometric observables from the AMBER raw data includes a frameselection step, which provides a set of “good” frames according to a given criterion, or a sequence of several selection criteria (Malbet et al. 2011). The criterion used the most is based on the Signal to Noise Ratio (SNR) of the fringecontrast in each frame. Two selection methods are provided in standard by amdlib:

retaining a given percentage of frames sorted according to their SNR values (called the “percentage” method),

retaining the frames with SNR values higher than a given threshold (called the “threshold” method).
With both methods, the choice of the optimal value (assumed not biasing the final results) of the selection criterion (in percentage of threshold) is done a posteriori, from a sequence of processings using different values of the selection criterion (see Millour et al. 2007, for details). Two drawbacks of this are: the need for a severalsteps procedure (to be performed manually), and the risk for biasing the final results if the selection criterion is illdetermined.
In our method, that we have implemented in a specific Yorick script added to amdlib, we apply the automatic procedure:

remove only the aberrant measurements (in squared visibility and triple product), for each spectral channel independently. Using this specific script results in amounts of useful data larger than the amount obtained with the standard procedure, without increasing biasing effect (only the lower and upper tails of the input data histograms are rejected);

assign a weight to each individual measurement. After many preliminary tests, we decided to choose weights based on the SNR of the coherent flux, and on the piston deviation;

compute the temporal weighted average of the squared visibility and triple product for each spectral channel, with uncertainties derived from the bootstrap technique invented by Efron (1979).
These automatic operations mentioned above, needing neither severalstep selection nor manual choice of the optimal value of the selection criterion, are described in more details in the following subsections.
3.1 Rejecting the aberrant data
To remove the aberrant measurements, we use a combination of physical and statistical basic criteria:

strictly positive sum and product of the fluxes in the photometric channels, since negative flux measurements result from shortcomings in the determination of the continuum in each spectrum;

coherentflux SNR higher than unity, since inclusion of data with low SNR values reduces the reliability of the estimators of the interferometric observables (Millour et al. 2008);

squared visibility greater than and lower than (“boxplot filtering”), where denotes the first quartile, the third quartile, and the interquartile range (Hoaglin et al. 1983). With triple product data, this third criterion applies to the squared visibility for each baseline, and boxplot filtering on is added to the list of criteria.
Based on the removal of aberrant measurements at the spectral channel level, rather than on the selection of the “best” frames, our approach provides an amount of useful data larger than the amount obtained with the standard frameselection procedure. For data shown in Fig. 1, obtained for a “good” seeing on the calibrator target Lib, the number of removed input data varies only slightly from one spectral channel to the other (typically less than 3%). Besides, various amount of data are rejected, according to the seeing conditions.
Using the seeing parameters (seeing angle) and (coherence time), the percentages of rejected data, respective to three baselines and two seeing conditions are given in the following lines (where “good seeing” refers to =0.5 , =13 , and “poor seeing” refers to =1.7 , =1.3 ):

shortest baseline A0–D0 (32 ): 2% (good seeing), to 9% (poor seeing);

intermediate baseline D0–H0 (64 ): 2% to 16%;

longest baseline A0–H0 (96 ): 2% to 25%;

in addition, for the baseline triplet A0–D0–H0, the rejection of data amounts 7% and 50%.
Figure 2 shows the variation of the final squared visibility (at 2.2 μ) produced with amdlib, w.r.t. the percentage of selected frames, based on the SNR and pertaining to the calibrator target Lib, observed with the shortest baseline A0–D0 (32 ) in good and poor seeing conditions (top panel and bottom panel respectively). The squared visibility produced with our method is shown for comparison.
In this example, the position of the plateau of each curve indicates a value of the selection threshold of 50%, giving the a posteriori determination with amdlib. There we find that our automatic method gives a squaredvisibility value close to the 50% value, without a need for the severalsteps procedure, in good and poor seeing conditions. Besides, Fig. 3 shows that the standard amdlib frameselection method keeps the extremal V values in poor seeing condition, which may bias the result of the calculation of the average over the frames contained in each exposure. On the contrary, the boxplot filtering used with our method rejects only the lower and upper tails of the V histogram, thus leading to estimates of the temporal average, more reliable than with the frameselection method.
Next subsection describes the weighting scheme used to average the nonaberrant data over each exposure. Note that selecting the frames with the standard amdlib procedure is equivalent to assigning unity weights to the selected frames, and null weights to the others.
3.2 Computing the raw observables
To reduce the instrumental effects on the data, at sthe pectral channel level, we compute in each exposure (300 frames in our example) the weighted average of the squared visibility, and of the complex triple product. For the visibility, as weight associated with each frame we use the ratio of the SNR to the relative excursion of the piston, where the relative excursion of the piston itself is the ratio of the piston excursion (absolute value) for a given frame, to the average over the whole exposure. For the triple product, each individual weight is given by the geometrical mean of the three weights associated with the three baselines. The final uncertainties are derived from the variance of the distribution of the weighted means, obtained by random sampling with replacement of the original series of data (directbootstrap method) (Efron & Tibshirani 1993). Note that our specific script computes the real and imaginary parts of the triple product, from which we extract the calibrated closure phase (see Sect. 4.2).
Figure 4 shows an example of squared visibility and triple product produced by the reduction script (before calibration) from one single exposure of 300 input frames, obtained with the calibrator target Lib, for the three baselines A0–D0 (32 ), D0–H0 (64 ), and A0–H0 (96 ).
In complement to these interferometric observables, we compute two other quantities used further in the calibration procedure (see Sect. 4.2.3):

the weight of each exposure, defined as the product of two ratios: one is the average of the SNR to its variance; the other is the average of the inverse of the piston excursion (absolute value) to its variance;

the final number of nonaberrant data, for each spectral channel.
4 Calibrating the data
The calibration process is the key point to obtain accurate estimates of the “true” (intrinsic) observables of the science targets. To correct the measurements for environmental and instrumental instabilities, we observe reference stars (the calibrators), with known angular diameters and independently welldefined brightness distributions.
To calibrate the wavelengthdependent measurements obtained with VLTI/AMBER, we have developed a library of IDL functions, included in the SPIDAST (SPectroInterferometric Data Analysis Software Tool) modular software suite (Cruzalèbes et al. 2008, 2010), which allows one to:

link each spectral channel to a wavelength value. We use a specific method based on the correlation of calibrator measured spectra with synthetic templates given by the MARCS model (Gustafsson et al. 2008);

measure and correct for the degradation of the spatial coherence. As calibrators, we use stars with brightness distribution of LimbDarkened (LD) discs, for which angular diameters are given by fits of synthetic spectra on the wideband spectrophotometric measurements in the infrared (Cruzalèbes et al. 2010).
4.1 Computing the spectral shifts
To reach a high precision in the angular diameter estimation using model fitting, spectral calibration is a critical point (see e.g., Wittkowski et al. 2008; Domiciano de Souza et al. 2008; Stefl et al. 2011). Given the lack of any internal instrumental module for wavelength calibration in the optical setup of AMBER (RobbeDubois et al. 2007; Petrov et al. 2007), amdlib provides calibrated wavelength tables, computed from a theoretical polynomial dispersion law (Mérand et al. 2010), but with coefficients still badly known, leading to wavelength shifts up to 10 in the detector plane, at medium spectral resolution (Malbet et al. 2011). To correct for this drawback, amdlib shifts the wavelength table using the correlation of the measured spectrum with a template table containing the telluric lines.
To improve the precision of the wavelength calibration, observers usually compute the coefficients of the dispersion law by identifying visually some prominent telluric lines in the measured spectrum (e.g., Kraus et al. 2009; Weigelt et al. 2011; Ohnaka et al. 2011). We propose an alternative automatic approach, based on the crosscorrelation of the calibrator measured spectra with synthetic templates produced by the MARCS + Turbospectrum codes (Alvarez & Plez 1998; Plez 2012), including the telluric lines as well. This method was previously suggested by Plez (2003) for Gaia, and by Decin et al. (2004) or Decin & Eriksson (2007) to calibrate the SPITZER spectrograms. It is particularly suitable for stellar spectra showing easily identifiable spectral features, as it is the case with our sample of cool giant calibrators, showing strong CO bands in the 2.126–2.474 μ spectral range (MartíVidal et al. 2011). Figure 5 shows the synthetic spectrum produced by the MARCS + Turbospectrum code, for the calibrator Lib. Some reference spectral lines are shown for identification of the corresponding stellar atmospheric elements.
For each exposure, obtained with each calibrator, and with the working instrumental setting, our automatic spectral calibration process contains the following steps:

apply the heliocentric and systemic radial velocity corrections, and multiply the synthetic spectrum by the atmospheric and instrumental transmittance profiles. Atmospheric transmission data for the southern sites (CTIO, Chile) are given by the USAF atmospheric code PLEXUS (Cohen, Wheaton & Megeath 2003);

remove the continuum parts of the raw and synthetic spectra, and normalize the resulting spectra. As a first approximation, we estimate the continuum by decreasing the spectral resolution to =40. The main drawback of this method is to produce apparent “pseudocontinua” that are lower than the real continuum levels, even in almost linefree regions (Rix et al. 2004). Since our goal is only to perform spectral calibration by correlation, this drawback has no effect on the final result;

divide the spectrum in contiguous subwindows of the same size, and compute the wavelength shift by phasecorrelation (Vera & Torres 2008) in each subwindow. If the spectrum is divided in subwindows, we apply a polynomial law of degree to correct for the wavelength shifts. Subpixel precision is obtained by embedding the crosspower spectrum in the middle of a twotimes larger array of zeros, before computing the phasecorrelation function, by inverse Fourier Transform of the final crosspower spectrum (GuizarSicairos, Thurman & Fienup 2008);

finds the position of the barycenter of the correlation peak for each subwindow, computes the coefficients of the interpolation law (assumed to be polynomial), and corrects the global wavelength table;

calculates the residuals of the wavelength shifts, using once again the phasecorrelation method on the wavelengthcorrected spectrum.
Name  Spec. Type  

Ret  G8IIIII  4780(230)    2.54  2.5(2) 
Ori  K0IIIb  4670(230)  2.20(2)  2.263  2.1(2) 
Col  K0III  4660(230)  2.48(3)  2.38  2.7(2) 
Hya  K0III  4660(230)    2.62   
Lib  K0III  4660(230)    2.31   
Sgr  K0III  4660(230)    2.50  2.4(2) 
Eri  K0.5IIIb  4600(220)  2.18(2)  2.35  2.5(2) 
Psc  K0.5III  4580(220)  2.00(2)  2.11  2.1(2) 
Scl  K1III  4510(220)  2.13(3)  2.039  2.2(2) 
HR 2113  K1.5III  4440(220)    2.48(2)  2.4(2) 
TrA  K1.5III  4440(220)  2.56(7)  2.43(2)   
Cet  K1.5III  4440(220)  3.44(4)  3.323  3.3(2) 
HR 3282  K2.5IIIII  4330(210)  2.54(4)  2.32(7)   
HR 2411  K3III  4260(210)  1.90(3)  1.76   
51 Hya  K3III  4260(210)  2.28(3)  2.23  2.4(2) 
For each observing night, the final wavelength table, associated with the working instrumental setting, is obtained by ensemble average of the corrected tables (of all exposures with all calibrators measured with the same instrumental setting), rejecting the tables showing wavelengthshift residuals higher than 8 (5 times the nominal spectral resolution). Figure 6 shows the final continuumcorrected spectrum, given by the spectralcalibration procedure, with the calibrator Lib. The model spectrum drawn with the dashed line is produced by the MARCS + Turbospectrum code.
To show the difference between the wavelengthcalibration provided by SPIDAST and by amdlib, we plot in Fig. 7 the corrections from the initial polynomial dispersion law of AMBER, given by the two softwares. While amdlib computes only a unique correction value applied over the whole spectrum, our method computes a polynomial correction (here parabolic from three subwindows), with subpixel deviations with respect to the polynomial law, caused by the final ensemble averaging of the corrected wavelength tables. The choice of the degree of the polynomial correction provided by SPIDAST allows the accuracy of the wavelength correction to be adapted to the scientific programme, reaching a subpixel precision if necessary.
4.2 Calibrating spectrointerferometric data
In order to assess and correct for the measurement defects on the science targets, proper interferometric calibration (not yet fully supported with amdlib) is needed. The Instrument Response Function (IRF), also called “transfer function” (Perrin 2003; van Belle & van Belle 2005; Boden 2007; Cruzalèbes et al. 2010), is derived from observations of calibrators as concomitant as possible to the ones of the science targets.
The standard calibration method is based on the assumption that a given “raw” (measured) quantity of any observable , is equal to the product of the “true” (intrinsic) quantity , multiplied by a global degradation factor (in other words the IRF). Such a factor cannot practically be modellised with sufficient reliability (numerous parameters and some of them remaining unknown).
To calibrate the quantity related to the science target, a second assumption applies, as soon as appropriate conditions (described later on) are satisfied and can be then determined by using the relation
(6) 
with self explanatory notations, stating that the degradation factor is identical for both the calibration source and the science target. This assertion is all the more valid that the required conditions are fullfilled. The relation here above yields an empirical determination of , via , relying on commonplace statistical tools. The main condition to satisfy is that the IRF is stable enough, and the standard approach to avoid too large variations of the IRF is to perform sequences, where observation of calibrator and science targets are interlaced and repeated several times, rapidly sampling the IRF. So, the actual stability requirement is that the IRF might be slowly variable over a given “calibratorsciencecalibrator” observation sequence. Besides, the calibrator and the science target should be close enough regarding time of observation, angular separation, and brightness (in the spectral domain of work), so that adjustments of the whole interferometer would not be significantly modified from one source to the other. In these conditions, only the atmospheric turbulence might cause substantial degradations of the measurements. For this latter case, servoloop systems are to be used to reduce turbulence effects.
Usually, calibrators are found in dedicated catalogs and the corresponding sources should meet the desired requirements. However, the number of calibrators is necessarily limited and somewhat depends on the type of science targets, so that the selected calibrators might stand out from the ideal conditions: pointlike source, small angular separation, brightness matching.
In practice, the limited choice of calibrators makes it necessary to accept “partially” resolved targets (van Belle & van Belle 2005), angular separations counted in degrees, and brightness mismatch amounting some units of magnitudes. For example, a bright calibrator is rarely found in a close angular vicinity of a given bright science target. Such a situation nevertheless remains convenient for the calibration procedure, since switching between science targets and calibrators is fast enough and does not affect the configuration of the interferometer, and its response as well.
This procedure starts with a first selection of calibrators via the SearchCal tool^{3}^{3}3www.mariotti.fr/searchcal_page.htm, created by the JMMC working group “Catalogue of calibration sources”, providing rough estimates for angular diameters, with comparatively large uncertainties. Following this first selection, the sources are searched in the calibrator catalogues of Bordé et al. (2002) and Mérand et al. (2005), in order to find angular diameters of better quality. Since some calibrators have not a diameter in the published catalogues, it is necessary to make our own determinations of those diameters (see Subsect. 4.2.1). Moreover, in order to control these determinations and to build a homogeneous set of diameters, we apply our own method to recalculate all the diameters for the selected calibrators. The good agreement found between our determinations and the ones of the Bordé’s Catalogue of Calibrator Stars for LongBaseline Stellar Interferometry (CCSL) (as shown in Table 2, and in Fig. 8), attests to the reliability of our specific determinations of angular diameters.
4.2.1 Determining calibrator angular diameters
To derive the IRF for each spectrointerferometric (SPI) observable, we firstly need to get a reliable estimate of the angular diameter of each calibrator, and its associated error. Fitting the stellar models of Kurucz (1979) on the absolutely calibrated spectrophotometric templates of Cohen et al. (1999), the CCSL contains 374 calibrators, with limbdarkened angular diameters ranging from 1 to 3 . Since our selected calibrators might have no angulardiameter estimate (one third of the calibrators of our dataset have no value found in the catalogues), or poorly estimated, we have added to SPIDAST various routines for determination of the angular diameter from indirect methods, presented in detail in Cruzalèbes et al. (2010). These routines compute the angular diameter:

by combining the linear diameter, deduced from absolute luminosity and effective temperature from the MorganKeenan (MK) spectral type, with the parallax;

by using experimental laws based on the interstellarcorrected colour index (surface brightness method);

by scaling synthetic spectra on broadband photometric measurements (infrared flux method);

by fitting synthetic spectra on infrared spectrophotometric measurements (see Sect. 6).
Table 2 gives the angular diameter values of our calibrators found in the CCSL, as well as our final estimates, derived from the fit of MARCS + Turbospectrum synthetic spectra on photometric or spectrophotometric measurements. In the last column, we also give the estimates found in the JMMC Stellar Diameter Catalogue (JSDC) of Lafrasse et al. (2010), with accuracies between 7 and 10% (Delfosse 2004; Bonneau et al. 2006).
Figure 8 shows the comparison of our estimates with the CCSL values. The mean difference between our estimates and the CCSL is 5%. The mean difference for the three targets with the available calibrated templates of Cohen (i.e. Ori, Scl, and Cet) is lower: 3.5%. This shows the satisfactory agreement between our results and the CCSL, which confirms the reliability of our approach for determination of the calibrator angular diameters, implemented in SPIDAST.
4.2.2 Modelling the calibrator visibility and triple product
Once we determined the angular diameter of each calibrator, and its associated uncertainty, we derive the IRF for each SPI observable, computing the ratio of the observable measured on the calibrator, to the synthetic observable, assumed to represent the “true” calibrator observable. Each synthetic SPI observable is derived from the synthetic coherent flux

the synthetic flux (Eq. 1) is given by the synthetic coherent flux computed for the null baseline (visibility unity);

the squared synthetic visibility (Eq. 2) by the ratio of the squared synthetic coherent flux to the synthetic squared flux;

the synthetic bispectrum (Eq. 3) by the product of the synthetic coherent flux for the three baselines composing each triplet;

and the synthetic triple product (Eq. 4) by the ratio of the synthetic bispectrum to the cubed synthetic flux.
Applying the vanCittertZernike theorem, the model coherent flux of each calibrator, assumed to emit a luminous intensity with circular symmetry, is derived from the Hankel transform of the source intrinsic spectral radiance (in )
(7) 
where is the stellar angular diameter, the impact parameter, i.e. the distance from the centre of the stellar disc ( at disc center), the length of the projected baseline of the (i,j)pair of apertures, and J the Bessel function of order zero.
Since compact photospheres are known to deviate from simple uniform discs (see Hajian et al. 1998, and references therein), radially decreasing from their photometric centre, we use the MARCS + Turbospectrum codes to produce models of CentretoLimb Variation (CLV) profiles, with input parameters derived from the MK spectral type. Figure 9 shows the – map of typical CLV profiles given by the model, where the reduced radial coordinate is given by = (Young 2003).
With “partially” resolved calibrators, the deviation from the UD model increases with the baseline length. Figure 10 shows the squared visibility and the real part of the triple product produced by the MARCS and the UD models, for the calibrator Ceti (=3.320.02 ), deduced from measurements obtained with the D0–I1–G1 baseline triplet (G1–I1=46 , D0–G1=69 , and D0–I1=79 ). In this example, the differences in visibility and triple product between the UD model and the LDMARCS model are larger than the uncertainty of the synthetic observables (derived by propagating the uncertainty on the angular diameter), hence the justification to use the LD model instead of the UD model.
4.2.3 Correcting for the degradations of the IRF
The main source of nonstationary instabilities which affects the fringe formation is the atmosphericphase turbulence (Roddier 1981). This issue is partially solved thanks to the use of fringetracking servoloops scanning the fringes with cycle periods smaller than the seeing coherence time, defined as the time over which the phase fluctuations remain coherent (Davis & Tango 1996; Kellerer & Tokovinin 2007). Thus, the residual instabilities are expected to cause only slow drifts of the IRF between observations of the science targets and their associated calibrators. In order to estimate the response at the time when each science target is measured, linear interpolation between the successive measurements on its calibrator(s) is legitimate (Perrin 2003), provided that the “calibratorsciencecalibrator” observation sequence, repeated several times, rapidly samples the evaluation of the IRF (Berger, Dumas & Kaüfer 2011).
To perform the calibration of the squared visibility and triple product data, SPIDAST applies the following procedure, working at spectral channel level:

Compute the measured IRF for each exposure, given by the ratio of the calibrator measured observable to the calibrator synthetic observable (at the same spectral resolution, for the same baseline or baseline triplet). Values and uncertainties are derived from the 4thorder approximation formulae of Winzer (2000). For triple product data, the approximation formulae apply separately to the real and imaginary parts of the complex ratio. Since the calibrator angulardiameter uncertainties are smaller than 10% (see Table 2), the uncertainty of a model observable is simply obtained from its first partial derivative w.r.t. the angular diameter, multiplied by the angular diameter uncertainty.

Average the IRF measurements over each OB. As weight associated with each exposure, we use the ratio of the associated weight (defined in Sect. 3.2) to the variance of the calibrator “raw” observable.

Determine the IRF at the mean time of each exposure obtained with the science target, from linear interpolation, or polynomial fit, on the averaged IRF measurements. Figure 11 shows the temporal variation of the measured and interpolated response in squared visibility and triple product, obtained during one observing night, with three successive calibrators.

Divide the observable measured on the science target for each exposure of a given OB, by each interpolated value of the IRF obtained for the same OB (using Winzer’s formulae). Applying a method similar to that used with the getCal tool of the Cal. Inst. Tech. (2008), a weight is associated with each calibrated ratio, which combines information on the observations of the calibrator and of the science target (science–calibrator angular separation and observationtime delay, coherence time, number of nonaberrant data).

Compute the weighted average over each OB of the ratios used in the calibration, which gives the final calibrated SPI observable (see Subsect. 4.2.4). The argument of the final calibrated triple product gives the final calibrated closure phase.
Although the calibration procedure used by SPIDAST is based on the standard calibration method, great care has been taken to provide reliable uncertainty estimates, thanks to the use of weights tracing the data quality at different levels of the processing. Thus, SPIDAST produces large uncertainties on final calibrated observables obtained with input data of poor quality. Estimating reliable error bars is of crucial importance to get calibrated data, usable for the fitting process of parametric models, even on data of uneven quality (see Sect. 6).
4.2.4 Comparing the final calibrated observables with amdlib
Although the calibration routine provided with the amdlib package is not yet officially validated (Malbet et al. 2011), we compute the differences in squared visibility and closure phase between SPIDAST and amdlib, and compare them with the associated uncertainties.
Figure 12 shows the results obtained on the science target Oph, with the baseline triplet A0–D0–H0, under good and poor seeing conditions. The histograms of the uncertainties show that the uncertainties in squared visibility estimated with amdlib are smaller than with SPIDAST (whatever the seeing), while it is the contrary for the uncertainties in closure phase.
To explain the disagreement between the uncertainties in squared visibility, we remind that SPIDAST takes into account instrumental and environmental conditions in the calculation of the uncertainties, while amdlib uses only unweighted variance. Thus, we consider as underestimated the uncertainties in squared visibility produced by amdlib.
Regarding the closure phase, it is more difficult to explain why the uncertainties from amdlib are larger than the ones from SPIDAST. Indeed, no details are found in Malbet et al. (2011), expected to describe the amdlib calculation. SPIDAST computes the uncertainties on the real and imaginary parts of the calibrated triple product, and derives the uncertainties on the final closure phase from them. A deep analysis based on the source code of the calibration routine provided by amdlib reveals that amdlib computes the uncertainties in the final closure phase, adding the uncertainties in closure phase of the science target and of the estimated transfer function. Since the closure phase is less stable than the triple product, the final uncertainties in closure phase handled with amdlib are larger than with SPIDAST.
To evaluate the agreement between the two procedures, we compare the differences between the final calibrated SPI observables they produce, to the associated uncertainties. Since the median differences are smaller than the median uncertainties (Fig. 12), we conclude that the agreement between the two softwares is satisfactory, SPIDAST giving more reliable uncertainties than amdlib.
5 Measuring the deviation from circular symmetry
The closure phase couples the phases of the Fourier transform of the source brightness distribution at the three spatial frequencies , , and , probed by the three baselines in the following way: . In the case of a centrosymmetrical brightness distribution onaxis, the phase of the Fourier transform is uniformly zero and is naturally equal to zero. If this source is offaxis, the phase is linear w.r.t. the spatial frequency, so that , and, subsequently, also here is equal to zero. Thus, a nonzero phase closure is a clear indication, at least qualitatively, for a deviation from circular symmetry.
In Sect. 5.2, we introduce a new parameter (called the centrosymmetry parameter), based on the triple product, more sensitive to deviation from centrosymmetry than the closure phase.
5.1 Global closure phase
The integration of the real and imaginary parts of the triple product , over the observation spectral band , leads to the global closure phase, defined by Ragland et al. (2006) and Tatebe et al. (2006) as
(8) 
A value of the global closure phase close to zero (less than 1 , for our sample) suggests a high degree of centrosymmetry of the brightness distribution, related to the observation spectral band. We compute the uncertainty on the global closure phase thanks to the directbootstrap method, using random sampling with replacement of the spectral dataset of the calibrated triple product. The values of the global closure phases in , and their associated uncertainties, are given in the third column of Table 3, for each science target.
Name  CSP (°)  (°)  SNR/SNR 

Oph  0.76(3)  0.16(4)  7.0 
Car  0.84(3)  0.22(5)  5.7 
L Pup  0.90(4)  0.39(5)  3.2 
Cet  0.98(4)  0.32(6)  4.6 
Ara  1.12(4)  0.22(6)  8.1 
TrA  1.21(5)  0.20(7)  8.3 
Hya  1.21(9)  0.16(11)  9.4 
TW Oph  2.94(5)  2.91(5)  1.1 
CE Tau  3.01(11)  0.85(17)  5.5 
Hyi  3.25(11)  2.83(13)  1.4 
CMa  4.41(197)  2.88(155)  1.2 
Lib  5.12(20)  1.90(33)  4.4 
Ret  5.16(37)  4.44(115)  3.6 
TX Psc  5.34(25)  2.17(32)  3.4 
Ori  8.10(82)  7.48(84)  1.1 
W Ori  10.58(46)  1.35(79)  13.3 
R Scl  11.63(59)  10.63(67)  1.2 
T Cet  27.88(51)  28.79(41)  0.8 
5.2 CentroSymmetry Parameter (CSP)
When the positive and negative values of along the spectral domain almost mutually compensate (as seen on Fig. 13 for TX Psc and W Ori), the spectral integration produces a nearly null global closure phase. As mentioned above, this can be taken as a hint for centrosymmetry. However, is clearly nonnull in some parts of the spectrum, and this is a hint for a deviation from centrosymmetry. To rule out the contradiction, and so as to consider this latter possibility, we introduce a new estimator, that we call the CentroSymmetry Parameter (CSP), similar to the global closure phase, but using instead the absolute value of in the numerator, and the modulus of in the denominator,
(9) 
As for the global closure phase, a small CSP value (less than 2 , for our sample) suggests a high degree of centrosymmetry. Significantly high CSP values, as judged from their uncertainties (bootstrap method again), require to use asymmetric brightness distribution models, for the fitting on the SPI data. The values of the global CSP in the band, and their uncertainties, are given in the second column of Table 3.
To illustrate the difference between the global closure phase and the CSP, we plot in Fig. 13 three quantities, for 3 scientific targets showing a global closure phase close to zero,

in the panels on the left is shown the imaginary part of the triple product , divided by the spectral mean^{4}^{4}4defined as the integral w.r.t. the wavelength, divided by the spectral bandwidth of the real part . Note that the spectral mean of this ratio is equal to the tangent of the global closure phase, as defined in Eq. (8);

in the panels forming the central column is shown the absolute value of the imaginary part of the triple product , divided by the spectral mean of the modulus . Note that the spectral mean of this ratio is equal to the sine of the global CSP, as defined in Eq. (9);

in the panels on the right is shown the absolute value of the real part of the triple product , divided by the spectral mean of the modulus .
The three panels on the left show three typical behaviours for :

uniformly close to zero ( Oph);

decreasing symmetrically around zero (TX Psc);

increasing symmetrically around zero (W Ori).
For each of these situations, the integration over the spectrum produces a nearly null result (hint for centrosymmetry). However the CSP allows to detect deviations from centrosymmetry (see Fig. 14).
If the real and imaginary parts of the triple product are wavelengthindependent, for the observation spectral range, one can show that . If not, no analytical relation linking the two quantities can be derived from Eqs. (8) and (9), because of the integrals in the numerators and the denominators. Figure 14 shows the CSP in displayed versus the global closure phase in the same spectral band, for each science target, choosing the OB which gives the smallest uncertainty on the CSP. Almost all the stars fall along the diagonals of the (, CSP) diagram, which trace the relation CSP = . One (W Ori), however, is flagged as asymmetrical with the CSP indicator, but not with . This is precisely why we favour the CSP over , as discussed above. T Cet has the largest values for CSP and as well, indicating large deviation from centrosymmetry (may be due to strong asymmetries at the wavelengths of the CO bands or even suggestive of the presence of a binary companion). The second largest CSP value pertains to R Scl, recently revealed as a wide binary by Maercker et al. (2012), using the ALMA array at millimetric wavelengths. In addition, the last column of Table 3 shows that, except for T Cet, the SNR for the CSP is comparable or larger than the SNR for the global closure phase.
6 Fitting parametric chromatic models
In order to interpret the final calibrated SPI observables, SPIDAST provides a fitting routine, based on the modified gradientexpansion algorithm, very similar to the algorithm invented by Levenberg (1944), and improved by Marquardt (1963). The fit applies on any calibrated SPI measurements (to be chosen between visibility, flux, coherent flux, bispectrum, triple product, and closure phase), using a library of singlecomponent or composite parametric chromatic models, characterised by the Fourier spectrum of their intensity distribution, and the associated firstorder partial derivatives w.r.t. the model parameters. Figure 15 shows examples of visibilityfitting results obtained with SPIDAST, using two models (UD and LDMARCS). This fitting routine, applied on spectrophotometric data, is used to determine the calibrator angular diameters from the MARCS + Turbospectrum synthetic spectra (see Subsect. 4.2.1).
Note that the LITpro software^{5}^{5}5www.jmmc.fr/litpro_page.htm, developed by the JMMC working group “Model fitting”, uses a set of elementary geometrical and centertolimb darkening functions, as well (TallonBosc et al. 2008). However, it does not offer, contrary to SPIDAST, the possibility to fit stellardisk models with synthetic tabulated radiance data (or exitance, for fits on spectrophotometric measurements), such those produced by the MARCS + Turbospectrum code.
Since the uncertainties of the final calibrated data are not normally distributed, the covariance matrix, that comes out of the chisquared fit, cannot be used to infer the parameter uncertainties (Press et al. 2007). To compute the uncertainties, SPIDAST uses the residualbootstrap method, described in detail in Cruzalèbes et al. (2010):

“synthetic” datasets are produced from random sampling with replacement of the Pearson residuals (difference between calibrated and fitted values divided by the uncertainty of the observed value), added to the initial fitted values;

fitting the model on these new datasets produces a set of chisquared minima, following a probability distribution, from which we extract the boundaries of the confidence interval, with a given confidence level;

the parameter values associated with these chisquared boundaries give the upper and lower limits of the parameter estimates, leading to asymmetric uncertainties.
The whole procedure allows the computing of reliable estimates of the parameter uncertainties.
7 Conclusion
In the present paper, we introduce the new SPIDAST software, developed since 2006 with the aim to reduce, calibrate, and interpret the visibility and triple product measurements obtained with the VLTI/AMBER facility. SPIDAST contains a whole set of modules, which can be launched separately or in an automatic batch file, and summarised herafter:

The raw data reduction used by SPIDAST computes the weighted average of nonaberrant data, at spectralchannel level, providing estimates of the SPI observables using an automated procedure, while the method presently used with amdlib selects the “best” frames according to a quality threshold determined a posteriori, after several trials.

The wavelength calibration procedure performed by SPIDAST provides spectral shifts following a polynomial law, tracing them at the channel level. This is done by computing the crosscorrelation of the measured spectra of the calibrators with their synthetic spectra produced by the MARCS model, while amdlib only provides a constant spectral shift over the spectrum, from the correlation with the telluric lines.

For selected calibrators not included in the calibrator catalogues, SPIDAST provides several routines for estimation of the angular diameter with indirect methods, the most accurate being the fit of stellaratmosphere model spectra given by MARCS on spectrophotometric data. The calibrator synthetic observables are derived using CLV functions produced by MARCS.

To obtain an accurate interferometric calibration (via an automated procedure), SPIDAST: (1) divides the calibrator raw data with the associated synthetic observables in each spectral channel, which gives the instrumental response function in squared visibility and triple product; (2) interpolates or fits the response at the time of each exposure on the science target; (3) divides the science raw data with the interpolated/fitted response, which gives the science calibrated observables for each exposure; (4) computes their weighted average over each OB. At each processing step, the uncertainties are computed thanks to the bootstrap method applied on the weighted means.

Using the real and imaginary parts of the calibrated triple product, SPIDAST measures the deviation from centrosymmetry of the brightness distribution of each science target in the observation spectral band, thanks to a new parameter, more sensitive to asymmetries than the global closure phase.

Finally, SPIDAST proposes a complete fitting tool, using a set of parametric and chromatic models, and accepting input tables of flux/intensity synthetic data.
Such a careful calibration process of SPI data is a crucial step for their trustworthy astrophysical exploitation, which is the topic of two associated papers (Cruzalèbes et al. 2013b, c). Parameter extraction using nonlinear fits of source model, as well as aperturesynthesis image reconstruction, need reliable estimates of the calibrated observables, with robust uncertainties. Our reduction, calibration, and fitting routines also apply to any other spectral datasets, including spectroscopic data. We made the SPIDAST software public^{6}^{6}6https://forge.oca.eu/trac/spidast: the source code of any program of our software suite can be obtained by sending an email to the first author of the present paper.
Acknowledgments
The authors thank the ESOParanal VLTI team for supporting their AMBER observations, especially the night astronomers A. Mérand, G. Montagnier, F. Patru, J.B. Le Bouquin, S. Rengaswamy, and W.J. de Wit, the VLTI group coordinator S. Brillant, and the telescope and instrument operators A. and J. Cortes, A. Pino, C. Herrera, D. Castex, S. Cerda, and C. Cid. They also thank the Programme National de Physique Stellaire (PNPS) for supporting part of this collaborative research. A.J. is grateful to B. Plez, K. Eriksson, and T. Masseron for their ongoing support on the use of the MARCS code. S.S. was (partly) supported by the Austrian Science Fund through FWF project P19503N16; A.C. is postdoctoral fellow from F.R.S.FNRS (Belgium; grant 2.4513.11); E.P. is supported by PRODEX. The data were partially reduced using the publicly available data reduction software package amdlib, kindly provided by the JeanMarie Mariotti Center. This study used the SIMBAD and VIZIER databases at the CDS, Strasbourg (France), and NASAs ADS bibliographic services.
References
 Alvarez & Plez (1998) Alvarez R., Plez B., 1998, A&A, 330, 1109
 Berger et al. (2011) Berger J.P., Dumas C., Kaüfer A., 2011, VLT Paranal Science Operations AMBER User Manual. ESO, Garching b. München, Germany
 Boden (2007) Boden A. F., 2007, NewA, 51, 617
 Bonneau et al. (2006) Bonneau D., Clausse J.M., Delfosse X., et al., 2006, A&A, 456, 789
 Bordé et al. (2002) Bordé P., Coudé du Foresto V., Chagnon G., Perrin G., 2002, A&A, 393, 183
 Cal. Inst. Tech. (2008) Cal. Inst. Tech. 2008, The getCal Interferometric Observation Planning Tool Suite. NASA Exoplanet Science Inst., 2.10.7 edn
 Cohen et al. (1999) Cohen M., Walker R. G., Carter B., et al., 1999, AJ, 117, 1864
 Cohen et al. (2003) Cohen M., Wheaton W. A., Megeath S. T., 2003, AJ, 126, 1090
 Cruzalèbes et al. (2013b) Cruzalèbes P., Jorissen A., Rabbia Y., et al., 2013b, in preparation
 Cruzalèbes et al. (2013c) Cruzalèbes P., Jorissen A., Rabbia Y., et al., 2013c, in preparation
 Cruzalèbes et al. (2010) Cruzalèbes P., Jorissen A., Sacuto S., Bonneau D., 2010, A&A, 515, A6
 Cruzalèbes et al. (2008) Cruzalèbes P., Spang A., Sacuto S., 2008, in Kaufer A., Kerber F., eds, The 2007 ESO Instrument Calibration Workshop Calibration of AMBER Visibilities at Low Spectral Resolution. Springer, pp 479–482
 Davis & Tango (1996) Davis J., Tango W., 1996, PASP, 108, 456
 Decin & Eriksson (2007) Decin L., Eriksson K., 2007, A&A, 472, 1041
 Decin et al. (2004) Decin L., Morris P. W., Appleton P. N., et al., 2004, ApJSS, 154, 408
 Delfosse (2004) Delfosse X., 2004, Technical report, Stellar diameter estimation from photospheric indices. JMMC
 Domiciano de Souza et al. (2008) Domiciano de Souza A., Bendjoya P., Vakili F., Millour F., Petrov R. G., 2008, A&A, 489, L5
 Efron (1979) Efron B., 1979, Ann. Statist., 7, 1
 Efron & Tibshirani (1993) Efron B., Tibshirani R., 1993, An Introduction to the Bootstrap, CRC Press edn. Vol. 57 of Monographs on Stat. and Appl. Prob., Chapman and Hall
 GuizarSicairos et al. (2008) GuizarSicairos M., Thurman S. T., Fienup J. R., 2008, Opt. Lett., 33, 156
 Gustafsson et al. (2008) Gustafsson B., Edvardsson B., Eriksson K., et al., 2008, A&A, 486, 951
 Hajian et al. (1998) Hajian A. R., Armstrong J. T., Hummel C. A., et al., 1998, ApJL, 496, 484
 Hoaglin et al. (1983) Hoaglin D. C., Mosteller F., Tukey J. W., 1983, Understanding Robust and Exploratory Data Analysis
 Kellerer & Tokovinin (2007) Kellerer A., Tokovinin A., 2007, A&A, 461, 775
 Kraus et al. (2009) Kraus S., Weigelt G., Balega Y. Y., et al., 2009, A&A, 497, 195
 Kurucz (1979) Kurucz R. L., 1979, ApJSS, 40, 1
 Lafrasse et al. (2010) Lafrasse S., Mella G., Bonneau D., et al., 2010, in W. C. Danchi et al. ed., Optical and infrared interferometry II Vol. 7734 of SPIE Conf. Ser., Building the ’JMMC Stellar Diameters Catalog’ using SearchCal. San Diego, CA
 Levenberg (1944) Levenberg K., 1944, Q. Appl. Math., 2, 164
 Maercker et al. (2012) Maercker M., Mohamed S., Vlemmings W. H. T., et al., 2012, Nature, 490, 232

Malbet et al. (2011)
Malbet F., Millour F., Absil O., Duvert G., 2011, User Manual 
Issue 3.0.3 JMMCMAN27200001, AMBER Data Reduction Software,
www.jmmc.fr/doc/approved/JMMCMAN27200001.pdf
. JMMC  Marquardt (1963) Marquardt D. W., 1963, Indian J. Indust. Appl. Math., 11, 431
 MartíVidal et al. (2011) MartíVidal I., Marcaide J. M., Quirrenbach A., et al., 2011, A&A, 529, A115
 Mérand et al. (2005) Mérand A., Bordé P., Coudé Du Foresto V., 2005, A&A, 433, 1155
 Mérand et al. (2010) Mérand A., Stefl S., Bourget P., et al., 2010, in W. C. Danchi et al. ed., Optical and infrared interferometry II Vol. 7734 of SPIE Conf. Ser., Perspectives for the AMBER Beam Combiner. San Diego, CA
 Millour et al. (2008) Millour F., Petrov R., Malbet F., et al., 2008, in Kaufer A., Kerber F., eds, The 2007 ESO Instrument Calibration Workshop AMBER on the VLTI: data processing and calibration issues. Springer, pp 461–470
 Millour et al. (2007) Millour F., Petrov R. G., Chesneau O., et al., 2007, A&A, 464, 107
 Ohnaka et al. (2011) Ohnaka K., Weigelt G., Millour F., et al., 2011, A&A, 529, A163
 Perrin (2003) Perrin G., 2003, A&A, 400, 1173
 Petrov et al. (2007) Petrov R. G., Malbet F., Weigelt G., et al., 2007, A&A, 464, 1
 Plez (2003) Plez B., 2003, in U. Munari ed., GAIA Spectroscopy: Science and Technology Vol. 298 of ASP Conf. Ser., Cool star atmospheres and spectra for GAIA: MARCS models. pp 189–198
 Plez (2012) Plez B., 2012, Astrophys. Source Code Libr., p. 1205.004
 Press et al. (2007) Press W. H., Teukolsky S. A., Vetterling W. T., et al., 2007, Numerical Recipes 3rd Edition  The Art of Scientific Computing. Cambridge Univ. Press
 Ragland et al. (2006) Ragland S., Traub W. A., Berger J., et al., 2006, ApJL, 652, 650
 Rix et al. (2004) Rix S. A., Pettini M., Leitherer C., et al., 2004, ApJL, 615, 98
 RobbeDubois et al. (2007) RobbeDubois S., Lagarde S., Petrov R. G., et al., 2007, A&A, 464, 13
 Roddier (1981) Roddier F., 1981, Prog. Opt., 19, 281
 Stefl et al. (2011) Stefl S., Mérand A., Le Bouquin J.B., et al., 2011, in The Origin and Fate of the Sun: Evolution of Solarmass Stars Observed with High Angular Resolution Problems of the AMBER wavelength calibration. to be published, Garching b. München, Germany
 TallonBosc et al. (2008) TallonBosc I., Tallon M., Thiébaut E., et al., 2008, in Optical and Infrared Interferometry Vol. 7013 of SPIE Conf. Ser., LITpro: a model fitting software for optical interferometry. Marseille, France, pp 70131J1–70131J9
 Tatebe et al. (2006) Tatebe K., Chandler A. A., Hale D. D. S., Townes C. H., 2006, ApJL, 652, 666
 Tatulli et al. (2007) Tatulli E., Millour F., Chelli A., et al., 2007, A&A, 464, 29
 van Belle & van Belle (2005) van Belle G. T., van Belle G., 2005, PASP, 117, 1263
 Vera & Torres (2008) Vera E., Torres S., 2008, in EURASIP ed., 16th ESP Conf. (EUSIPCO 2008) Subpixel accuracy analysis of phase correlation registration methods using aliased imagery. Lausanne, Switzerland
 Weigelt et al. (2011) Weigelt G., Grinin V. P., Groh J. H., et al., 2011, A&A, 527, A103
 Weigelt (1977) Weigelt G. P., 1977, Opt. Commun., 21, 55
 Winzer (2000) Winzer P., 2000, Rev. Sci. Instrum., 71, 1447
 Wittkowski et al. (2008) Wittkowski M., Boboltz D. A., Driebe T., et al., 2008, A&A, 479, L21
 Young (2003) Young J., 2003, in Perrin G., Malbet F., eds, Observing with the VLTI Vol. 6 of EAS Publ. Ser., From Visibilities to Science with Simple Models. pp 181–190