% mn2esample.tex
%
% v2.1 released 22nd May 2002 (G. Hutton)
%
% The mnsample.tex file has been amended to highlight
% the proper use of LaTeX2e code with the class file
% and using natbib cross-referencing. These changes
% do not reflect the original paper by A. V. Raveendran.
%
% Previous versions of this sample document were
% compatible with the LaTeX 2.09 style file mn.sty
% v1.2 released 5th September 1994 (M. Reed)
% v1.1 released 18th July 1994
% v1.0 released 28th January 1994

\documentclass[useAMS,usenatbib,onecolumn]{mn2e}
\usepackage{graphicx}

% If your system does not have the AMS fonts version 2.0 installed, then
% remove the useAMS option.
%
% useAMS allows you to obtain upright Greek characters.
% e.g. \umu, \upi etc.  See the section on "Upright Greek characters" in
% this guide for further information.
%
% If you are using AMS 2.0 fonts, bold math letters/symbols are available
% at a larger range of sizes for NFSS release 1 and 2 (using \boldmath or
% preferably \bmath).
%
% The usenatbib command allows the use of Patrick Daly's natbib.sty for
% cross-referencing.
%
% If you wish to typeset the paper in Times font (if you do not have the
% PostScript Type 1 Computer Modern fonts you will need to do this to get
% smoother fonts in a PDF file) then uncomment the next line
% \usepackage{Times}

%%%%% AUTHORS - PLACE YOUR OWN MACROS HERE %%%%%

% OBJECTS MACROS
%%%%%%%%%%%%%%%%
\def\hd{HD~152}

% GREEK LETTERS
%%%%%%%%%%%%%%%%
\def\a{\alpha}
\def\b{\beta}
\def\g{\gamma}
\def\t{\theta}
\def\o{\omega}
\def\s{\sigma}
\def\sb{\sigma^2}

% UNITS MACROS
%%%%%%%%%%%%%%
\def\kms{km\,s$^{-1}$}

% BIB MACROS
%%%%%%%%%%%%%
\def\aap{A\&A}
\def\aaps{A\&AS}
\def\aj{AJ}
\def\apj{ApJ}
\def\apjs{ApJS}
\def\apjl{ApJ}
\def\aplett{Astrophys. Lett.}
\def\mnras{MNRAS}
\def\nat{Nature}
\def\pasp{PASP}

% citation ALIASES
%%%%%%%%%%%%%%%%%%
\defcitealias{WHS67}{WHS67}


% STATS MACROS
%%%%%%%%%%%%%%%
\def\x{ \overline{x} }
\def\xb{ \overline{x^2} }
\def\y{ \overline{y} }
\def\yb{ \overline{y^2} }
\def\xy{ \overline{xy} }


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title[]{The Li\`ege SB2 Orbital Solution Algorithm}
\author[H.~Sana]{H.~Sana$^{1}$\thanks{E-mail:
sana@astro.ulg.ac.be}\\
$^{1}$Institut d'Astrophysique, Li\`ege University, All\'ee du 6 Ao\^ut 17, B-4000 Li\`ege, Belgium\\}
\begin{document}

\date{Accepted 1988 December 15. Received 1988 December 14; in original form 1988 October 11}

\pagerange{\pageref{firstpage}--\pageref{lastpage}} \pubyear{2002}

\maketitle

\label{firstpage}

\begin{abstract}
We present a semi-analytical algorithm to derive self-consistent SB2 orbital solutions for both components of a spectroscopic binary system. The method combines an orthogonal regression technique and an adaptation of the Wolfe et al. (1967) code, which allows to derive SB1 orbital solutions with only a limited {\it a priori} assumption on the period value and none on the other orbital parameters. The present algorithm preserves much of the advantage of the original Wolfe et al.\ method, in particular no fastidious exploration of the 8-dimension parameter space is required. As an illustrating case, the method is applied to the well known HD~152248 O+O colliding wind binary.
\end{abstract}

\begin{keywords}
Binaries: spectroscopic -- Stars: fundamental parameters -- Methods: numerical
\end{keywords}

\section{Motivation and original idea}


In an eccentric binary, the radial velocity $v$ of a component can be written as a function of the true anomaly $\t$ through the equation \citep[e.g.][]{Ait35} : 
\begin{equation}
 v(\t)=\gamma + K \sin i (\cos(\o+\t)+e \cos\t) \label{eq: rv1}
\end{equation}
or similarly 
\begin{equation} 
v(\t)=\gamma + \sin i \sin \o \frac{d}{dt}(r \cos \t) + \sin i \cos \o \frac{d}{dt}(r \sin \t) \label{eq: rv2}
\end{equation}
where $r$ is the distance to the center of mass. The other notations used are defined in Table~\ref{tab: orbparam}.\\

For an SB1 system, thus for which only the radial velocity (RV) of one component out of the two can be measured, we (i.e.\ the Li\`ege team) usually adopt the \citet[hereafter \citetalias{WHS67}]{WHS67} semi-analytical method. Briefly the method is that of Wilsing-Russell \citep{Wilsing1894, Russell1902}, followed by a differential correction (D.C.). From first guess of the period, the algorithm uses a development of  $r \cos \t$ and $r \sin \t$ in harmonic series of the mean anomaly, that uses Bessel functions of the eccentricity. The Bessel functions are then developed in asymptotic series of the eccentricity to obtain a set of equations that linearly depends on parameters $C_k$'s and $S_k$'s. These latter are themselves functions of the orbital parameters. Thanks to the linearity of the system, the $C_k$'s and $S_k$'s can be derived by a least square technique. In a second step, the non linear equations giving the orbital elements as a function of the $C_k$'s and $S_k$'s are solved iteratively. For the case of moderate or large eccentricities, the preliminary orbit so obtained can be improved  by means of a D.C. method. At this stage, the period can also be adjusted. \\


\begin{table}
\centering
\caption{Adopted notations for the orbital parameters of a massive binary system.}
\label{tab: orbparam} 
\begin{tabular}{r l}
\hline
\hline
 Notation & Description\\
\hline
$P$            & Orbital period \\
$i$            & System inclination \\
$e$            & Orbit eccentricity [0,1]\\
$T_0$          & Time of periastron passage ($e\ne0$) or time of primary conjunction ($e=0$) \\
$\omega$       & Longitude of the periastron ($e\ne0$)\\
$K_{1,2}$      & Radial velocity curve semi-amplitude for component 1, 2\\
$\gamma_{1,2}$ & Apparent systemic radial velocity associated with component 1, 2\\
$a_{1,2}$      & Semi-major axis of the orbit of component 1, 2. We also note $a=a_1+a_2$\\
$M_{1,2}$      & Mass of component 1, 2 \\
$q=M_1/M_2$    & Binary component mass ratio \\
$R_\mathrm{RL_{1,2}}$ & Roche lobe radius associated with component 1, 2 \\
\hline
\end{tabular}
\end{table}


In an SB2 system, one can write Eq. \ref{eq: rv1} or Eq. \ref{eq: rv2} for each of the two components of the system, so that e.g.:

\begin{equation} \left\{ \begin{array}{lll}
 v_1(\t)&=&\gamma_1 + K_1 \sin i (\cos(\o_1+\t)+e \cos\t) \\
 v_2(\t)&=&\gamma_2 + K_2 \sin i (\cos(\o_2+\t)+e \cos\t) \end{array} \right. \label{eq: rv1ps} 
\end{equation}

 Because of the coupling between these two equations (through the variables $e$, $\o_2=\pi+\o_1$ and $P$), the \citetalias{WHS67} method can not directly be used to provide a common orbital solution. Alternatively one could consider both the primary and secondary data sets as independent SB1 sets and one could derive a separate solution for both components. However, as seen in Table \ref{tab: hd248}, this provides two  distinct values for the common parameters, i.e.\ those parameters that induce a coupling between the Eqs.~\ref{eq: rv1ps}. 

 Of course one could implement appropriate numerical methods to solve
 non-linear systems of equations. However, to our knowledge, these
 methods rely on an accurate first guess of the solution or on
 complicated parameter space exploration algorithms (such as the
 so-called genetic algorithms). The latter can be quite heavy to
 implement in the present 8-dimension ($P$, $e$, $\o$, $T_0$, $K_1$,
 $K_2$, $\g_1$ and $\g_2$) parameter space. One of the advantages of
 the \citetalias{WHS67} algorithm is, indeed, the fact that no {\it a
 priori} assumption on the probable location of the exact solution in
 the parameter space is needed, except for the period (which can
 easily be obtained thanks to e.g.\ a Fourier analysis technique). It
 is thus efficient in terms of computer resources and it is robust for
 a large range of eccentricities (up to $e \approx 0.8$ according to
 WHS67). These characteristics make very attractive the WHS67
 method. In the present work, our aim has been to build an algorithm able to consistently deal with SB2 system while preserving the advantages of the \citetalias{WHS67} method.

To do so, we adopted the following idea as a guideline. The input SB2 data set is converted into an equivalent SB1 data set by applying an appropriate transformation on the measured radial velocities. The SB1 set so created is then used as an input for the \citetalias{WHS67} method. As a last step, the obtained equivalent SB1 solution  undergoes an inverse transformation that finally yields the desired SB2 orbital solution. The sketch of  Fig. \ref{fig: WHS2} summarizes the basic reasoning of the proposed method. The next section describes into more details the transformation adopted.  \\

\section{The transformation in the RV space}

The transformation from the SB2 to the equivalent SB1 data set relies on the following observation. 
Given the Eqs.~\ref{eq: rv1ps} and the fact that, in a Keplerian binary, $\o_2=\o_1+\pi$, one can write: 
\begin{equation}
 \frac{v_1(\t)-\gamma_1}{K_1}= -\frac{v_2(\t)-\gamma_2}{K_2} \label{eq: rv3}
\end{equation}
or equivalently 
\begin{equation} 
v_2(\t) =  b + c v_1(\t) \label{eq: rv6}
\end{equation}
with $b=\gamma_2-\frac{K_2}{K_1}\gamma_1$ and
$c=-\frac{K_2}{K_1}$. This latter equation is linear in the parameters
$b$ and $c$  and, given a set of $k=1...N$ observation couples
$\left(v_1(\t_k),v_2(\t_k)\right)$, the system can be solved using a
linear least-square technique. However, by opposition to the case
usually encountered, both variables $v_1$ and $v_2$ are spoiled by
errors of similar magnitudes. We therefore apply a so-called
orthogonal linear regression, the details of which are discussed in a
dedicated section (Sect.~\ref{sect: orth}). For the present time, let
us assume that, given an appropriate technique, the  $b$  and $c$
parameters and an estimate of their related errors, can be
determined. At this stage, we note that these parameters already
provide valuable information on the physical properties of the binary
system as $c=-\frac{K_2}{K_1}=-\frac{M_1}{M_2}$, so that the mass
ratio is directly obtained  without any {\it a priori} assumption on
the system parameters.\\

\begin{figure}
\centering
\includegraphics[width=15cm]{whs2}
\caption{Sketch of the underlying idea for adapting the \citetalias{WHS67} method to SB2 system.}
\label{fig: WHS2}
\end{figure}

At this stage, one could think that our aim has been reached. Indeed
according to equation Eq. \ref{eq: rv6}, it is a child game, knowing
the values of $b$ and $c$, to convert the secondary RVs into
equivalent primary velocities. Starting with an SB2 set of $N$ couples
$\left(v_1(\t_k),v_2(\t_k)\right)_{k=1}^{k=N}$, we result with an SB1
set containing $2N$ velocity points: the $N$ primary RV measurements
and $N$ points resulting from the conversion of the secondary RVs to
equivalent primary data points. As schematized in Fig. \ref{fig:
  WHS2}, this equivalent SB1 set can now be used as an input of the
\citetalias{WHS67} algorithm. This provides the orbital elements for,
in the considered case, the primary component: $e$, $\o$, $T_0$, $K_1$
and $\gamma_1$ as well as their respective errors. We then apply the
``inverse'' transformation to derive the remaining secondary
parameters $K_2$ and $\gamma_2$ from the expressions of  $b$ and
$c$. Error estimates follow from the error propagation theory. As an
illustrating case, Table \ref{tab: hd248} presents to orbital solution
of the \hd248 binary (using the data set of \citealt{SRG01}) that was
computed using the here described method (column: $v_2 \to v_1$). It
also illustrates one of the drawback of this way of converting the SB2
data into an equivalent SB1 set. \hd248 is a binary with two almost
identical components. The strengths of their spectral lines are thus
very similar and their RVs are measured with a similar
accuracy. However, as seen from Table \ref{tab: hd248}, the orbital
elements associated with the primary are much better constrained than
those associated with the secondary. This results from the additional
conversion applied to recover the secondary parameters once the
primary solution has been obtained. This unfortunately yields an additional propagation of the errors, leading larger estimates for the secondary parameter uncertainties. Proceeding in the opposite way (i.e. converting the primary velocities into equivalent secondary ones, computing the secondary orbital solution and finally deriving the remaining elements for the primary) yields the symmetric situation. The secondary parameters are now much better constrained compared to the primary ones (see Table \ref{tab: hd248}, column $v_1 \to v_2$). Note that this problem only concerns the error estimates. The values of the orbital parameters are of course not affected. \\

To circumvent this drawback, we propose now an alternative approach
which, in some way, arrange things so that both components ``meet
halfway''. Rather than transforming the RVs of one component  into
equivalent RVs of the other component, while the actual measurements
of the RVs associated to the later components remain unchanged, we now
propose to transform both components velocities, in order to create a
{\it fake} SB1 data set that is equivalent to the SB2 system. The
orbital solution of this equivalent SB1 system can be obtained through
a modified version of the \citetalias{WHS67}
algorithm\footnote{Actually, the \citetalias{WHS67} algorithma had to
  be slightly updated because of the sign difference between the
  primary and secondary RV equations. This results from the fact that $\o_2=\o_1+\pi$ (see e.g. Eqs.~\ref{eq: rv1ps}).} and allows,  in a last step, to deduce the true SB2 solution.
This equivalent SB1 system will be noted by an asterisk ($^{*}$) in the following and is described by the equation :
\begin{equation}
v^*(\t)= \Gamma + \mathcal{K} \sin i (\cos(\o+\t)+e \cos\t). \label{eq: rv7}
\end{equation}
It has the same eccentricity $e$, the same longitude of periastron $\o$ and  the same time of periastron passage $T_0$ than the real SB2 system. In Eq.~\ref{eq: rv7}, $\Gamma$ and $\mathcal{K}$ are the systemic velocity and the semi-amplitude of the RV curve of the fake SB1 system. They  respectively correspond to the difference between the systemic velocities of the two components of the real SB2 system and to the geometric average of their RV curve semi-amplitudes :
\begin{eqnarray}
\Gamma     &=&\g_1-\g_2 \label{eq: rv8}\\
\mathcal{K}&=&\sqrt{K_1 K_2}\label{eq: rv9}
\end{eqnarray}

The corresponding velocity transformation is given by :
\begin{equation} \left\{ \begin{array}{lll}
v_1^*&=&\sqrt{-c} \left( v_1-\frac{b}{1-c}\right)  \vspace*{1mm} \\
v_2^*&=&\frac{1}{\sqrt{-c}} \left( v_2-\frac{b}{1-c} \right) \end{array} \right. \label{eq: SB1*}
\end{equation}
The two components are now symmetrically handled. When computing the primary and secondary RV curve semi-amplitudes and their related errors from Eq. \ref{eq: rv9}, one can show that : 
$\frac{\sigma_{K_1}}{K_1} = \frac{\sigma_{K_2}}{K_2}\label{eq: rv12}$.
The relative errors on the primary and secondary RV-curve semi-amplitudes are now equal and the previous asymmetric situation is avoided. Physically speaking, for a system with two identical components such as \hd248, the parameters specific to the individual components will now be determined with the same accuracy. In the case of two different stars in a binary system, their respective RV-curve semi-amplitude will be measured with a relative accuracy related to their mass-ratio: $\frac{\sigma_{K_2}}{\sigma_{K_1}}=\frac{K_2}{K_1}=\frac{M_1}{M_2}$. Though it can not be appropriate for all binary systems, it has been successfully apply to various SB2 massive binaries. %Indeed, for SB2 system, we observe the spectral signature of the two binary components. With the current instrumental facilities used, this implies that their brightness should not differ by more than about one order of magnitude.  The two component of a SB2 massive binary are thus most of the time of the same luminosity class. The mass-luminosity relation thus indicates that, compared to its companion,  the more massive the component is, the brightest it is, the more intense its spectral lines are and thus the more accurately they can be measured. At the first order this is well expressed by Eq. \ref{eq: rv12}. Other effects, such as rotation, should not much disturb this sketch much as, in massive binaries, the star rotation is usually synchronized with the orbital rotation rate, so that none of the component should present lines very much broader (and thus of which the centroid is more difficult to find) than the other.

%*********************************************************
\begin{table}
\centering
\caption{Orbital and physical parameters of the \hd248 binary computed
  using various techniques. The data were taken from
  \citet{SRG01}. The same weights have been given to the primary and
  secondary points.  The notations of Table \ref{tab: orbparam} have
  been used. The period has been kept constant to render a direct comparison more easy. $x \to y $ has the meaning of `conversion from $x$ to $y$'. }
\label{tab: hd248}
\label{tab: all_orb}
\begin{tabular}{c | c c | c c c}
\hline
\hline
                      & \multicolumn{2}{c}{Individual solutions} & \multicolumn{3}{c}{Common solutions}\\
                      & $v_1$             &         $ v_2$     & $v_2 \to v_1$         & $v_1 \to  v_2$       &  $\left. \begin{array}{l} v_1 \\ v_2 \end{array} \right\} \to v^*$ \\
\hline		       	   	       	       	   	       	    	   	       	       	   	       	       	   	       
$P$ (d)               & 5.816005          & 5.816005           &   5.816005            &   5.816005            &   5.816005            \\
$s_y/s_x$             & n.a.              & n.a.               &   1.0                 &   1.0                 &   1.0                 \\
$m_1/m_2$             & n.a.              & n.a.               &   0.989 $\pm$ 0.011   &   0.989 $\pm$ 0.011   &   0.989 $\pm$ 0.011   \\
$e      $             & 0.133 $\pm$ 0.006 & 0.134 $\pm$ 0.007  &   0.134 $\pm$ 0.006   &   0.134 $\pm$ 0.006   &   0.134 $\pm$ 0.005   \\
$\omega_1 $ (\degr)   &   80.6 $\pm$ 4.2  &  85.622 $\pm$ 4.8  &    83.2 $\pm$ 4.2     &    83.2 $\pm$ 4.2     &    83.1 $\pm$ 3.2     \\
$T_0$ (HJD            &  2003.827         &  2003.863          &    2003.846           &    2003.846           &    2003.846           \\
$-$2\,450\,000)       & $\pm$ 0.064       & $\pm$ 0.073        &   $\pm$ 0.064 	       &   $\pm$ 0.064         &   $\pm$ 0.050 	      \\
$K_1$\ (\kms)         & 216.6 $\pm$ 1.2   & n.a.               &   216.7 $\pm$ 1.3     &   216.7 $\pm$ 4.2     &   216.7 $\pm$ 1.5     \\
$K_2$\ (\kms)         & n.a.              & 214.3 $\pm$ 1.5    &   214.2 $\pm$ 4.2     &   214.2 $\pm$ 1.3     &   214.2 $\pm$ 1.5     \\
$\gamma_1$\ (\kms)    &  -0.3 $\pm$ 1.4   & n.a.               &     0.7 $\pm$ 1.4     &     0.7 $\pm$ 3.4     &     0.7 $\pm$ 1.5     \\
$\gamma_2$\ (\kms)    & n.a.              &  -0.6 $\pm$ 1.6    &     0.4 $\pm$ 3.4     &     0.4 $\pm$ 1.4     &     0.4 $\pm$ 1.5     \\
r.m.s. (\kms)         & 7.0               & 8.2                &   8.0                 &   8.0                 &   8.0                 \\
\hline				  
\end{tabular}
\end{table}
%*********************************************************


\section{The orthogonal linear regression } \label{sect: orth}

This section summarizes the main properties of the orthogonal linear regression technique applied to fit the $v_2$ vs. $v_1$ relation (Eq.~\ref{eq: rv6}). The detailed calculation are provided in App.~\ref{app: orth}.
Let us consider a data set formed by $N$ couples of measurements $(x_k,y_k)$ to which one want to fit the linear model  $y=cx+b$. The usual linear least-square technique minimizes the merit function :
\begin{equation}
\chi^2=\sum_{k=1}^N \Delta_k^2, \hspace{0.5cm}\mathrm{ with } \hspace{0.5cm} \Delta_k=y_k-cx_k-b, \label{eq: LS} \end{equation}
eventually taking into account various weights or error estimates for the measurements. Eq.~\ref{eq: LS} actually minimizes the sum of the square deviations ($\Delta_k$) to the best fit model, these deviations being measured along the $y$-axis. This is an intuitive approach when $\sigma_y>>\sigma_x$. Alternatively, when $\sigma_y<<\sigma_x$, one could use $\Delta_k=x_k-\frac{y_k+b}{c}$ (as long as  $c\ne0$). However, when the two variables have uncertainties of the same order of magnitude, one may rather prefer to measure the $\Delta_k$'s along a direction which is orthogonal to the best fit model rather than along one given axis or along the other. This approach falls under the generic terms of ``orthogonal linear regression''. When all measurements have the same uncertainties, the latter principle yields $\Delta_k^2=\frac{(y_k-cx_k-b)^2}{1+c^2}$. When dealing with various uncertainties, one can write $\Delta_k^2=\frac{(y_k-cx_k-b)^2}{\sigma_{y_k}^2+c^2\sigma_{x_k}^2}$. When the ratios of the secondary to primary uncertainties are all equal whatever the considered measurement couple $k$, this can be rewritten as: $\Delta_k^2=\frac{1}{\sigma_{x}^2} \frac{(y_k-cx_k-b)^2}{s^2+c^2}$ where $s=\sigma_{y}/\sigma_{x}=\sigma_{y_k}/\sigma_{x_k}$\ $\forall k$. \\

The latter situation is often encountered while dealing with SB2 RV measurements. Indeed one can often estimate that, given the considered stars, the primary RVs are $s$ times more accurate that the secondary ones. However, it is often more difficult to obtain an accurate estimate of the amplitude of these uncertainties. Experimental scientists have since long been facing such kind of problems, which they usually circumvent by adopting a proper relative weighting for the measurements rather than an absolute weighting. In this latter case, the appropriate merit function can be written :
\begin{equation}
\chi^2=\sum_{k=1}^N w_k \frac{(y_k-cx_k-b)^2}{s^2+c^2}, \label{eq: OLS} 
\end{equation}
where $w_k$ is the weight given to the measurements of the couple
$(x_k,y_k)$. The present approach is equivalent to the one that adopts
the primary weight $w_{p_k}=w_k$ while the secondary ones are given by $w_{s_k}=w_{p_k}/s^2$. In this configuration, one can show that the orthogonal least square solution is given by (App.~\ref{app: orth}) :
\begin{eqnarray}
a &=& s (F+\sqrt{F^2+1}), \hspace{0.5cm} \mathrm{ with } \hspace{0.5cm} F=\frac{S_y-s^2 S_x}{2sC_{xy}} \label{eq: a}\\
b &=& \y - c \x. \label{eq: b}
\end{eqnarray} 

In the latter equations, one have used the following notations : \begin{equation} \begin{array}{lll}
S_x = \xb-\x^2     & \hspace*{1.0cm}\x = \frac{ \sum_{k=1}^N w_k x_k}{\sum_{k=1}^N w_k } & \hspace*{1.0cm} \xb=\frac{ \sum_{k=1}^N w_k x_k^2}{\sum_{k=1}^N w_k } \vspace*{1mm}\\
S_y = \yb-\y^2     & \hspace*{1.0cm} \y = \frac{ \sum_{k=1}^N w_k y_k}{\sum_{k=1}^N w_k } & \hspace*{1.0cm} \yb=\frac{ \sum_{k=1}^N w_k y_k^2}{\sum_{k=1}^N w_k }\vspace*{1mm}\\
C_{xy}=\xy-\x\ \y    & \hspace*{1.0cm} \xy = \frac{ \sum_{k=1}^N w_k x_k y_k}{\sum_{k=1}^N w_k }\\
\end{array} \end{equation}


Under these hypotheses and notations, it is possible to show that the theory of error propagation yields an analytical expressions for the related uncertainties :
\begin{eqnarray}
\sigma_c^2 & = & \frac{\chi^2_\nu}{N} c^2 \mathcal{F}  \\
\sigma_b^2 & = & \frac{\chi^2_\nu}{N} (s^2+c^2+  c^2 \x^2 \mathcal{F} )
\end{eqnarray}
where 
\begin{equation} 
\mathcal{F}=\frac{S_y+s^2 S_x}{C_{xy}^2} \hspace*{0.5cm}\mathrm{ and }\hspace*{0.5cm} \chi^2_\nu=\frac{N \chi^2}{(N-2) \sum_{k=1}^N w_k}. \end{equation}


\section{Final remark}

The present method requires a preliminary guess of the orbital method
as well as an estimate of the ratio  $s$  between the secondary and
primary uncertainties. While the first can be obtained through e.g. a
Fourier analysis, the latter might be more difficult to accurately
estimate. The present method can provide an original way to circumvent
this problem. Indeed the method is based on the construction of a fake
SB1 data set that is equivalent to the SB2 physical system. The
construction of this fake set relies on the best fit parameters (see
Eqs. \ref{eq: a} and \ref{eq: b}) of the $v_2$ vs. $v_1$ relation and given. The latter solutions depend to some extent on the adopted value of $s$. As a consequence, if the adopted $s$ value does not reflect the true ratio between the primary and secondary uncertainty, the conversion of the primary and secondary RVs will provide inconsistent contributions to the {\it fake} SB1 system. Thus, the final fit will not be good and large differences between the data and the final SB2 solution  will be observed. Hence, one can naturally use an exploration of the parameter space to estimate the $s$ ratio while minimizing the square deviations between the observations and the computed SB2 solutions:
\begin{equation}
\chi^2_\mathrm{SB2}=\frac{\sum_{k=1}^{k=N}w_{p_k} \left( v_{p_k}-\gamma_1-K_1 \left( \cos\left(\o+\t_k\right) + e \cos\o \right)\right)^2 + \sum_{k=1}^{k=N}w_{s_k} \left( v_{s_k}-\gamma_2+K_2 \left( \cos\left(\o+\t_k\right) + e \cos\o \right)\right)^2}{\sum_{k=1}^{k=N}w_{p_k}+\sum_{k=1}^{k=N}w_{s_k}}.
\end{equation}
 As a conclusion and assuming that $s$ is also a free parameter of the
 model, the method  presented here allows to restrict the exploration
 of the 9-dimension (actually 8 dimensions + one for $s$) parameter space to a single direction corresponding to $s$.


\section{Summary}
We have shown that, using an appropriate transformation, it is possible to adapt the \citetalias{WHS67} algorithm to derive consistent orbital solutions for the two parameters of an SB2 system. The propose method relies on a linear orthogonal regression technique that yields the best fit parameters of the relation $v_2=cv_1+b$, where $c=-K_2/K_1=-M_1/M_2$ and $b=\gamma_2+c\gamma_1$. Hence the mass ratio is obtained with no {\it a priori} assumption on the system, not even a first guess of the orbital period. We have shown that, if the ratio of the primary and secondary RV measurements can be taken as constant for all the measurements, it is possible to derive a simple analytical expression to estimate the uncertainties on the best fit parameters $b$ and $c$ and, hence, on the system mass-ratio. We outline that the present algorithm presents much  of the advantages of the original \citetalias{WHS67} method. In particular, it only requires a first guess of the orbital period. All the other orbital parameters are obtained without any {\it a priori} assumption on their location in the parameter space. Hence, the method does not need to explore the 8-dimension parameters and, as such, is pretty efficient in terms of computation time. Finally, as the \citetalias{WHS67} method, the present algorithm is adapted to a wide range of eccentricities although a detailed study of the break up point is still lacking.

%_____________ ACKNOWLEDGMENTS __________________________________________


\section*{Acknowledgments}
The original idea of transforming the RVs of one component into equivalent RVs for the other component was initially suggested by Dr. G. Rauw.

%_____________ BIBLIOGRAPHY _____________________________________________

\bibliographystyle{mn2e}
\bibliography{/datas6/XMM_CAT_PAPER/ngc6231_Xcat}



\appendix
\section{Orthogonal regression: detailed calculations} \label{app: orth}
Let the merit function given by Eq.~\ref{eq: OLS} : 
\begin{equation}
 \chi^2=\sum_{k=1}^N w_k\frac{\left( y_k-cx_k-b\right) ^2}{s^2+c^2}
\end{equation}
  be minimized with respect to the parameters of the model $c$ and $b$. We thus set the derivative along $b$ and $c$ equal to zero:

\begin{eqnarray}
\frac{\partial \chi^2}{\partial b} & = & -\frac{2}{s^2+c^2}\sum_k w_k\left( y_k-cx_k-b\right)  \iff b   =   \y - c \x 
\end{eqnarray}

and

\begin{eqnarray}
 \frac{\partial \chi^2}{\partial c}  &=  &-\frac{2c}{\left( s^2+c^2\right) ^2} \sum_k w_k \left( y_k-cx_k-b\right) ^2 - \frac{2}{s^2+c^2}\sum_k w_k x_k \left( y_k-cx_k-b\right)   = 0  \vspace*{1mm}\\
\iff && c\sum_k w_k\left( y_k^2+c^2x_k^2+b^2-2cx_ky_k-2by_k+2cbx_k\right) +\left( s^2+c^2\right) \sum_k w_k \left( x_ky_k-cx_k^2-bx_k\right)  = 0 \vspace*{1mm}\\
\iff &&  c \yb + \underline{c^3 \xb} + cb^2-\underline{2}c^2\xy-2cb\y+\underline{2}c^2b\x + \left( s^2+\underline{c^2}\right) \xy - \left( s^2+\underline{c^2}\right) c\xb - \left( s^2+\underline{c^2}\right) b\x  = 0\vspace*{1mm}\vspace*{1mm}\\
\iff && c\yb + c \left( \underline{\y^2}-\underline{2c\x\ \y}+\underline{c^2\x^2}\right) -c^2\xy-\underline{2}c\y^2+\underline{2c^2\x\ \y}+c^2\x\ \y-\underline{c^3\x^2}+s^2\xy-cs^2\xb-s^2\x\ \y+cs^2\x^2 = 0\\
\iff && c^2\left( \x\ \y-\xy\right) +c\left( \yb-\y^2-s^2\xb+s^2\x^2\right) +\left( s^2\xy-s^2\x\ \y\right)  =  0 \vspace*{1mm}\\
\iff && -c^2C_{xy}+c\left( S_y-s^2S_x\right) +s^2C_{xy} = 0 \vspace*{1mm}\\
\iff && c  =  \frac{\left( S\y-s^2S_x\right) \pm\sqrt{\left( S_y-s^2S_x\right) ^2+4C_{xy}s^2}}{2C_{xy}}\vspace*{1mm}\\
\iff && c   =  s \left( F\pm\sqrt{F^2+1}\right)  
\end{eqnarray}

with $F = \frac{S_y-s^2S_x}{2sC_{xy}}$ as given by Eq.~\ref{eq: a}. One can show that the $+$ solution is indeed associated to a minimum of the $\chi^2$ function. Before computing the errors, let us first estimate the partial derivatives of $c$:

\begin{eqnarray}
\frac{\partial c}{\partial F} & = & s \left( 1+\frac{1}{2} \frac{2F}{\sqrt{F^2+1}}\right)  =  s\frac{\sqrt{F^2+1}+F}{\sqrt{F^2+1}} = \frac{c}{\sqrt{F^2+1}}\vspace*{1mm}\\
F & = & \frac{S_y-s^2S_x}{2sC_{xy}} = \frac{\yb-\y^2-s^2\left( \xb-\x^2\right) }{2s\left( \xy-\x\ \y\right) }\vspace*{1mm}\\
\frac{\partial F}{\partial y_k} & = & \frac{\left( \frac{2w_k y_k}{\sum_l w_l}-\frac{2w_k\y}{\sum_l w_l}\right) 2sC_{xy} - \left( S_y-s^2S_x\right) 2s\left( \frac{w_kx_k}{\sum_l w_l}-\frac{w_k\x}{\sum_l w_l}\right)}{4s^2C^2_{xy}}\vspace*{1mm}\\
 & = &  \frac{w_k}{\sum_l w_l}\frac{1}{2sC^2_{xy}}\left( \left( y_k-\y\right) 2C_{xy}-\left( S_y-s^2S_x\right) \left( x_k-\x\right) \right)\vspace*{1mm}\\
 & = &  \frac{w_k}{\sum_l w_l}\frac{1}{sC_{xy}}\left( \left( y_k-\y\right) -\frac{S_y-s^2S_x}{2C_{xy}}\left( x_k-\x\right) \right)\vspace*{1mm}\\
 & = &  \frac{w_k}{\sum_l w_l}\frac{1}{C_{xy}}\left( \frac{1}{s}\left( y_k-\y\right) -F\left( x_k-\x\right) \right)\vspace*{1mm}\\
\frac{\partial F}{\partial x_k} & = & \frac{-s^2\left( \frac{2w_k x_k}{\sum_l w_l}-\frac{2w_k\x}{\sum_l w_l}\right) 2sC_{xy} - \left( S_y-s^2S_x\right) 2s\left( \frac{w_ky_k}{\sum_l w_l}-\frac{w_k\y}{\sum_l w_l}\right) }{4s^2C^2_{xy}}\vspace*{1mm}\\
 & = &  \frac{w_k}{\sum_l w_l}\frac{1}{C_{xy}}\left( -s\left( x_k-\x\right) -\frac{S_y-s^2S_x}{2sC_{xy}}\left( y_k-\y\right) \right)\vspace*{1mm}\\
 & = &  \frac{w_k}{\sum_l w_l}\frac{-1}{C_{xy}}\left( s\left( x_k-\x\right) +F\left( y_k-\y\right) \right)
\end{eqnarray}

Now let us compute the uncertainty on $c$ by mean of the error propagation theory \citep[see e.g.][]{Bev69}. In the following we note the absolute uncertainty $\s_{x_k}=\s / w_k$ and we assume that the $w_k$ have been normalized. Similarly, we obtained  $\s_{y_k}=s^2 \s / w_k$. Hence :

\begin{eqnarray}
\sb_c &=& \sum_k \left( \left( \frac{\partial c}{\partial y_k}\right)^2 \sb_{y_k} + \left( \frac{\partial c}{\partial x_k}\right) ^2 \sb_{x_k} \right)\vspace*{1mm}
\\
 &=& \sb \sum_k \frac{1}{w_k}  \left( \frac{\partial c}{\partial F}\right)^2 \left(  \left(\frac{\partial F}{\partial x_k}\right) ^2 +  \left(\frac{\partial F}{\partial y_k}\right)^2 s^2 \right) \vspace*{1mm}
\\
&=& \frac{\sb c^2}{1+F^2} \frac{1}{C^2_{xy}} \frac{1}{\left( \sum_l w_l\right)^2} \sum_k w_k \Bigg( s^2 \left( x_k-\x\right) ^2+F^2\left( y_k-\y\right) ^2+ \underline{2sF\left( y_k-\y\right) \left( x_k-\x\right) } \vspace*{1mm}\\ 
&& \hspace*{6cm} +  s^2\left( \frac{1}{s^2}\left( y_k-\y\right) ^2+F^2\left( x_k-\x\right)^2- \underline{ \frac{2}{s}F\left( y_k-\y\right) \left( x_k-\x\right) } \right) \Bigg)\vspace*{1mm}
 \\ 
&=& \frac{\sb c^2}{1+F^2}\frac{1}{C^2_{xy}} \frac{1}{\left( \sum_l w_l\right)^2} \sum_k w_k \left( \left( x_k-\x\right) ^2s^2\left( 1+F^2\right) +\left( y_k-\y\right) ^2\left( F^2+1\right) \right) \vspace*{1mm}
\\
&=& \frac{\sb c^2}{C^2_{xy}} \frac{1}{\left( \sum_l w_l\right) ^2} \left( \sum_k  w_k s^2 \left( x_k-\x\right) ^2 +  \sum_k  w_k  \left( x_y-\y\right) ^2 \right) \vspace*{1mm}
\\
&=& \frac{\sb c^2}{\sum_l w_l}\frac{s^2S_x+S_y}{C^2_{xy}}
\end{eqnarray}

 and for $b$ :
\begin{eqnarray}
\frac{\partial b}{\partial y_k} & = &\frac{w_k}{\sum_l w_l} -\x \frac{\partial c}{\partial y_k}\vspace*{1mm}
 \\
 & = &\frac{w_k}{\sum_l w_l} -\x \frac{c}{\sqrt{F^2+1}}\frac{\partial F}{\partial y_k} \vspace*{1mm}
\\
 & = &\frac{w_k}{\sum_l w_l} -\x \frac{c}{\sqrt{F^2+1}}\frac{w_k}{\sum_l w_l}\frac{1}{C_{xy}}\left( \frac{1}{s}\left( y_k-\y\right)-F\left( x_k-\x\right)\right) \vspace*{1mm}
\\
 & = &\frac{w_k}{\sum_l w_l} \left( 1- \frac{c\x}{\sqrt{F^2+1}}\frac{1}{C_{xy}}\left( \frac{1}{s}\left( y_k-\y\right)-F\left( x_k-\x\right)\right)\right)
\vspace*{1mm} \\
\frac{\partial b}{\partial x_k} & = &- \frac{w_k}{\sum_l w_l} c -\x \frac{\partial c}{\partial x_k}\vspace*{1mm}
 \\
 & = &-\frac{w_k}{\sum_l w_l} \left( c- \frac{c\x}{\sqrt{F^2+1}}\frac{1}{C_{xy}}\left( s\left( x_k-\x\right)-F\left( y_k-\y\right)\right)\right) 
\end{eqnarray}

\begin{eqnarray}
\sb_b & = & \sb \sum_k\frac{1}{w_k}\left( \left( \frac{\partial b}{\partial x_k} \right)^2 + \left( \frac{\partial b}{\partial y_k} \right)^2 s^2\right) 
\vspace*{1mm}
\\
&=& \frac{\sb}{\left(\sum_l w_l\right)^2} \sum_k w_k \Bigg( c^2+\frac{c^2\x^2}{C^2_{xy}}\frac{1}{F^2+1}\left(s\left(x_k-\x\right)+F\left(y_k-\y\right)\right)^2 + s^2+\frac{c^2\x^2}{C^2_{xy}}\frac{s^2}{F^2+1}\left(\frac{1}{s}\left(y_k-\y\right)-F\left(x_k-\x\right)\right)^2 \Bigg)
\vspace*{1mm}
\\
& = & \frac{\sb}{\left(\sum_l w_l\right)^2} \left( \sum_k w_k \left(c^2+s^2\right)+\sum_k w_k \frac{c^2\x^2}{C^2_{xy}}\frac{1}{F^2+1}\left( s^2 \left( x_k-\x \right)^2+F^2\left( y_k-\y \right)^2 +\left( y_k-\y \right)^2 +s^2 F \left( x_k-\x \right)^2\right) \right)
\vspace*{1mm} 
\\
& = & \frac{\sb}{\sum_l w_l} \left( c^2+s^2 + \frac{c^2\x^2}{C^2_{xy}}\frac{1}{\sum_l w_l} \sum_k w_k \left( s^2 \left( x_k-\x \right)^2 + \left( y_k-\y \right)^2 \right) \right)
\vspace*{1mm}
\\
& = & \frac{\sb}{\sum_l w_l} \left( c^2+s^2 + c^2\x^2\frac{s^2S_x+S_y}{C^2_{xy}} \right)
 \end{eqnarray}

Finally, one obtains the absolute magnitude of the errors by estimating $\sigma$. To do so, we assume that the deviation between the data and the model is only due to the measurements uncertainties; then :
\begin{equation}
\sb \approx \frac{\sum_k w_k \left(y_k-cx_k-b\right)^2}{\left(N-2\right) \left(s^2+c^2\right)}=\frac{\chi^2}{N-2}
\end{equation}
where, as previously, the $w_k$ are assumed to be normalized.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The Li\`ege SB2 Orbital Solution Algorithm: an exhaustive derivation}

This section presents the details of the Li\`ege   algorithm used to compute SB2 orbital solution. Let us first write Eq.~\ref{eq: rv2} for each star :
\begin{equation} 
\left\{ \begin{array}{ll}
v_1\left(\t\right)=\gamma_1 + \sin i \sin \o_1 \frac{d}{dt}\left(r_1 \cos \t\right) + \sin i \cos \o_1 \frac{d}{dt}\left(r_1 \sin \t\right) \vspace*{1mm}\\
v_2\left(\t\right)=\gamma_2 + \sin i \sin \o_2 \frac{d}{dt}\left(r_2 \cos \t\right) + \sin i \cos \o_2 \frac{d}{dt}\left(r_2 \sin \t\right)  \end{array} \right. \label{eq: v1v2}\end{equation}

Let us note $\o_1=\o$ and $\o_2=\o+\pi$ (hence $\sin \o_2=-\sin \o$ and  $\cos \o_2=-\cos \o$) and let us adopt the notation $M$ for the mean anomaly and $E$ for the eccentric anomaly. It is well known that : \begin{equation} M=E-e\sin E=\frac{2\pi}{P} \left(t-T_0\right)=\mu\left(t-T_0\right)\end{equation}

\subsection{The harmonic analysis}
Following the theory of elliptic motion, let us now rewrite the terms $r \cos \t$ and $r \sin \t$ using harmonic series of the mean anomaly \citep[see e.g.][p.~73]{BC61}:

\begin{eqnarray}
\frac{r}{a}\cos \t & = & \cos E-e \\ 
                   & = & \frac{-3e}{2}-2\sum_{s=1}^{\infty}\frac{1}{s}J'_s\left(se\right) \cos\left(sM\right) \label{eq: BesselC}\\
\frac{r}{a}\sin \t & = & \sqrt{1-e^2}  \sin E \\
                   & = & \sqrt{1-e^2} \left(\frac{2}{e} \sum_{s=1}^{\infty}\frac{1}{s}J_s\left(se\right) \sin\left(sM\right)\right) \label{eq: BesselS}
\end{eqnarray}
where the $J_s(x)$ are the Bessel functions of the first kind and $J'_s(x)$, their first derivatives. In a next step, one develops the Bessel functions in asymptotic series of the eccentricity. Following \citet[pp.~79-80]{BC61}, one can write :

\begin{equation} \begin{array}{lllll}
\frac{r}{a}\cos \t & \approx & -\frac{3}{2}e+\left(1-\frac{3}{8}e^2+\frac{5}{192}e^4-\frac{7}{9216}e^6\right) \cos M  
                               +\left(\frac{1}{2}e-\frac{1}{3}e^3+\frac{1}{16}e^5-\frac{1}{180}e^7\right) \cos2M
                               +\left(\frac{3}{8}e^2-\frac{45}{128}e^4+\frac{567}{5120}e^6\right) \cos3M  \vspace*{1mm}\\    
                          & &  +\left(\frac{1}{3}e^3-\frac{2}{5}e^5+\frac{8}{45}e^7\right) \cos4M  
                               +\left(\frac{125}{384}e^4-\frac{4375}{9216}e^6\right) \cos5M  
                               +\left(\frac{27}{80}e^5-\frac{81}{140}e^7\right) \cos6M  
                               +\frac{16807}{46080}e^6 \cos7M     \vspace*{1mm}      \\     
& = & -\frac{3}{2}e 
      + X_1 \cos M  
      +\frac{e}{2} X_2 \cos 2M\vspace*{1mm} 
      +\frac{e^2}{3} X_3 \cos 3M  
      +\frac{e^3}{4} X_4 \cos 4M 
      +\frac{e^4}{5} X_5 \cos 5M 
      +\frac{e^5}{6} X_6 \cos 6M 
      +\frac{e^6}{7} X_7 \cos 7M   
\end{array}\end{equation}

\begin{equation} \begin{array}{lllll}
\frac{r}{a}\sin \t & \approx & \left(1-\frac{5}{8}e^2+\frac{11}{192}e^4-\frac{457}{9216}e^6\right) \sin M  
                              +\left(\frac{1}{2}e-\frac{5}{12}e^3+\frac{1}{24}e^5-\frac{1}{45}e^7\right) \sin2M  
                              +\left(\frac{3}{8}e^2-\frac{51}{128}e^4+\frac{543}{5120}e^6\right) \sin3M  \vspace*{1mm}\\  
                   &  &       +\left(\frac{1}{3}e^3-\frac{13}{30}e^5+\frac{13}{72}e^7\right) \sin4M 	           
                              +\left(\frac{125}{384}e^4-\frac{4625}{9216}e^6\right) \sin5M 		           
                              +\left(\frac{27}{80}e^5-\frac{135}{224}e^7\right) \sin6M 		   	   
                              +\frac{16807}{46080}e^6 \sin7M			 \vspace*{1mm}\\ 
& = &   Y_1 \sin M    		
      +\frac{e}{2} Y_2 \sin 2M  	
      +\frac{e^2}{3} Y_3 \sin 3M     
      +\frac{e^3}{4} Y_4 \sin 4M     
      +\frac{e^4}{5} Y_5 \sin 5M      
      +\frac{e^5}{6} Y_6 \sin 6M
      +\frac{e^6}{7} Y_7 \sin 7M     
\end{array}\end{equation}

where the $X_s$ and $Y_s$  are defined so that the equation is verified and depend only on $e$. In the present development, the harmonic analysis given by Eqs. \ref{eq: BesselC} and \ref{eq: BesselS} have been limited to $s=7$, which corresponds to neglect powers of $e$ above the seventh. Both developments can be rewritten in short by summing on $k$ :
\begin{eqnarray}
r \cos \t &=& -\frac{3}{2}e +\sum_{k=1}^7 a \frac{e^{k-1} X_k}{k}\cos kM \\  
r \sin \t &=& \sum_{k=1}^7 a \frac{e^{k-1} Y_k}{k}\sin kM 
\end{eqnarray}
The derivatives with respect the time $t$ can now be easily evaluated :
\begin{eqnarray}
\frac{d}{dt}\left(r \cos \t\right) &=& -a\mu\sum_{k=1}^7 e^{k-1} X_k \sin kM \\  
\frac{d}{dt}\left(r \sin \t\right) &=&  a\mu\sum_{k=1}^7 e^{k-1} Y_k \cos kM 
\end{eqnarray}
so that the system \ref{eq: v1v2} becomes :\begin{equation}
\left\{ \begin{array}{l}
v_1=\gamma_1 + a_1 \mu \sin i \left( - \sin \o\sum_{k=1}^7 e^{k-1} X_k \sin kM+ \cos\o \sum_{k=1}^7 e^{k-1} Y_k \cos kM \right) \vspace*{1mm}\\
v_2=\gamma_2 - a_2 \mu \sin i \left( - \sin \o\sum_{k=1}^7 e^k X_k \sin kM +\cos\o \sum_{k=1}^7 e^k Y_k \cos kM \right) \end{array} \right. \label{eq: step2} \end{equation}

However the time of periastron passage $T_0$ is not known so that one can not compute $M$ directly from the observations. Adopting an arbitrary time $t_0$, one can write 
\begin{equation}
M=\mu \left(t-T_0\right)=\mu \left(t-t_0\right)+\mu \left(t_0-T_0\right)=M'+M_0.
\end{equation}
Since 
 \begin{eqnarray}
\cos k\left(M'+M_0\right)=\cos kM' \cos kM_0 - \sin kM' \sin kM_0 \\
\sin k\left(M'+M_0\right)=\sin kM' \cos kM_0 + \cos kM' \sin kM_0 
\end{eqnarray}
the system \ref{eq: step2} becomes : 

\begin{equation} \left\{ \begin{array}{lll}
v_1&=&\gamma_1 + a_1 \mu \sin i \bigg( - \sin \o \sum_{k=1}^7 e^{k-1} X_k \sin kM' \cos kM_0  - \sin \o \sum_{k=1}^7 e^{k-1} X_k \cos kM' \sin kM_0\\
   & & \hspace{2.2cm}  + \cos\o \sum_{k=1}^7 e^{k-1} Y_k \cos kM' \cos kM_0  - \cos\o \sum_{k=1}^7 e^{k-1} Y_k \sin kM' \sin kM_0 \bigg) \\
v_2&=&\gamma_2 - a_2 \mu \sin i \bigg( - \sin \o \sum_{k=1}^7 e^{k-1} X_k \sin kM' \cos kM_0  - \sin \o \sum_{k=1}^7 e^{k-1} X_k \cos kM' \sin kM_0\\
   & & \hspace{2.2cm}  + \cos\o \sum_{k=1}^7 e^{k-1} Y_k \cos kM' \cos kM_0  - \cos\o \sum_{k=1}^7 e^{k-1} Y_k \sin kM' \sin kM_0 \bigg)
 \end{array} \right. \end{equation}

Grouping the coefficients of the unknowns $\sin kM'$ and $\cos kM'$, one writes :

\begin{equation} \left\{ \begin{array}{lll}
v_1&=&\gamma_1 +  \sum_{k=1}^7 \sin kM' e^{k-1} a_1 \mu \sin i \left( - \sin \o\  X_k  \cos kM_0 - \cos\o \sin kM_0 \right) \vspace*{1mm}\\
   & &  \hspace{2.2cm}   + \sum_{k=1}^7 \cos kM' e^{k-1} a_1 \mu \sin i \left( - \sin \o X_k  \sin kM_0  + \cos\o  \cos kM_0 \right) \vspace*{1mm}\\
v_2&=&\gamma_2 - \sum_{k=1}^7 \sin kM' e^k a_2 \mu \sin i \left( - \sin \o\  X_k  \cos kM_0 - \cos\o \sin kM_0 \right) \vspace*{1mm}\\
   & &  \hspace{2.2cm}   - \sum_{k=1}^7 \cos kM' e^k a_2 \mu \sin i \left( - \sin \o X_k  \sin kM_0  + \cos\o  \cos kM_0 \right)\\ \end{array} \right. \label{eq: step3} \end{equation}

Adopting 
\begin{eqnarray} 
S_k^{(j)}&=&-\mu e^{k-1}a_j \sin i \left( \sin \o X_k \cos kM_0 + \cos \o Y_k \sin kM_0 \right) \label{eq: Sk}\\
C_k^{(j)}&=&+\mu e^{k-1}a_j \sin i \left( - \sin \o X_k \sin kM_0 + \cos \o Y_k \cos kM_0 \right)\label{eq: Ck}
\end{eqnarray}
 
the system \ref{eq: step3} takes a much shortened appearance:

\begin{equation} \left\{ \begin{array}{lll}
v_1&=&\gamma_1 + \sum_{k=1}^7 S_k^{(1)} \sin kM' +\sum_{k=1}^7 C_k^{(1)} \cos kM' \vspace*{1mm}\\
v_2&=&\gamma_2 - \sum_{k=1}^7 S_k^{(2)} \sin kM' +\sum_{k=1}^7 C_k^{(2)} \cos kM' \\ \end{array} \right. \label{eq: step4} \end{equation}

Before transforming the SB2 system into a fake but equivalent SB1 system, it is useful to note that
\begin{equation}
\frac{S_k^{(1)}}{a_1}=\frac{S_k^{(2)}}{a_2}\hspace{0.5cm}\mathrm{and}\hspace{0.5cm}\frac{C_k^{(1)}}{a_1}=\frac{C_k^{(2)}}{a_2} \\
\end{equation}
Let us note by an asterisk ($^*$) the parameters relative to the fake SB1 system.   Using the transformations given by     Eqs.~\ref{eq: SB1*} applied to the system \ref{eq: step3}, one can write :

\begin{equation} \left\{ \begin{array}{lll}
v_1-\frac{b}{1-c}&=&\gamma_1-\frac{b}{1-c} + \sum_{k=1}^7 S_k^{(1)} \sin kM' + \sum_{k=1}^7 C_k^{(1)} \cos kM' \vspace*{1mm}\\
&=&\frac{\gamma_1-c\gamma_1-\gamma_2+c\gamma_2}{1-c} + \sum_{k=1}^7 S_k^{(1)} \sin kM' +\sum_{k=1}^7 C_k^{(1)} \cos kM' \vspace*{1mm}\\
&=&\frac{\gamma_1-\gamma_2}{1-c} + \sum_{k=1}^7 S_k^{(1)} \sin kM' +\sum_{k=1}^7 C_k^{(1)} \cos kM' \vspace*{1mm}\\
v_2-\frac{b}{1-c}&=&\gamma_2-\frac{b}{1-c} - \sum_{k=1}^7 S_k^{(2)} \sin kM' +\sum_{k=1}^7 C_k^{(2)} \cos kM' \vspace*{1mm}\\
&=&\frac{\gamma_2-c\gamma_2-\gamma_2+c\gamma_1}{1-c} - \sum_{k=1}^7 S_k^{(2)} \sin kM' +\sum_{k=1}^7 C_k^{(2)} \cos kM' \vspace*{1mm}\\
&=&c\frac{\gamma_2-\gamma_1}{1-c} - \sum_{k=1}^7 S_k^{(2)} \sin kM' +\sum_{k=1}^7 C_k^{(2)} \cos kM' \\  
\end{array} \right. \label{eq: step5} \end{equation}

\begin{equation} \left\{ \begin{array}{lll}
v_1^*&=&\left(v_1-\frac{b}{1-c}\right)\sqrt{-c}= \frac{\sqrt{-c}}{1-c}\left(\gamma_1 -\gamma_2\right) + \sum_{k=1}^7 \sqrt{-c} \left( S_k^{(1)} \sin kM' + C_k^{(1)} \cos kM'  \right)\vspace*{1mm}\\
v_2^*&=&\left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}=-\frac{\sqrt{-c}}{1-c}\left(\gamma_1 -\gamma_2\right) - \sum_{k=1}^7 \frac{1}{\sqrt{-c}} \left( S_k^{(2)} \sin kM' +  C_k^{(2)} \cos kM'  \right)\\ \end{array} \right. \label{eq: step6} \end{equation}


However, from Eqs. \ref{eq: Sk} and \ref{eq: Ck} 
\begin{eqnarray} 
\sqrt{-c}S_k^{(1)}&=&-\mu e^{k-1} \sqrt{\frac{a_2}{a_1}}a_1 \sin i \left( \sin \o X_k \cos kM_0 + \cos \o Y_k \sin kM_0 \right) 
\vspace*{1mm}
 \\
&=&-\mu e^{k-1} \sqrt{a_1a_2} \sin i \left( \sin \o X_k \cos kM_0 + \cos \o Y_k \sin kM_0 \right)
\vspace*{1mm}
 \\
&=&-\mu e^{k-1} a^* \sin i \left( \sin \o X_k \cos kM_0 + \cos \o Y_k \sin kM_0 \right) = S_k^*
\vspace*{1mm}
 \\
\sqrt{-c}C_k^{(1)}&=&+\mu e^{k-1} \sqrt{\frac{a_2}{a_1}} a_1 \sin i \left( - \sin \o X_k \sin kM_0 + \cos \o Y_k \cos kM_0 \right)
\vspace*{1mm}
 \\
 &=&+\mu e^{k-1} \sqrt{a_1 a_2} \sin i \left( - \sin \o X_k \sin kM_0 + \cos \o Y_k \cos kM_0 \right)
\vspace*{1mm}
 \\
 &=&+\mu e^{k-1} a^* \sin i \left( - \sin \o X_k \sin kM_0 + \cos \o Y_k \cos kM_0 \right)= C_k^*
\end{eqnarray}

where the $S_k^*$ and  $C_k^*$ functions are identical to the $S_k^{(j)}$ and $C_k^{(j)}$ except that they concern a system of separation $a^*=\sqrt{a_1 a_2}$ instead of one of separation $a_j$. Similarly :
\begin{eqnarray} 
\frac{1}{\sqrt{-c}}S_k^{(2)} & = & -\mu e^{k-1}\sqrt{a_1a_2} \sin i \left( \sin \o X_k \cos kM_0 + \cos \o Y_k \sin kM_0 \right) = S_k^*\vspace*{1mm}\\
\frac{1}{\sqrt{-c}}C_k^{(2)} & = & +\mu e^{k-1}\sqrt{a_1a_2} \sin i \left( - \sin \o X_k \sin kM_0 + \cos \o Y_k \cos kM_0 \right)= C_k^*
\end{eqnarray}

Therefore, adopting :
\begin{equation}
\Gamma = \frac{\sqrt{-c}}{1-c}\left(\gamma_1 -\gamma_2\right) \end{equation}
one can finally achieve :
 
\begin{equation} \left\{ \begin{array}{lll}
v_1^*&=&\Gamma + \sum_{k=1}^7 S_k^{*} \sin kM' +\sum_{k=1}^7 C_k^{*} \cos kM'  \vspace*{1mm}\\
v_2^*&=&-\Gamma - \sum_{k=1}^7 S_k^{*} \sin kM' -\sum_{k=1}^7 C_k^{*} \cos kM'  \\ \end{array} \right. \label{eq: step8} \end{equation}
that forms a linear system of equations with unknowns $\Gamma$, $S_k^{*}$ and $C_k^{*}$. The coefficient $\sin kM'$ and $\cos kM'$ can indeed be easily computed from the observations date and a first guess of the period.
System \ref{eq: step8} can be rewritten in matrix form :

\begin{eqnarray} 
\left( \begin{array}{c} | \\ v_1^* \\ | \\ | \\ v_2^* \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{ccc}
| &   |      & | \\
1 & \sin kM' &\cos kM' \\ 
| &   |      & | \\
| &   |      & | \\
-1 & -\sin kM' & -\cos kM' \\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \Gamma \\ S_k^{*} \\ C_k^{*} \\\end{array} \right)\label{eq: step9}
\end{eqnarray} 
The above matrix is of dimension $2N \times (2k+1)$ where $N$ is the number of $(v_1,v_2)$ measurements and $k$ is the number of terms adopted in the harmonic analysis. Here we have presented the harmonic development with $k=7$, however, in practice, it is often reasonable to limit the analysis to $k =3$.
One can now solve the system \ref{eq: step9} and recover the values of
the  unknowns $\Gamma$, $S_k^{*}$ and $C_k^{*}$. From the latter
values, one has then to extract the information on the fake SB1
system. This is adressed in the next section.

\subsection{The approximate SB1 solution}
To derived the orbital parameters of the fake SB1 solution, we can
proceed exactly as explained by \citetalias{WHS67}. The demonstration
is provided here for the sake of completeness, although it is
identical to the one presented by \citetalias{WHS67}. However, one
should keep in mind that the signification of the symbols used are
different. Let us define:
\begin{equation} \begin{array}{lll}
b_1\cos\beta_1=Y_1\cos\o & \hspace*{1cm}& b_2\cos\beta_2=Y_2\cos\o \\
b_1\sin\beta_1=X_1\sin\o & \hspace*{1cm}& b_2\sin\beta_2=X_2\sin\o \\
\alpha_1=M_0+\beta_1     & \hspace*{1cm}& A_1^{*}=b_1\mu a^{*}\sin i\\
\alpha_2=2M_0+\beta_2    & \hspace*{1cm}& A_2^{*}=b_2\mu a^{*}e\sin i\\
\end{array}\end{equation} 
Since $X_1$, $X_2$, $Y_1$ and $Y_2$ are positive, $b_1$, $b_2$, $A^*_1$ and $A^*_2$ can also be chosen positive. One thus obtains :
\begin{eqnarray}
S_1^{*} & = & - a^{*} b_1 \mu \sin i \left( \sin\beta_1 \cos M_0 + \cos\beta_1 \sin M_0 \right) \nonumber\\
        & = & -A_1^{*}\sin\left(\beta_1+M_0\right) \nonumber\\
        & = & -A_1^{*}\sin\alpha_1 \\
S_2^{*} & = & -a^{*} b_2 e \mu \sin i \sin\left(\beta_2+2M_0\right) \nonumber\\
        & = & -A_2^{*}\sin\alpha_2 \\
C_1^{*} & = &   a^{*} b_1 \mu \sin i \left(-\sin\beta_1 \sin M_0 + \cos\beta_1 \cos M_0 \right) \nonumber\\
        & = &  A_1^{*}\cos\left(\beta_1+M_0\right) \nonumber\\
        & = &  A_1^{*}\cos\alpha_1 \\
C_2^{*} & = &  a^{*} b_2 e \mu \sin i \cos\left(\beta_2+2M_0\right) \nonumber\\
        & = &  A_2^{*}\cos\alpha_2 
\end{eqnarray}
Hence :
\begin{equation} \begin{array}{lllllll}
\tan \alpha_1 & = & - S_1^{*} / C_1^{*} &\hspace*{1cm} & \left(A_1^{*}\right)^2 & = & \left(C_1^{*}\right)^2 + \left(S_1^{*}\right)^2 \\
\tan \alpha_2 & = & - S_2^{*} / C_2^{*} & \hspace*{1cm}& \left(A_2^{*}\right)^2 & = & \left(C_2^{*}\right)^2 + \left(S_2^{*}\right)^2 
\end{array} \end{equation}
which yields $A^*_1$, $A^*_2$, $\alpha_1$ and $\alpha_2$.
The quantities $b_1$, $b_2$, $\beta_1$ and $\beta_2$ are computed by iteration, beginning with the following initial value for $\o$ and $e$:
\begin{equation}
\o^{(0)}=2\a_1-\a_2 \hspace*{0.5cm}\mathrm{and}\hspace*{0.5cm} e^{(0)}=A^*_2/A^*_1
\end{equation}
The preliminary $b_1^{(0)}$, $b_2^{(0)}$, $\beta_1^{(0)}$ and $\beta_2^{(0)}$  are computed from :
\begin{eqnarray}
 b_1^{(0)} & = & \left(Y_1^{(0)} \right)^2 \cos^2\o^{(0)} + \left(X_1^{(0)} \right)^2 \sin^2\o^{(0)} \\
 b_2^{(0)} & = & \left(Y_2^{(0)} \right)^2 \cos^2\o^{(0)} + \left(X_2^{(0)} \right)^2 \sin^2\o^{(0)} \\
\tan \b_1^{(0)} & = & X_1^{(0)} \tan \o^{(0)} / Y_1^{(0)}\\
\tan \b_2^{(0)} & = & X_2^{(0)} \tan \o^{(0)} / Y_2^{(0)}
\end{eqnarray}
The improved value $\o^{(1)}$ and $e^{(1)}$ are then found from :
\begin{eqnarray}
\o^{(1)} & = & 2\a_1-\a_2-2\left( \b_1^{(0)} - \o^{(0)} \right) + \left( \b_2^{(0)} - \o^{(0)} \right)\\
e^{(1)} & = & \frac{b_1^{(0)}A^*_2}{b_2^{(0)}A^*_1}
\end{eqnarray}
The process can be repeated starting from $\o^{(1)}$ and $e^{(1)}$ until the final values of $\o$, $e$, $b_1$, $b_2$, $\b_1$ and $\b_2$ are found. The physical parameters of the SB1 system can then be computed thanks to :
\begin{eqnarray}
a^* \sin i &=& \frac{A^*_1}{b_1 \mu} \\
M_0 &=& \a_2-\a_1+\b_1-\b_2 \\
T_0&=&t_0-\frac{M_0}{\mu} \\
\mathcal{K} &=&  \frac{\mu a^*\sin i}{\sqrt{1-e^2}} = \frac{A^*_1}{b_1\sqrt{1-e^2}}
\end{eqnarray}
According to \citetalias{WHS67}, {\it ``the method is very accurate for small eccentricities, but can be used for moderate and high eccentricities provided that the series \ref{eq: BesselC} and \ref{eq: BesselS} are carried further than just $k=2$. Hence, more unknowns than just $\Gamma$, $S_1$, $S_2$, $C_1$ and $C_2$ are involved. The only reason for including more terms is to obtain more accurate values for the first five. For eccentricities greater than about 0.6, the series \ref{eq: BesselC} and \ref{eq: BesselS} diverge. The problem is that the elements are being determined from a set of converging coefficient by a method using diverging series. It was found that reasonable values of the elements could be obtained from eccentricities as high as 0.83 solving for five or six S's and C's. Although the accuracy of such results are difficult to judge, these elements can nevertheless be used as the basis for a differential correction''.}

\subsection{The differential correction}
The previous stage yields approximate values for the parameters $e$, $\o$, $T_0$, $\Gamma$ and $\mathcal{K}$ of the fake SB1 system. For moderate or high eccentricity, it is necessary to apply a differential correction (D.C.) to improve the preliminary orbit. At this stage, the orbital period $P$ (or more precisely $\mu=2\pi/P$) can also be adjusted. The radial velocity equation of the SB1 system considered here can be written:
\begin{equation}
v^*=\pm\Gamma \pm \mathcal{K} \left( e \cos\o+ \cos(\o+\t) \right)
\end{equation}
where the ($+$) correspond to the primary transformed velocities and the ($-$) to the secondary ones. This equation hides five independent variables: $\Gamma$, $\mathcal{K}$, $\o$, $e$, $\mu$ and $T_0$. Let us derive $v^*$ with respect to each of them:
\begin{eqnarray}
\frac{\partial v^*}{\partial \Gamma }      & = & \pm 1
\vspace*{1mm}
\\
\frac{\partial v^*}{\partial \mathcal{K} } & = & \pm \left( e \cos\o+ \cos(\o+\t) \right)
\vspace*{1mm}
\\
\frac{\partial v^*}{\partial \o } & = & \mp \mathcal{K} \left( e \sin\o+ \sin(\o+\t) \right)
\vspace*{1mm}
\\
\frac{\partial v^*}{\partial \mu } & = & \mp \mathcal{K} \sin(\o+\t) \frac{d \t}{d E}\frac{d E}{d \phi}\frac{d \phi}{d \mu}
 \vspace*{1mm}
\\
                                   &   & \mathrm{since} \hspace*{0.5cm} \tan\frac{\t}{2}=\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2},  \hspace*{0.5cm}\frac{d\t}{d E}=\sqrt{\frac{1+e}{1-e}}\frac{\cos^2\t/2}{\cos^2E/2} = \sqrt{\frac{1+e}{1-e}}\frac{1-e}{1-e \cos E}
 \vspace*{1mm}
\\
                                   &   & \mathrm{since}\hspace*{0.5cm} \phi=E-e\sin E, \hspace*{0.5cm} \frac{d E}{d \phi} = \left( 1-e\cos E\right)^{-1} 
 \vspace*{1mm}
\\
                                   &   & \mathrm{since}\hspace*{0.5cm} \phi=\mu\left(t-T_0\right), \hspace*{0.5cm} \frac{d \phi}{d \mu} = \left(t-T_0\right)
 \vspace*{1mm}
\\
\frac{\partial v^*}{\partial \mu } & = & \mp \mathcal{K} \sin(\o+\t) \sqrt{\frac{1+e}{1-e}} \frac{1-e}{1-e \cos E} \frac{t-T_0}{1-e \cos E}
 \vspace*{1mm}
\\
                                   &   & \pm \mathcal{K} \sin(\o+\t) \sqrt{1-e^2} \left(T_0-t\right)\frac{\left(1+e\cos\t\right)^2}{\left(1-e^2\right)^2}
 \vspace*{1mm}
\\
                                   &   & \pm \mathcal{K} \sin(\o+\t)\left(1+e\cos\t\right)^2\left(1-e^2\right)^{-3/2}
 \vspace*{1mm}
\\
\frac{\partial v^*}{\partial e }   & = &  \pm \mathcal{K} \left( \cos\o-\sin(\o+\t)\frac{d\t}{de}\right)
 \vspace*{1mm}
\\
                                   &   & \mathrm{since}\hspace*{0.5cm} \tan\frac{\t}{2}=\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2},\hspace*{0.5cm}\frac{d\t}{de}\frac{1}{\cos^2\frac{\t}{2}}=\left( \frac{1+e}{1-e}\right)^{-1/2}\frac{1-e+1+e}{\left(1-e\right)^2}\tan\frac{E}{2}+\sqrt{\frac{1+e}{1-e}}\frac{1}{\cos^2\frac{E}{2}}\frac{dE}{de}
 \vspace*{1mm}
\\
                                   &   & \hspace*{4.0cm} = \sqrt{\frac{1-e}{1+e}}\frac{2}{\left( 1-e\right)^2}\tan\frac{E}{2}+\sqrt{\frac{1+e}{1-e}}\frac{1}{\cos^2\frac{E}{2}}\frac{dE}{de}
 \vspace*{1mm}
\\
                                   &   & \mathrm{since}\hspace*{0.5cm} \phi=E-e\sin E, \hspace*{0.5cm} \frac{d \phi}{d e}=0=\frac{d E}{d e}-\sin E -e \cos E\frac{dE}{de}
 \vspace*{1mm}
\\
                                   &   & \hspace*{4.0cm} \frac{dE}{de}=\frac{\sin E}{1-e\cos E}
 \vspace*{1mm}
\\
                                   &   & \mathrm{hence,} \hspace*{0.5cm}\frac{d\t}{de}=\frac{1-e}{1+e}\frac{2\cos^2\frac{\t}{2}}{\left(1-e\right)^2}\frac{\sin\frac{\t}{2}}{\cos\frac{\t}{2}}+\sqrt{\frac{1+e}{1-e}}\frac{\cos^2\frac{\t}{2}}{cos^2\frac{E}{2}}\frac{\sin E}{1-e\cos E}
 \vspace*{1mm}
\\
                                   &   & \hspace*{2.0cm} =\frac{\sin \t}{1-e^2}+\frac{\sin\t}{1-e\cos E}
 \vspace*{1mm}
\\
                                   &   & \hspace*{2.0cm} =\frac{\sin \t}{1-e^2}\left(1+\frac{1-e^2}{1-e^2}\left(1+e\cos \t\right)\right)
 \vspace*{1mm}
\\
                                   &   & \hspace*{2.0cm} =\frac{\sin \t}{1-e^2}\left(2+e\cos \t\right)
 \vspace*{1mm}
\\
\frac{\partial v^*}{\partial e }   & = &  \pm \mathcal{K}\left( \cos\o-\left(1-e^2\right)^{-1} \sin\t \sin(\t+\o)\left( 2+e\cos\t\right) \right)
 \vspace*{1mm}
\\
\frac{\partial v^*}{\partial T0 }  & = &  \mp \mathcal{K} \frac{d\t}{dE}\frac{dE}{d\phi}\frac{d\phi}{dT_0}
  \vspace*{1mm}
\\
                                   &   & \mathrm{since}\hspace*{0.5cm}\frac{d\phi}{dT_0}=-\mu 
 \vspace*{1mm}
\\
\frac{\partial v^*}{\partial T0 }  & = &  \pm \mathcal{K} \sin(\o+\t)\sqrt{\frac{1+e}{1-e}}\frac{1-e}{\left(1-e\cos E\right)^2}\mu
 \vspace*{1mm}
\\
                                   & = &  \pm \mathcal{K}  \left(1-e^2\right)^{-3/2} \sin(\o+\t) \mu \left(1+e\cos\t\right)^2
\end{eqnarray}
Hence one can write the linear equation system :
\begin{eqnarray} 
\left( \begin{array}{c} | \\ (O_k-C_k)_1* \\ | \\ | \\ (O_k-C_k)_2^* \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{cccccc}
   |   &   |    &   |    &   |    &   |    &   |   \\
q_1(k) & q_2(k) & q_3(k) & q_4(k) & q_5(k) & q_6(k)\\ 
   |   &   |    &   |    &   |    &   |    &   |   \\
     | &   |      & | \\
-q_1(k) & -q_2(k) & -q_3(k) & -q_4(k) & -q_5(k) & -q_6(k)\\ 
   |   &   |    &   |    &   |    &   |    &   |   \\
\end{array} \right)
\left(\begin{array}{c} \Delta\Gamma \\ \Delta\mathcal{K} \\ \Delta e  \\ \Delta \o  \\ \Delta T_0  \\ \Delta \mu  \\  \end{array} \right)
\end{eqnarray} 
with
\begin{eqnarray}
q_1(k)&=&1 
\vspace*{1mm}
\\
q_2(k)&=&e\cos\o+cos(\o+\t_k) 
\vspace*{1mm}
\\
q_3(k)&=&\mathcal{K}\left(\cos\o-\sin\t_k \sin(\o+\t_k)\frac{2+e\cos\t_k}{1-e^2} \right)
\vspace*{1mm}
\\
q_4(k)&=&-\mathcal{K} \left(e \sin\o+ \sin(\o+\t_k) \right)
\vspace*{1mm}
\\
q_5(k)&=&\mathcal{K}\mu \left(1+e\cos\t_k \right)^2\frac{\sin(\o+\t_k)}{(1-e^2)^{3/2}}
\vspace*{1mm}
\\
q_6(k)&=& \mathcal{K} \left(1+e\cos\t_k \right)^2\frac{\sin(\o+\t_k)}{(1-e^2)^{3/2}}  \left(T_0-hjd_k \right)
\vspace*{1mm}
\end{eqnarray}
that can be solved thanks to a linear least-square technique. This yields the corrections $\Delta\Gamma$, $\Delta\mathcal{K}$, $\Delta e$, $\Delta \o$, $\Delta T_0$ and eventually $\Delta \mu$ to be added to the previous determination of the orbital parameters. One can proceed recursively until sufficient convergence is reached. The errors on the orbital elements are simply the errors on the last step of the D.C. . This yields the final orbital solution for the fake SB1 system.  

\subsection{From the fake SB1 system to the real SB2 system} \label{finalel}
The parameters of the fake SB1 system obtained through the modified \citetalias{WHS67} method are $\Gamma$, $\mathcal{K}$, $e$, $\o$, $T_0$ and, eventually a new value for $\mu=2\pi/P$. An estimation of their 1-$\sigma$ uncertainties is also obtained. Hence it is a simple exercise to derive the parameters of the real SB2 system and of the associate uncertainties. 

\subsubsection{The common parameters}
Both the fake SB1 and the real SB2 system are characterized by the same eccentricity, periastron longitude and time of periastron passage so that the values of $e$, $\o$, $T_0$ (as well as of $\sigma_e$, $\sigma_\o$ and  $\sigma_{T_0}$) can be directly taken from the fake SB1 solution. The period $P$ is obtained from $P=2\pi/\mu$ while the associated error is given by 
\begin{eqnarray} 
\s_P&=&\frac{2\pi}{\mu^2}\s_\mu
\end{eqnarray}

\subsubsection{The systemic velocities}
From the Eqs. \ref{eq: SB1*}, one easily finds that :
\begin{eqnarray}
\gamma_1&=&\frac{\Gamma}{\sqrt{-c}}+\frac{b}{1-c}
\vspace*{1mm}
\\
\gamma_2&=&-\sqrt{-c}\Gamma+\frac{b}{1-c}
\end{eqnarray}
The errors are computed using the theory of error propagation:
\begin{equation} \begin{array}{lll}
      \hspace*{1.5cm} 
\frac{\partial \gamma_{1,2}}{\partial b}=\frac{1}{1-c} 
      & \hspace*{1.5cm}
\frac{\partial \gamma_{1}}{\partial \Gamma} =\frac{1}{\sqrt{-c}}
       &\hspace*{1.5cm}
\frac{\partial \gamma_{2}}{\partial \Gamma} = -\sqrt{-c} 
      \vspace*{2mm}
      \\
      \hspace*{1.5cm}
\frac{\partial \gamma_{1}}{\partial c} = \frac{\Gamma}{2\left(-c\right)^{3/2}}+\frac{b}{\left(1-c\right)^2} 
      & \hspace*{1.5cm}
\frac{\partial \gamma_{2}}{\partial c} = \frac{\Gamma}{2\sqrt{-c}}+\frac{b}{\left(1-c\right)^2} 
\end{array} \end{equation}
Hence
\begin{eqnarray}
\sb_{\gamma_1}&=&\frac{1}{-c}\sb_\Gamma+\frac{\sb_b}{\left(1-c\right)^2} + \left( \frac{\Gamma}{2\left(-c\right)^{3/2}}+\frac{b}{\left(1-c\right)^2}\right)^2\sb_c 
\vspace*{1mm}
 \\
\sb_{\gamma_2}&=&-c\sb_\Gamma+\frac{\sb_b}{\left(1-c\right)^2} + \left( \frac{\Gamma}{2\sqrt{-c}}+\frac{b}{\left(1-c\right)^2}\right)^2\sb_c 
\end{eqnarray}

\subsubsection{The semi-amplitudes of the RV curves}
\begin{eqnarray}
K_1&=&\frac{\mathcal{K}}{\sqrt{-c}}\\
K_2&=&\mathcal{K}\sqrt{-c}
\end{eqnarray}

\begin{equation}
\begin{array}{ll}
      \hspace*{1.5cm}
\frac{\partial K_1}{\partial \mathcal{K}}=\frac{1}{\sqrt{-c}} 
     & \hspace*{1.5cm}
\frac{\partial K_2}{\partial \mathcal{K}}=\sqrt{-c} 
     \vspace*{2mm}
     \\
      \hspace*{1.5cm}
\frac{\partial K_1}{\partial c}=\frac{\mathcal{K}}{2\left(-c\right)^{-3/2}} 
     & \hspace*{1.5cm}
\frac{\partial K_2}{\partial c}=-\frac{\mathcal{K}}{2\sqrt{-c}} 
 \end{array} \end{equation}

\begin{eqnarray}
\sb_{K_1}&=&\frac{\sb_\mathcal{K}}{-c}+\frac{\mathcal{K}^2}{4\left(-c\right)^3}\sb_c\\
\sb_{K_2}&=&-c\sb_\mathcal{K}+\frac{\mathcal{K}^2}{4\left(-c\right)}\sb_c
\end{eqnarray}
Note also that 
\begin{eqnarray}
\sb_{K_1}&=&\frac{1}{\left( -c \right)^2}\left( -c\sb_\mathcal{K}+\frac{\mathcal{K}^2}{4\left(-c\right)}\sb_c\right)=\frac{\sb_{K_2}}{c^2}
\end{eqnarray}
As $c=K_2/K_1$, one obtains the previously announced results:
\begin{equation}
\frac{\sb_{K_1}}{K_1}=\frac{\sb_{K_2}}{K_2}
\end{equation}

\subsubsection{Other interesting parameters}
 
\begin{itemize}
 \item $a_1 \sin i$ and $a_2 \sin i$ 
\begin{equation}
a_1\sin i=\frac{K_1\sqrt{1-e^2}}{\mu}=\frac{\left(\mathcal{K}/\sqrt{-c}\right)\sqrt{1-e^2}}{\mu}
\end{equation}
\begin{equation}\begin{array} {lll}
      \hspace*{1.5cm}
\frac{\partial a_1\sin i}{\partial e}=-\frac{e}{1-e^2}a_1\sin i 
     & \hspace*{1.5cm}
\frac{\partial a_1\sin i}{\partial c}=\frac{a_1\sin i}{-2c}
      \hspace*{1.5cm}
     \frac{\partial a_1\sin i}{\partial \mathcal{K}}=\frac{a_1\sin i}{\mathcal{K}} 
\end{array} \end{equation}
\begin{equation}
\sb_{a_1\sin i}=\left(\frac{e^2}{\left(1-e^2\right)^2}\sb_e + \frac{\sb_{K_1}}{K_1^2} \right) \left( a_1\sin i\right)^2
\end{equation}
Similarly for $ a_2\sin i=\frac{K_2\sqrt{1-e^2}}{\mu}$ :

\begin{equation}
\sb_{a_2\sin i}=\left(\frac{e^2}{\left(1-e^2\right)^2}\sb_e + \frac{\sb_{K_2}}{K_2^2} \right) \left( a_2\sin i\right)^2
\end{equation}
If $\mu$ has been included in the D.C., one can have to add an additional contribution to the above error estimate: $\left(\frac{\s_\mu}{\mu}a_j\sin i\right)^2$.\vspace*{.4cm}\\

 \item $M_1 \sin^3 i$ and $M_2 \sin^3 i$ 
\begin{equation}
M_1 \sin^3 i = \frac{P}{2\pi G}K_2^3\left(1-e^2\right)^{3/2}\left(1-\frac{1}{c}\right)^2=\frac{P}{2\pi G}\left(1-e^2\right)^{3/2}\frac{\mathcal{K}^3}{\sqrt{-c}}\left(c-1\right)^2
\end{equation}
\begin{equation}
M_2 \sin^3 i = \frac{P}{2\pi G}K_1^3\left(1-e^2\right)^{3/2}\left(1-\frac{1}{c}\right)^2=\frac{P}{2\pi G}\left(1-e^2\right)^{3/2}\frac{\mathcal{K}^3}{\left(-c\right)^{3/2}}\left(c-1\right)^2
\end{equation}
\begin{equation}\begin{array}{lll}
      \hspace*{1.2cm}
\frac{\partial M_1 \sin^3 i}{\partial c} = \left(\frac{1}{-2c}+\frac{2}{c-1}\right)M_1 \sin^3 i &
      \hspace*{1.2cm}
\frac{\partial M_j \sin^3 i}{\partial\mathcal{K}} = \frac{3}{\mathcal{K}}M_j \sin^3 i 
     & \hspace*{1.2cm}
\frac{\partial M_j \sin^3 i}{\partial e} = \frac{3e}{1-e^2}M_j \sin^3 i 
\\
     \hspace*{1.2cm}
\frac{\partial M_2 \sin^3 i}{\partial c} = \left(\frac{-3}{2c}+\frac{2}{c-1}\right)M_2 \sin^3 i 
     & \hspace*{1.2cm}
\frac{\partial M_j \sin^3 i}{\partial \mu} = \frac{M_j \sin^3 i}{\mu} 
 \end{array}\end{equation}
\begin{equation}
\sb_{ M_1 \sin^3 i}=\left( \frac{9\sb_\mathcal{K}}{\mathcal{K}^2}+ \frac{9e^2\sb_e}{\left(1-e^2\right)^2} +  \left( \frac{1}{-2c}+\frac{2}{c-1}\right)^2 \right) \left( M_1 \sin^3 i\right)^2
\end{equation}
\begin{equation}
\sb_{ M_2 \sin^3 i}=\left( \frac{9\sb_\mathcal{K}}{\mathcal{K}^2}+ \frac{9e^2\sb_e}{\left(1-e^2\right)^2} +  \left(\frac{-3}{2c}+\frac{2}{c-1}\right)^2\right)\left( M_2 \sin^3 i\right)^2
\end{equation}
plus, eventually, the period term: $\left( \frac{M_j \sin^3 i}{\mu}\s_\mu\right)^2$.\vspace*{.4cm}\\

 \item $\frac{R_\mathrm{RL,1}}{a_1+a_2}$ and $\frac{R_\mathrm{RL,2}}{a_1+a_2}$ \\

Let us write $q_1=-c$ and $q_2=-1/c$, hence $\s_{q_1}=\s_c$ and $\s_{q_2}=\s_c/c^2$. Using \citet{Egg83}, one can write :
\begin{equation}
\frac{R_\mathrm{RL,1}}{a_1+a_2}=\frac{0.49q_1^{2/3}}{0.6q_1^{2/3}+\ln\left(1+q_1^{1/3}\right)} 
\end{equation}
\begin{equation}
\frac{R_\mathrm{RL,2}}{a_1+a_2}=\frac{0.49q_2^{2/3}}{0.6q_2^{2/3}+\ln\left(1+q_2^{1/3}\right)} 
\end{equation}
 The respective error are given by :
\begin{equation}
\s_{\frac{R_\mathrm{RL,j}}{a_1+a_2}}=\frac{49}{300}\frac{|2q_j^{-1/3}\ln\left(1+q_j^{1/3}\right)-\left(1+q_j^{1/3}\right)^{-1}|}{\left(0.6q^{2/3}+\ln\left(1+q_j^{1/3}\right)\right)^2} \s_{q_j}
\end{equation}\vspace*{.4cm}


 \item $f_\mathrm{mass,1}$ and $f_\mathrm{mass,2}$

\begin{equation}
f_\mathrm{mass,j}=\frac{P}{2\pi G}\left(1-e^2\right)^{3/2}K_j^3
\end{equation}
\begin{eqnarray}
f_\mathrm{mass,1}&=&\frac{P}{2\pi G}\left(1-e^2\right)^{3/2}\mathcal{K}^3\left(-c\right)^{-3/2}
\\
f_\mathrm{mass,2}&=&\frac{P}{2\pi G}\left(1-e^2\right)^{3/2}\mathcal{K}^3\left(-c\right)^{3/2}
\end{eqnarray}
\begin{equation}\begin{array}{ll}
\hspace*{1.5cm}\frac{\partial f_\mathrm{mass,j}}{\partial P}=\frac{f_\mathrm{mass,j}}{P}
&
\hspace*{1.5cm}\frac{\partial f_\mathrm{mass,j}}{\partial \mathcal{K}}=3\frac{f_\mathrm{mass,j}}{\mathcal{K}}
\\
\hspace*{1.5cm}\frac{\partial f_\mathrm{mass,j}}{\partial e}=3\frac{f_\mathrm{mass,j}}{1-e^2}
&
\hspace*{1.5cm}\frac{\partial f_\mathrm{mass,j}}{\partial c}=-\frac{3}{2}\frac{f_\mathrm{mass,j}}{\left(-1\right)^j c}\\
\end{array}\end{equation}
\begin{equation}
\sb_{f_\mathrm{mass,j}}=\left( \frac{\sb_P}{P^2}+\frac{9e^2\sb_e}{\left(1-e^2\right)}+\frac{9\sb_\mathcal{K}}{\mathcal{K}^2}+\frac{9\sb_c}{4c^2}\right) f^2_\mathrm{mass,j}
\end{equation}
where the period term should only be included if this parameter is improved by the D.C. .

\end{itemize}
 
\section{A synopsis of the code}
The code to derive the eccentric SB2 orbital solution is structured in several important subroutines. These are :
\begin{enumerate}
\item[-] READDATA : reads the input data file,
\item[-] ORBIT : computes the SB2 eccentric orbital solution,
\item[-] MODEL : computes the phase for each measurements as well as the $O-C$ deviation and the final $\chi^2$ corresponding to the best fit solution derived in ORBIT,
\item[-] VELCURVE : computes the best fit RV-curves according to the best fit solutions found in ORBIT.
\end{enumerate}

The ORBIT subroutine is itself structured in various elements :
\begin{enumerate}
\item[-] MASSRATIO : for SB2 systems, performs the orthogonal regression and builds the {\it fake} SB1 system,
\item[-] INIT : builds the linear system \ref{eq: step9} to be inverted,
\item[-] APPROXEL : computes the approximate orbital solution,
\item[-] DC : given a criterion, performs the D.C. recursively until convergence is reached, then computes the improved orbital elements of the handled SB1 system,
\item[-] FINALEL :  derives the final SB2 parameters from the {\it fake} SB1 final solution.\end{enumerate}

\subsection{Reading the input file and input parameters}
The following parameters are necessary inputs. They are read either from the main program or from the READDATA routine.
\begin{tabular}{rl}
$isw1$ & Allows for an adjustment of the period ? (yes=1) \\
$isw2$ & In the case no D.C. is needed, forced it anyway ? (yes=1)\\
$isw4$ & Number of indirect (or additional) parameters to be considered in the fit ?\\
$mlabel$ & Maximum number of D.C. steps \\
$m$ & Number of coefficient in the harmonic analysis \\
$id$ & Object name \\
$P$ & Initial guess for the orbital period \\
$hjdl$ & Approximate epoch near which $T_0$ will be computed\\
$s$ & $ s = \s_y / \s_x = \s_{v_2} / \s_{v_1}$ : ratio of the secondary to primary uncertainties \\
$np$ & Number of couples $\left(v_1,v_2\right)$\\
$\left(hjd,vp,wtp,vs,wts \right)_{k=1}^{np}$ & Julian date, $v_1$ and weight , $v_2$ and weight \\
\end{tabular}

\subsection{Building the fake SB1 system: MASSRATIO}
The routine is performed for SB2 system only. 
It first selects all the points for which $wtp(k)=wts(k)$ in the data list and passes them to the orthogonal regression routine. Let us assume that there are $n_\mathrm{eff}$ such couple $(v_1,v_2)_k$. The points for which $wtp(k)\ne wts(k)$ are not considered in the orthogonal regression (QUESTION: Should we also reject all points with a null weight?).  

\subsubsection{The orthogonal regression by itself}
Let us note in this section $w_k=wtp(k)=wts(k)$. The weight are first normalized so that $\sum_{k=1}^{n_\mathrm{eff}}w_k=n_\mathrm{eff}$. 

\begin{eqnarray}
\rho&=&\left(S_y-s^2S_x\right)^2+4s^2C_{xy}^2\\
sq1=c&=&\frac{\left(S_y-s^2S_x\right)+\sqrt{\rho}}{2C_{x_y}} \\
sq2=b&=&\y-c\x \\
\mathcal{F}&=&\frac{S_y+s^2S_x}{C_{xy}^2}\\
\sb_c&=&\frac{\chi^2_\nu}{n_\mathrm{eff}}c^2\mathcal{F}\\
\sb_b&=&\frac{\chi^2_\nu}{n_\mathrm{eff}}\left(s^2+c^2+c^2\x^2\mathcal{F}\right)\\
\chi^2_\nu&=&\frac{n_\mathrm{eff} \chi^2}{(n_\mathrm{eff}-2)\sum_k w_k} \\
\chi^2&=&\sum_k w_k \frac{\left(y_k-cx_k-b\right)^2}{s^2+c^2}
\end{eqnarray}

The code also computes the linear correlation coefficient
\begin{equation}
r=\frac{C_{xy}}{S_xS_y}
\end{equation}

\subsubsection{An absolute scale for the weights} \label{ssect: absscale}
From here all the data points are again considered, even those with $wtp(k) \ne wts(k)$.
We now use the best fit relation to bring the weight to an absolute scale:
\begin{eqnarray}
 \sb&\approx&\chi^2_\nu\\
 wtp(k)&=&\frac{wtp(k)}{\sb} \\
 wts(k)&=&\frac{wts(k)}{\sb}/s^2
\end{eqnarray}
so that $\sb_{v_p(k)}=1/wtp(k)$ and $\sb_{v_s(k)}=1/wts(k)$

\subsubsection{Initialization of the fake SB1 system}
Let us write $n = 2np$. The transformation of the SB2 velocities  into equivalent SB1$^*$ velocities yields, for $k=1...n/2$ :
\begin{eqnarray}
v^*(k)&=&\left(v_1(k)-\frac{b}{1-c}\right)\sqrt{-c} \\
\left(w^*(k)\right)^{-1}&=&-c\sb_{v_p(k)}-\frac{c}{\left(1-c\right)^2}\sb_b+\left(\frac{v^*(k)}{2c}-\frac{b\sqrt{-c}}{\left(1-c\right)^2}\right)^2\sb_c
\end{eqnarray}
 For $k=(n/2+1)...n$ :
\begin{eqnarray}
v^*(k)&=&\left(v_2(l)-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}} \hspace*{0.5cm}\mathrm{with }\hspace*{0.5cm} l=k-n/2\\
\left(w^*(k)\right)^{-1}&=&-\frac{\sb_{v_s(l)}}{c}-\frac{\sb_b}{c(1-c^2)}+\left(\frac{v^*(k)}{2c}+\frac{b}{\sqrt{-c}\left(1-c\right)^2}\right)^2\sb_c
\end{eqnarray}
QUESTION: Should wee account for $\s_c$ and $\s_b$ while converting $v_1$ and $v_2$ into $v^*$ (as done here above). If not, then the previous equations can be rewritten :
\begin{eqnarray}
\left(w^*(k)\right)^{-1}&=&-c\sb_{v_p(k)} \hspace*{0.5cm}  \mathrm{for}\hspace*{0.5cm} k=1...n/2 \vspace*{1mm} \\
\left(w^*(k)\right)^{-1}&=&-\frac{\sb_{v_s(l)}}{c} \hspace*{0.5cm} \mathrm{for}\hspace*{0.5cm} k=(n/2+1)...n  \hspace*{0.5cm} \mathrm{and}\hspace*{0.5cm}l=k-n/2
\end{eqnarray}
so that no absolute scaling (Sect.~\ref{ssect: absscale}) needs to be defined. Hence, one could simply use :
\begin{eqnarray}
w^*(k) &=&\frac{ wtp(k)}{-c} \hspace*{0.5cm} \mathrm{for}\hspace*{0.5cm} k=1...n/2 \vspace*{1mm} \\
w^*(k) &=&-c\ wts(l)/s^2 \hspace*{0.5cm} \mathrm{for}\hspace*{0.5cm} k=(n/2+1)...n\hspace*{0.5cm} \mathrm{and}\hspace*{0.5cm}l=k-n/2
\end{eqnarray}

\subsection{Initializing the Wilsing-Russell method: INIT}
The weights are first renormalized accounting this time for all the (transformed) RVs, so that $\sum_{k=1}^{2np} wt(k)=2np=n$. Then, from the $\{v^*(k),w^*(k)\}$, $k=1...n$, one constructs the matrix form of the linear system 

\begin{eqnarray} 
\left( \begin{array}{c} | \\ v^*(k) \\ | \\ | \\ v^*(k) \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{ccc}
| &   |      & | \\
1 & \sin hM'(k) &\cos hM'(k) \\ 
| &   |      & | \\
| &   |      & | \\
-1 & -\sin hM'(k) & -\cos hM'(k) \\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \Gamma \\ S_h^{*} \\ C_h^{*} \\\end{array} \right)
\end{eqnarray} 
with $h=1...m$ ($m$ being the number of coefficients in the harmonic analysis). The separation between the upper and lower part of the matrix happens between the lines $n$p and $np+1$. 
The system  is then solved using the Numerical Recipes routine GAUSSJ. This yields the best fit parameters $\Gamma$, $S^*_k$ and $C^*_k$.
The code then computes the QUALITY of the least-square fit :
\begin{eqnarray}\chi^2&=&\sum_{k=1}^{np} w^*(k) \left( v^*(k) - \Gamma - \sum_{h=1}^{h=m} \sin hM'(k)\  S^*_k -\sum_{h=1}^{h=m} \cos hM'(k)\  C^*_h \right)^2 
\\
& + & \sum_{k=np+1}^{2np} w^*(k) \left( v^*(k) + \Gamma + \sum_{h=1}^{h=m} \sin hM'(k)\  S^*_h + \sum_{h=1}^{h=m} \cos hM'(k)\  C^*_h \right)^2  
\end{eqnarray}
Hence 
\begin{eqnarray}
\chi^2_\nu&=&\frac{\chi^2}{2np-(2m+1)} \label{eq: quality}  \\
r_1^2&=&\frac{\chi^2}{\left(\sum_k w_k\right) \left(1-\frac{2m+1}{2np}\right)} = \chi^2_\nu
\end{eqnarray}



\subsection{Approximated SB1 orbital elements: APPROXEL}
The approximate elements $\mathcal{K}$, $e$, $\o$ and $T_0$ are then estimated from  $S^*_1$ and $C^*_1$, $S^*_2$ and $C^*_2$. If $e\le5\ 10^{-7}$, the program assumes a circular orbit (case not checked). If $e$ is larger, the code performs an iterative process to improve $e$ and $\o$ (stop criterion: $\Delta^2e \le 2.5\ 10^{-13}$) following exactly WHS67. Finally, it computes the parameters related to the SB1 system as well as their errors. The latter are computed by propagation from the diagonal elements of the COVAR matrix properly scales using the $\chi^2_\nu$ (Eq. \ref{eq: quality}) of the fit .

\begin{equation}
A_j^2=C^2_j+S^2_j
\end{equation}
\begin{equation}
\a_j=\left\{ \begin{array} {lll}   \arctan(-S_j/C_j) & if & C_j > 0 \vspace*{1mm}\\
                                   \arctan(-S_j/C_j)+\pi  & if & C_j < 0 \vspace*{1mm}\\
                                   \frac{\pi}{2}  & if & C_j = 0 \hspace*{2mm}\mathrm{and}\hspace*{2mm} S_j \le 0 \vspace*{1mm}\\
                                   \frac{3\pi}{2}  & if & C_j = 0 \hspace*{2mm}\mathrm{and}\hspace*{2mm} S_j > 0 
\end{array} \right. 
\end{equation}

\begin{eqnarray}
\sb_{A_j}&=&\frac{C_j^2 \sb_{C_j}+S_j^2 \sb_{S_j}}{A^2_j}
\vspace*{1mm}\\
\sb_{\a_j} &=&\frac{C_j^2 \sb_{S_j}+S_j^2 \sb_{C_j}}{A^4_j}
\vspace*{1mm}\\
a^* \sin i &=& \frac{A_1}{b_1\mu}\vspace*{1mm}
\vspace*{1mm}\\
T_0 &=& hjdl-\frac{\a_2-\a_1+\b_1-\b_2}{\mu} 
\vspace*{1mm} \\
\mathcal{K} &=& \mu(a^* \sin i)/\sqrt{1-e^2} 
\vspace*{1mm}\\
\sb_e&=&\frac{e^2 \sb_{A_1}+\sb_{A_2}}{A^2_1}\vspace*{1mm} \\
\sb_\o &=& 4\sb_{\a_1}+\sb_{\a_2}\vspace*{1mm} \\
\sb_\mathcal{K}&=&\sb_{A_1} \vspace*{1mm}\\
\sb_{a^* \sin i}&=&\sb_\mathcal{K}/\mu^2 \vspace*{1mm}\\
\sb_{T_0} &=& \frac{\sb_{\a_1}+\sb_{\a_2}}{\mu}
\end{eqnarray}

\subsection{The differential correction: DC}
If both $e>2\s_e$ and $ e> 0.03$, the program enters a differential correction (D.C.) phase. 
If the option $isw2=1$ has been chosen, the D.C. is applied whatever the result of the above criteria.
If allows by the user ($isw1=1$), the D.C. will also improve the guess value of the period (more precisely, it will improve $\mu$).
First the true anomaly $\t_k$ corresponding to the approximate solution is computed for each observation as well as $(O-C)_k=v^*_k-v_\mathrm{th}(\t_k)$. Then the following linear system is built :

\begin{eqnarray} 
\left( \begin{array}{c} | \\ (O-C)_k^* \\ | \\ | \\ (O-C)_k^* \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{cccccc}
| &   |      & | \\
q_1(k) & q_2(k) & q_3(k) & q_4(k) & q_5(k) & q_6(k)\\ 
| &   |      & | \\
| &   |      & | \\
-q_1(k) & -q_2(k) & -q_3(k) & -q_4(k) & -q_5(k) & -q_6(k)\\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \Delta\Gamma \\ \Delta\mathcal{K} \\ \Delta e  \\ \Delta \o  \\ \Delta T_0  \\ \Delta \mu  \\  \end{array} \right)
\end{eqnarray} 
with
\begin{eqnarray}
q_1(k)&=&1 
\vspace*{1mm}
\\
q_2(k)&=&e\cos\o+cos(\o+\t_k) 
\vspace*{1mm}
\\
q_3(k)&=&\mathcal{K}\left(\cos\o-\sin\t_k \sin(\o+\t_k)\frac{2+e\cos\t_k}{1-e^2} \right)
\vspace*{1mm}
\\
q_4(k)&=&-\mathcal{K} \left(e \sin\o+ \sin(\o+\t_k) \right)
\vspace*{1mm}
\\
q_5(k)&=&\mathcal{K}\mu \left(1+e\cos\t_k \right)^2\frac{\sin(\o+\t_k)}{(1-e^2)^{3/2}}
\vspace*{1mm}
\\
q_6(k)&=& \mathcal{K} \left(1+e\cos\t_k \right)^2\frac{\sin(\o+\t_k)}{(1-e^2)^{3/2}}  \left(T_0-hjd_k \right)
\vspace*{1mm}
\end{eqnarray}
The separation between the upper and lower part of the matrix occurs again between the lines $np$ and $np+1$. 
The matrix is inverted by means of the GAUSSJ routine and the
differential corrections are found.  The quality of the fit is
estimated by the QUALITY2 routine which computes the individual
$(O-C)_k$ as well as the merit function :
\begin{eqnarray}
\chi^2&=&\sum_{k=1}^{k=np} w_k \left((O-C)_k- \Delta\Gamma - q_2(k) \Delta\mathcal{K}  - q_3(k)\Delta e + q_4(k)\Delta \o - q_5(k) \Delta T_0 - q_6(k)\Delta \mu \right)^2  \\
&&+\sum_{k=np+1}^{k=2np} w_k \left((O-C)_k+ \Delta\Gamma + q_2(k) \Delta\mathcal{K}  + q_3(k)\Delta e - q_4(k)\Delta \o + q_5(k)  \Delta T_0 + q_6(k)\Delta \mu \right)^2 
\end{eqnarray} It follows that
\begin{eqnarray}
r_1^2 \hspace*{2mm} = &\chi^2_\nu \hspace*{2mm} =&\frac{\chi^2}{\nu}
\end{eqnarray}
where $\nu=N-6$ if the period (actually, $\mu=2\pi/P$) is included in the D.C.; $\nu=N-5$ otherwise (and $\Delta \mu$ is then ignored). 
The improved parameters $\Gamma$, $\mathcal{K}$, $e$, $\o$, $T_0$ and
eventually $\mu$ are simply computed by adding the various corrections found. 
The uncertainties on the improved parameters are taken as the square-root of the diagonal elements of the COVAR matrix times the square root of the $\chi^2_\nu$. 
The D.C. is performed again and again as long as the maximum number of  iterations allowed ($mlabel$) has not been reached and as long as the convergence criterion is not matched. The latter requires that $\Delta^2 \mathcal{K}>\sb_\mathcal{K}/9$,
where $\Delta \mathcal{K}$ is the D.C. step along the $\mathcal{K}$ direction while $\s_\mathcal{K}$ is the uncertainty $\Delta \mathcal{K}$.


\subsection{Computing the final elements: FINALEL}

 The FINALEL routine computes the parameters of the true SB2 and their 1-$\sigma$ dispersion.  It also computes various interesting parameters such as the minimal masses and the relative Roche lobe radii. The adopted formula have been previously given in Sect.~\ref{finalel}.


\subsection{Estimating the quality of the final SB2 system: MODELE}
The routine MODELE computes the $O-C$ and the $\chi^2$ with respect to the final SB2 solution:
\begin{eqnarray}
\chi_\mathrm{p}^2&=& \sum_{k=1}^{np}wtp(k)\left(vp(k)-\gamma_1-K_1\left(e\cos\o+\cos\left(\o+\t_k\right)\right)\right)^2 
\vspace*{1mm}
\\
\chi_\mathrm{s}^2&=& \sum_{k=1}^{np}wts(k)\left(vs(k)-\gamma_2+K_2\left(e\cos\o+\cos\left(\o+\t_k\right)\right)\right)^2 
\vspace*{1mm}
\\
\chi^2_\mathrm{SB2}&=&\chi_\mathrm{p}^2+\chi_\mathrm{s}^2
\end{eqnarray}
\begin{eqnarray}
r_\mathrm{SB2}^2&=&\frac{n\chi_\mathrm{SB2}^2 }{\nu_\mathrm{SB2} \sum_k(wtp_k+wts_k) }
\vspace*{1mm}
\\
r_\mathrm{p}^2&=&\frac{np\chi_\mathrm{p}^2 }{\nu_\mathrm{sgl} \sum_k wtp_k }
\vspace*{1mm}
\\
r_\mathrm{s}^2&=&\frac{np\chi_\mathrm{s}^2 }{\nu_\mathrm{sgl} \sum_k wts_k }
\end{eqnarray}
where $n=2np$ and $\nu_\mathrm{SB2}=n-8-isw4$ if the period has been adjusted; $\nu_\mathrm{SB2}=n-7-isw4$ otherwise. Finally, $\nu_\mathrm{sgl}=np-6-isw4$ if the period has been adjusted; $\nu_\mathrm{sgl}=np-5-isw4$ otherwise.
The $sv2$'s written in the result file correspond this time to $\chi^2$ (by opposition to $\chi^2_\nu$ for the results of the Wilsing-Russell method and of the D.C.).


\subsection{Computing the theoretical radial velocity curve(s): VELCURVE}
Finally, the VELCURVE routine computes the best fit RV curve(s) for the real SB1 or SB2 system with a step of 0.01 in phase.




\bsp

\label{lastpage}

\end{document}
