5 Model fitting

  The aperture synthesis imaging attempts to reconstruct the brightness distribution on the sky, in principle without making any a priori assumption about the object. This is excellent for exploring cases where the nature of the object is uncertain, e.g. concerning the number of resolved components in a multiple star or their approximate positions. The method is less useful for accurate quantitative evaluation, in particular because the offsets in parallax and proper motion merely produce a blurring of the image. Since the objects of interest here consist of a small number of point sources, direct modeling of the TD in terms of simple superposed signal components is usually possible. Such model fitting provides the most direct and accurate estimates of specific object parameters such as the trigonometric parallax or orbital elements. In this section we outline the fitting procedure and give an example of its practical realization by means of a publicly available computer program.

5.1 General method

  Let us assume that the object consists of n point sources with intensities A_j and positions x_j, y_j relative the reference point ($j=1 \dots n$). In general A_j, x_j, y_j vary with time, and so may be different for the different transits of the same object. Given a specific model of the object we express A_j, x_j, y_j as functions of time t and a set of model parameters $\vec{a}$. For instance, in the case of a non-variable orbital binary, $\vec{a}$ would consist of 15 parameters, viz. the five astrometric parameters of the mass centre, the magnitude of each component, the mass ratio, and seven elements for the relative orbit. Generally speaking, the object model is thus completely specified by n and the functions $A_j(t,\vec{a})$, $x_j(t,\vec{a})$, $y_j(t,\vec{a})$ for $j=1 \dots n$. In the equations below we suppress, for brevity, the explicit dependence on t and $\vec{a}$.

For a given transit the expected signal is modeled as the sum of the signals from the individual components, using Eqs. (4) and (5). Thus,

$$I_k = \sum_{j=1}^n A_j [ 1 + M_1\cos(p_k+f_xx_j+f_yy_j) \\ \quad+ M_2\cos 2(p_k+f_xx_j+f_yy_j) ].$$ (13)

Expanding the trigonometric functions and equating the terms with those in Eq. (1) yields

$$b_1 = \phantom{-M_1}\sum_j A_j, \\ \nonumber$$

$$b_2 = \phantom{-}M_1\sum_j A_j \cos(f_xx_j+f_yy_j), \\ \nonumber$$

$$b_3 = -M_1\sum_j A_j \sin(f_xx_j+f_yy_j), \\ \nonumber$$

$$b_4 = \phantom{-}M_2\sum_j A_j \cos 2(f_xx_j+f_yy_j), \\ \nonumber$$

 

$$b_5 = -M_2\sum_j A_j \sin 2(f_xx_j+f_yy_j).$$ (14)

Recall that A_j, x_j, y_j depend on the model parameters $\vec{a}$.The general procedure is then to adjust $\vec{a}$ in such a way that, for the whole set of transits, the calculated signal parameters b₁-b₅ from Eq. (14) agree, as well as possible, with the observed values. The adjustment may use the weighted least-squares method, using the standard errors of the observed signal parameters to set the weights; but other (and more robust) metrics can also be used. In general the problem can be formulated as a constrained minimization problem in the multi-dimensional model parameter space.

The trigonometric functions in Eq. (14) mean that the signal parameters b_i depend in a highly non-linear manner on the model parameters which affect x_j and y_j. For instance, in terms of a displacement of one of the point sources, the effect on b₄ and b₅ is approximately linear only for displacements less than about $1/2f \simeq 0.1$ arcsec, corresponding to 1 rad change in the modulation phase. In the aperture synthesis imaging this non-linearity is manifest in the complex structure of the "dirty beam" (Fig. 4) at all spatial scales larger than about 0.1 arcsec. Additional non-linearities in the complete object model may result from the geometrical description of the source positions, e.g. in terms of orbital elements.

The non-linearity of the object model has two important consequences for the model fitting. Firstly, it is usually necessary to use a non-linear, iterative adjustment algorithm, such as the Levenberg-Marquardt method (Press et al. 1992). Secondly, a good initial guess of the model parameters is usually required. In particular the parameters directly affecting the positions of the point sources need to be specified to within (what corresponds to) a few tenths of an arcsec. Without a good initial guess, the adjustment algorithm is likely to converge on some local minimum, typically resulting in positional errors of (approximately) an integer number of grid periods. The correct solution, corresponding to the global minimum, may in principle always be found through sufficiently extensive searching of the parameter space. Alternatively, sufficiently good initial guesses of the point source positions can often be obtained from the aperture synthesis imaging.

Various least-squares model fitting procedures were used for the reduction of double and multiple stars during the construction of the Hipparcos Catalogue (see Mignard et al. 1995 and references therein). The double-star processing of the NDAC data reduction consortium (Söderhjelm et al. 1992) essentially used the technique outlined above, taking the so-called Case History Files (a precursor to the TD) as input.

Perhaps the greatest potential of the TD lies in the possibility to combine the Hipparcos data with independent observations from other instruments and epochs. For instance, full determination of a binary orbit generally requires data covering at least a whole period. Ground-based speckle observations can sometimes provide this, constraining the geometry of the relative orbit much better than the Hipparcos data alone, and in turn leading to a better-determined space parallax. In some favourable cases the location of the mass centre in the relative orbit (and hence the mass ratio) can be determined (Söderhjelm et al. 1997; Söderhjelm 1999).

One complication of the Hipparcos double star processing has been the wide variety of applicable object models, and the consequent need to experiment and interact with the solutions. This process may be much facilitated by using general and flexible software for the model fitting, rather than highly specialized routines. An example of this is given below.

5.2 Model fitting using GaussFit

GaussFit (Jefferys et al. 1988a, 1988b) is a general program for the solution of least squares and robust estimation problems, developed as a platform to facilitate astrometric reduction of data from the Hubble Space Telescope. It is written in the C programming language and may thus be run under a variety of operating systems. In this section we outline the use of GaussFit for model fitting to the TD, again using the binary HIP 97237 as illustration.

GaussFit was used by Söderhjelm (1999) in a systematic re-examination of the solutions for several double and multiple objects, through a combination of TD with ground-based observations. Although not illustrated in the example below, the introduction of additional data (e.g. relative positions from speckle observations) is quite straightforward by means of GaussFit.

To run GaussFit, the user must supply several input files. During execution these files are read (and sometimes modified) by GaussFit, and additional output files generated. For application to the TD model fitting the following input files are required.

• The data file: this contains the observational data, in our case the TD. A special program (td2uv.f) is available (Sect. 6) to extract the TD for a given HIP number and format them as required by GaussFit. The resulting data file consists of 16 columns and one data line per transit. The columns contain a sequential number for the transit, the target position index (JT1), the time of the transit, the spatial frequencies f_x, f_y, $f_{\rm p}$, the signal parameters b₁-b₅ and their variances. The header of the data file defines the name of the variable associated with each column.
• The model file (cf. Fig. 8): this is a mathematical description of the object model written in the GaussFit programming language. This language is modeled on C, but includes some specific constructs. For instance, the declaration of variables distinguishes between "observations'' (input data with random errors that need to be taken into account in the fitting), "data'' (error-free input data), "parameters'' (to be adjusted by the program), and ordinary "variables''. The special function $\tt import()$ reads one line of data from the data file. The function $\tt export(x)$ sends the equation of condition x=0 to the estimation algorithm, taking into account the uncertainties of the observational data that went into calculating x.
• The parameter file: this contains the initial guesses of all the model parameters to be estimated. On output it contains the estimated parameter values and estimated errors.
• The environment file contains general information needed for the model fitting, such as the names of the data, parameter and output files; the type of estimation algorithm to be used (standard least squares or a robust method), and stopping rules for the iterations.

The reader is referred to the GaussFit User's Manual (Jefferys et al. 1988b) for detailed information.

  

Figure 8: An example of a double-star model defined in the GaussFit programming language

  

Figure 9: Part of the GaussFit output (slightly edited) obtained while fitting the double star model in Fig. 8 to the TD for HIP 97237

Figure 8 is an example of a GaussFit model file. It describes a binary with a fixed positional offset between the components (i.e. a long-period binary). The model parameters are thus the astrometric parameters of the primary relative to the reference point (${\tt x1}=\Delta\alpha*_1$, ${\tt y1}=\Delta\delta_1$, ${\tt par}=\Delta\pi$,${\tt xdot}=\Delta\mu_{\alpha*}$, ${\tt ydot}=\Delta\mu_\delta$), the position of the secondary relative to the reference point (${\tt x2}=\Delta\alpha*_2$, ${\tt y1}=\Delta\delta_2$); and the intensities of the components, ${\tt A1}=A_1$, ${\tt A2}=A_2$.The components are assumed to have the same parallax and proper motion. The expressions within the $\tt export()$ functions are easily recognized as the equations of condition, Eq. (14), written in terms of the model parameters. The five $\tt export()$ statements are divided among two $\tt import()$ loops (which means that the data file is forced to be read twice in each iteration): the reason is that GaussFit in its standard distribution version cannot handle more than four simultaneous equations of condition.

The model in Fig. 8 was applied to the TD of HIP 97237, using as starting approximation (3600, 3600, 0, 0, 0, 0.04, 4000, 4300, 0.02) for the variables in the parameter list (cf. Sect. 4.3). The "fair'' metric with an asymptotic relative efficiency of 0.95 was used for robust estimation of the parameters (Jefferys et al. 1988a). Part of the output file, containing the results of the final (10th) iteration, is shown in Fig. 9. It should be noted that the estimated standard errors (sigma values) given in the output file have already been scaled by $(\chi^2/\nu)^{1/2}$, using the chi-square ($\chi^2$) and degrees of freedom ($\nu$) given at the end of the file. Adding the results of the model fitting to the reference point data (Sect. 4.3) and using the magnitude conversion formula $Hp=-2.5\log(A/K)$ we obtain the following estimated parameters of the binary HIP 97237 (ICRS, epoch J1991.25): 0pt

$$\alpha_1 = 296.43961803~\mbox{deg}\pm \phantom{0}4.41~\mbox{mas}, \\ \nonumber$$

$$\alpha_2 = 296.43974889~\mbox{deg}\pm \phantom{0}8.65~\mbox{mas}, \\ \nonumber$$

$$\delta_1 = \phantom{0}27.12836895~\mbox{deg}\pm \phantom{0}7.00~\mbox{mas}, \\ \nonumber$$

$$\delta_2 = \phantom{0}27.12856546~\mbox{deg}\pm 13.37~\mbox{mas}, \\ \nonumber$$

$$\pi = \phantom{-00}80.74 \pm 8.80~\mbox{mas} \\ \nonumber$$

$$\mu_{\alpha*} = \phantom{00}\mbox{$-74.94$}\pm 5.24~\mbox{mas~yr}^{-1}, \\ \nonumber$$

$$\mu_\delta = -1203.93 \pm 7.68~\mbox{mas~yr}^{-1}, \\ \nonumber$$

$$Hp_1 = 12.88 \pm 0.01~\mbox{mag}, \\ \nonumber$$

$$Hp_2 = 13.60 \pm 0.02~\mbox{mag}. \nonumber$$

These data are in reasonable agreement with the values derived by Söderhjelm (1999) in an orbital solution combining the TD with ground-based speckle observations.