% mn2esample.tex
%
% v2.1 released 22nd May 2002 (G. Hutton)


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


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

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

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

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

\def\rsol{R$_{\sun}$}
\def\msol{M$_{\sun}$}

\def\rstar{R$_{\star}$}
\def\mstar{M$_{\star}$}

% 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}
\def\WHS67{\citetalias{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 Orbital Solution Package}
\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 the Li\`ege Orbital Solution Package (LOSP), a numerical
package that aims at computing the orbital parameters of spectroscopic binaries. The package deals with SB1 and SB2 systems and is able to adjust either circular or eccentric orbits through a weighted fit. The SB1 eccentric fit relies on the Wolfe et al.\ method. The SB2 eccentric orbits are adjusted using the Li\`ege SB2 Orbital Solution Algorithm (LOSA), that uses an orthogonal regression technique and a modified version of the Wolfe at al.\ algorithm to derive self-consistent solutions for both components of the system. The SB1 circular orbits are fitted through a linear least square technique. Finally, the same approach as the one adopted in LOSA as been adapted to SB2 circular systems. It allows to generate self-consistent solutions for both components. As an additional capability, LOSP can perform an exploration of the parameter space along the period axis and, eventually, along the ratio of the secondary to primary radial velocities uncertainties. Beyond the standard error computations, LOSP further provides the opportunity to perform Monte-Carlo simulations on the basis of the best-fit solution. This option allows an independent and robust check of the accuracy of the determined parameters in the four considered cases (SB1 vs. SB2, circular vs. eccentric orbit). 
\end{abstract}

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

\section{Introduction}

It has long been known that most stars in the Galaxy belong to  binary
 or multiple systems. Among the various kinds of systems,
 spectroscopic binaries are considered as privileged laboratories to
 constrain the stars fundamental parameters. Through their
 gravitational interaction, they potentially offers to accurately
 measure the masses. In practice, the situation is unfortunately less
 simple as accurate masses and radii can only be obtained for SB2
 eclipsing binaries. SB2 non eclipsing systems only yield the minimal
 masses $M~\sin^3 i$ of each components (where $i$ is the orbital
 inclination). Still the mass ratio of these systems can be measured
 independently of the inclination. Finally SB1 systems only provided a
 limited information on the stars weight through the so-called mass
 function $f_\mathrm{mass}=\frac{P K^3}{2\pi G}
 \left(1-e^2\right)^{3/2}$, where the  notations used are defined in Table \ref{tab: orbparam}. However, even SB1 systems bring important constraints on the formation and evolution of binaries through e.g. the period vs. eccentricity diagram. Beside the intrinsic quality of the data, an important aspect of orbital solution fitting lies in an accurate estimation of the uncertainties that spoil the derived orbital parameters. 

In this work, we present a package that allows to compute orbital solutions of spectroscopic binaries. Developed at the Li\`ege University (Belgium), the package indifferently addresses SB1 or SB2 systems and deals with both circular or eccentric orbits. 
It only requires a first guess of the location of the orbital period in the parameter space, which can  easily be obtained through an appropriate Fourier analysis \citep[ see also \citealt{GRR01} for comments]{HMM85}. The preliminary period value may then be improved either by a differential correction method (for eccentric systems) or through a global minimization of the $\chi^2$ along this particular degree of freedom  (both for eccentric or non eccentric systems). Finally, a particular attention has been given to the estimates of the errors. Beside a first estimate obtained through the theory of error propagation, the package also allows to perform Monte-Carlo simulations to estimate the confidence limits on the best-fit parameters. 

\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 for 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 \\
$f_\mathrm{mass_{1,2}}$ & Mass function associated with component 1, 2 \\
\hline
\end{tabular}
\end{table*}

\section{An overview of the package}

As mentioned above, LOSP mainly handles four cases: SB1 and SB2 eccentric orbits as well as  SB1 and SB2 circular orbits. We now briefly present the underlining methods adopted in each of the four cases. 

\subsection{SB1 eccentric orbit} \label{ssect: sb1}
The adopted algorithm is that of \citet{WHS67}, that relies on the Wilsing-Russell method \citep{Wilsing1894,Russell1902} followed by a differential correction (D.C.). It requires a first guess of the orbital period that can, on request, be improved by the D.C. As an additional capability, LOSP proposes to explore the period axis in a given interval and with a given step and return the period value that minimized the $\chi^2$ in the considered interval.  The \WHS67\ method offers the advantage that no a priori estimate of the location of the solution in the parameter space is required nor is an exploration of the 5-dimension  parameter space. \citet{WHS67} stated that the method is robust up to eccentricities about 0.8 though no detailed testing was presented.

\subsection{SB2 eccentric orbit} \label{ssect: sb2}

A dedicated algorithm for SB2 eccentric orbit has been developed and
is known as the Li\`ege SB2 orbital Solution Algorithm. Roughly
speaking, it uses an orthogonal regression technique to fit the $v_2$
vs. $v_1$ relationship (which directly provides the mass ratio) qnd to transform the SB2 RV data set into a {\it fake} but equivalent SB1 RV data set. The latter is characterized by the same period, eccentricity, longitude of periastron and time of periastron passage as the real SB2 system. The methods then applies a modified version of the \WHS67\ method to provide the orbital elements of the equivalent SB1 system from which the the final SB2 solution can be derived. The Li\`ege SB2 algorithm offers the advantage to provide self-consistent SB2 solutions while preserving much of the properties of the  \WHS67\ method. In particular, only a first guess of the period is needed and no fastidious exploration of the 8-dimension parameters space is required. Details of the algorithm are presented in \citet{LOSA}.

\subsection{SB1 circular orbit} \label{ssect: sb1ne}

For SB1 circular orbits, LOSP performs a linear least square fit of
the system RV equation:
\begin{equation}
v(t_i)=\g+K\sin \p_i, \label{eq: SB1ne}
\end{equation}
where $\p_i=2\pi \frac{t_i-T_0}{P}$ is the phase angle measured from the time of primary conjunction. This problem has three unknowns: $\g$, $K$ and $T_0$. Let us define $\p'=2\pi\frac{t_i-t_\mathrm{arb}}{P}$ where $t_\mathrm{arb}$ is an arbitrary time. The phase angle can now be rewritten: $\p_i=\p_i'+\o =2\pi\frac{t_i-t_\mathrm{arb}}{P}+2\pi\frac{t_\mathrm{arb}-T_0}{P}$. Hence Eq.~\ref{eq: SB1ne} can be rewritten:
\begin{equation}
v(t_i)=\g+K \cos \o \sin \p'_i +K \sin \o \cos \p'_i 
\end{equation}
that can now be solved using a linear least-square techniques with unknowns: $\g$, $K \cos \o$ and $K \sin \o$. The values of $K$ and $T_0$ follows from the latter determination. As in the eccentric case, LOSP offers in addition the option to scan a given period interval to search for the best period values.

\subsection{SB2 circular orbit} \label{ssect: sb2ne}

SB2 circular orbits are handle pretty much the same way as eccentric orbits, though the equation system is this time:
\begin{equation} \left\{ \begin{array}{l}
v_1=\g_1+K_1\sin \p\\
v_2=\g_2+K_2\sin (\p+\pi)
\end{array}\right. \end{equation}
which, for a given value of the period, has five unknowns: $\g_1$,  $\g_2$, $K_1$, $K_2$ and $T_0$.  It is easy to show that µ
\begin{equation}
 \frac{v_1(\p)-\gamma_1}{K_1}= -\frac{v_2(\p)-\gamma_2}{K_2} 
\end{equation}
or equivalently 
\begin{equation} 
v_2(\p) =  b + c v_1(\p)
\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 the same orthogonal (because both $v_1$ and $v_2$ are spoiled by errors) regression technique as in the SB2 eccentric case. On can apply the following transformation in the radial velocity space :
\begin{equation} \left\{ \begin{array}{l}
v^*_1=\sqrt{-c}\left(v_1-\frac{b}{1-c} \right)\\
v^*_2=\frac{1}{\sqrt{-c}}\left(v_2-\frac{b}{1-c} \right)
\end{array}\right. \end{equation}
to create the following equation system:
\begin{equation} \left\{ \begin{array}{l}
v^*_1= \G + \mathcal{K} \cos \o +  \mathcal{K} \sin \o \\
v^*_2= - \G - \mathcal{K} \cos \o -  \mathcal{K} \sin \o
\end{array}\right. \end{equation}
with $\G=\frac{\sqrt{-c}}{1-c}\left(\g_1-\g_2\right)$  and
$\mathcal{K}=\sqrt{K_1K_2}$. The latter system is now linear in the
parameters $\G$,  $\mathcal{K} \cos \o$ and $ \mathcal{K} \sin \o$ and
can be solved using a classical linear least square. The orbital
solution of the SB2 system is then easily found from the best-fit
values $\G$,  $\mathcal{K} \cos \o$ and $ \mathcal{K} \sin \o$ and from the results of the orthogonal regression that yields $b$ and $c$. 

\section{Where to get the package ?}

The  LOSP package consists of a FORTRAN-77 coded set of
programs/routines  handled transparently (from the user point of view)
through a single BASH script. As a last capability, phase diagrams of
the data and of the best-fit RV-curves are displayed using SUPER
MONGO. This allows a visual inspection of the adjustment results. LOSP
is quick, light ($\sim200$~Ko), easy to install and to use. It is
available to the community and can be downloaded as a single archive
file at ftp://arachnos.astro.ulg.ac.be/usr/sana/LOSP ($\sim30$~Ko).
More detailed documentation is also available at the same address ($\sim100$~Ko).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Summary}


%_____________ ACKNOWLEDGMENTS __________________________________________


\section*{Acknowledgments}


%_____________ BIBLIOGRAPHY _____________________________________________

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



\appendix




\section{Transformation in the RV space: circular orbit}

% 1
\begin{equation} \left\{ \begin{array}{lll}
 v_1&=&\g_1+K_1\cos\left(\t+\phi\right)  \\
 v_2&=&\g_2-K_2\cos\left(\t+\phi\right)  \end{array} \right.
\end{equation}

% 2
\begin{equation} \left\{ \begin{array}{lll}
 \left(v_1-\frac{b}{1-c}\right)\sqrt{-c}&=&\left(\g_1-\frac{\g_2-c\g_1}{1-c}\right)\sqrt{-c}+K_1\sqrt{-c}\cos\left(\t+\phi\right)  \\
 \left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}&=&\left(\g_2-\frac{\g_2-c\g_1}{1-c}\right)\frac{1}{\sqrt{-c}}-K_2\frac{1}{\sqrt{-c}}\cos\left(\t+\phi\right)  \end{array} \right.
\end{equation}

% 3
\begin{equation} \left\{ \begin{array}{lll}
 \left(v_1-\frac{b}{1-c}\right)\sqrt{-c}&=&\frac{\g_1-c\g_1-\g_2+c\g_1}{1-c}\sqrt{-c}+\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \\
 \left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}&=&\frac{\g_2-c\g_2-\g_2+c\g_1}{1-c}\frac{1}{\sqrt{-c}}-\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \end{array} \right.
\end{equation}

% 4
\begin{equation} \left\{ \begin{array}{lll}
 \left(v_1-\frac{b}{1-c}\right)\sqrt{-c}&=&\frac{\g_1-\g_2}{1-c}\left(\frac{1}{\sqrt{-c}}\right)^{-1}+\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \\
 \left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}&=&c\frac{\g_1-\g_2}{1-c}\frac{1}{\sqrt{-c}}-\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \end{array} \right.
\end{equation}

% 5
\begin{equation} \left\{ \begin{array}{lll}
 \left(v_1-\frac{b}{1-c}\right)\sqrt{-c}&=&\left(\g_1-\g_2\right)\frac{\sqrt{-c}}{1-c}+\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \\
 \left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}&=&-\left(\g_1-\g_2\right)\frac{\sqrt{-c}}{1-c}-\sqrt{K_1K_2}\cos\left(\t+\phi\right)  \end{array} \right.
\end{equation}

Linearizing the system:

\begin{equation} \left\{ \begin{array}{lllll}
 v_1^{*}&=&\left(v_1-\frac{b}{1-c}\right)\sqrt{-c}&=&\left(\g_1-\g_2\right)\frac{\sqrt{-c}}{1-c}+\sqrt{K_1K_2}\cos\t\sin\phi+\sqrt{K_1K_2}\sin\t\cos\phi \\
 v_2^{*}&=&\left(v_2-\frac{b}{1-c}\right)\frac{1}{\sqrt{-c}}&=&-\left(\g_1-\g_2\right)\frac{\sqrt{-c}}{1-c}-\sqrt{K_1K_2}\cos\t\sin\phi-\sqrt{K_1K_2}\sin\t\cos\phi
 \end{array} \right.
\end{equation}

\begin{eqnarray} 
\left( \begin{array}{c} | \\ v_1^* \\ | \\ v_2^* \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{ccc}
| &   |      & | \\
\frac{\sqrt{-c}}{1-c} & \sin\phi &\cos\phi \\ 
| &   |      & | \\
-\frac{\sqrt{-c}}{1-c} & -\sin\phi & -\cos\phi \\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \g_1-\g_2 \\ \sqrt{K_1K_2}\cos\t \\ \sqrt{K_1K_2}\sin\t \\\end{array} \right)
\end{eqnarray} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\section{Error propagation: SB2 Circular orbit}

\begin{itemize}
 \item $\g_1$ and $\g_2$ \\
    Definitions:
   \begin{center}
   \begin{eqnarray}
    \Gamma &=& \gamma_1 - \gamma_2 \\
    C &=& \sqrt{K_{1}K_{2}}\cos{\t} \\
    S&=&\sqrt{K_1K_2}\sin{\t} \\
    c&=&-\frac{K_2}{K_1} \\
    b&=&\gamma_2 - a\gamma_1
   \end{eqnarray}
   \end{center}

   \begin{center}
   \begin{equation}
    \g_1=\frac{\Gamma+b}{1-c}
   \end{equation}
   \end{center}
   \begin{center}
   \begin{eqnarray}
    \frac{\partial\g_1}{\partial\Gamma} &=& \frac{1}{1-c} \\
    \frac{\partial\g_1}{\partial b} &=& \frac{1}{1-c} \\
    \frac{\partial\g_1}{\partial c} &=& -\frac{\g_1}{1-c} \\
    \sb_{\g_1} &=& \frac{1}{\left(1-c\right)^2}\sb_\Gamma + \frac{1}{\left(1-c\right)^2}\sb_b + \frac{\left(\Gamma+b\right)^2}{\left(1-c\right)^4}\sb_c
   \end{eqnarray}
   \end{center}

   \begin{equation}
    \g_2=\frac{b+c\Gamma}{\left(1-c\right)}
   \end{equation}
   \begin{eqnarray}
    \frac{\partial\g_2}{\partial b} &=& \frac{1}{1-c} \\
    \frac{\partial\g_2}{\partial\Gamma} &=& \frac{c}{1-c} \\
    \frac{\partial\g_2}{\partial c} &=& \frac{\g_1}{\left(1-c\right)} \\
    \sb_{\g_2} &=& \frac{c^2}{\left(1-c\right)^2}\sb_\Gamma + \frac{1}{\left(1-c\right)^2}\sb_{\g_2} + \frac{\g^2_1}{\left(1-c\right)^2}\sb_c
   \end{eqnarray}

 \vspace{.5cm}
 \item $K_1$ and $K_2$\\

  \begin{equation}
   K_1 = \sqrt{\frac{C^2+S^2}{-c}}
  \end{equation}
  \begin{eqnarray}
   \frac{\partial K_1}{\partial C} &=& \frac{C}{S^2+C^2}K_1 \\
   \frac{\partial K_1}{\partial S} &=& \frac{S}{S^2+C^2}K_1 \\
   \frac{\partial K_1}{\partial c} &=& \frac{-1}{2c}K_1 \\
   \sb_{K_1} &=& \frac{C^2K_1^2}{\left(S^2+C^2\right)^2}\sb_C + \frac{S^2K_1^2}{\left(S^2+C^2\right)^2}\sb_S + 2\frac{SC}{\left(S^2+C^2\right)^2}K_1^2\sb_{CS} + \frac{K_1^2}{4a^2}\sb_c 
  \end{eqnarray}   

  \begin{equation}
   K_2=\sqrt{-c\left(S^2+C^2\right)}
  \end{equation}
  \begin{eqnarray}
   \frac{\partial K_2}{\partial S} &=& \frac{S}{S^2+C^2}K_2 \\
   \frac{\partial K_2}{\partial C} &=& \frac{C}{S^2+C^2}K_2 \\
   \frac{\partial K_2}{\partial c} &=& \frac{1}{-2c}\sqrt{-c\left(S^2+C^2\right)} \\
   \sb_{K_2} &=& \frac{K_2^2}{\left(S^2+C^2\right)^2}\left(S^2\sb_S+C^2\sb_C+2CS\sb_{CS}\right) + \frac{K_2^2}{4c^2}\sb_c
  \end{eqnarray}

 \vspace{.5cm}
 \item $\t=\arctan{\left(\frac{S}{C}\right)}$

  \begin{eqnarray}
   \frac{\partial\t}{\partial S} &=& \frac{C}{S^2+C^2} \\
   \frac{\partial\t}{\partial C} &=& -\frac{S}{S^2+C^2} \\
   \sb_\t &=& \frac{C^2}{\left(S^2+C^2\right)^2}\sb_S + \frac{S^2}{\left(S^2+C^2\right)^2}\sb_C - 2\frac{CS}{\left(S^2+C^2\right)^2}\sb_{SC}
  \end{eqnarray}

 \vspace{0.5cm}
 \item $T_0=t_{arb}-\frac{P}{2\pi}\arctan{\frac{S}{C}}; \qquad \sb_{T_0}=\frac{P}{2\pi}\s_\t$

 \vspace{0.5cm}
 \item $a_1\sin{i}=\frac{K_1}{2\pi}P(\textrm{sec}); \qquad \s_{a_1\sin{i}}=\frac{P}{2\pi}\s_{K_1}$

 \vspace{0.5cm}
 \item $a_2\sin{i}=\frac{K_2}{2\pi}P(\textrm{sec}); \qquad \s_{a_2\sin{i}}=\frac{P}{2\pi}\s_{K_2}$

 \vspace{0.5cm}
 \item $m_1\sin^3{i}=\frac{P}{2\pi G}\left(S^2+C^2\right)^{3/2}\left(v-1\right)^2\left(-c\right)^{-1/2}$

   \begin{eqnarray}
    \frac{\partial m_1\sin^3i}{\partial C} &=& 3\frac{m_1\sin^3i}{S^2+C^2}C \\
    \frac{\partial m_1\sin^3i}{\partial S} &=& 3\frac{m_1\sin^3i}{S^2+C^2}S \\
    \frac{\partial m_1\sin^3i}{\partial c} &=& \frac{m_1sin^3i}{2c\left(c-1\right)}\left(3c+1\right)
   \end{eqnarray}

  \vspace{5mm}
  \item $m_2\sin^3i=\frac{P}{2\pi G}\left(1-c\right)^2\left(\frac{C^2+S^2}{-c}\right)^{3/2}$

   \begin{eqnarray}
    \frac{\partial m_2\sin^3i}{\partial C} &=& 3\frac{m_2sin^3i}{S^2+C^2}C \\
    \frac{\partial m_2\sin^3i}{\partial S} &=& 3\frac{m_2sin^3i}{S^2+C^2}S \\
    \frac{\partial m_2\sin^3i}{\partial c} &=& -\frac{m_2sin^3i}{2c\left(1-c\right)}\left(c+3\right)
   \end{eqnarray}

\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 : read the input data file,
\item ORBIT : computes the SB1 or 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 only, performs the orthogonal regression and builds the {\it fake} SB1 system,
\item INIT : builds the linear system to be inverted ...
\item APPROXEL : computes the approximate orbital solution 
\item DC : given that the DC criterion is encountered, performed the DC recursively until sufficient convergence is reached, compute the improved orbital element of the handled SB1 system,
\item FINALEL : for SB2 systems, derived the final SB2 parameters from the {\it fake} SB1 final solution; for SB1 system, computes the mass function.\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)\\
$isw3$ & Is it an SB1 (isw3=1) or SB2 (isw3=2) system ? \\
$isw4$ & Number of indirect (or additional) parameters to be considered in the fit ?\\
$isw5$ & Is it an circular (isw5=0) or eccentric orbit (isw5=1) ? \\
$isw6$ & Operation to be considered : 0: fit the orbit  \\
$mlabel$ & Maximum number of D.C. steps \\
$m$ & Number of coefficient in the harmonic analysis \\
$id$ & object name \\
$P$ & initial guess of 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)_{i=1}^{np}$ & Julian date, $v_1$ and weight , $v_2$ and weight \\
\end{tabular}

\subsection{The orthogonal regression: MASSRATIO}
The routine is performed for SB2 system only (isw3=2). 
It selects all the points for which wtp(i)=wts(i) in the data lists and passes them to the orthogonal regression routine. The points for which wtp(i)$\ne$wts(i) are not considered in the orthogonal regression (QUESTION: Should we also reject all points with a null weight?).  Let us note in this section $w_i=wtp(i)=wts(i)$. The weight are first normalized so that $\sum_{i=1}^{np}w_i=np$. 

\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}{np}c^2\mathcal{F}\\
\sb_b=\frac{\chi^2_\nu}{np}(s^2+c^2+c^2\x^2\mathcal{F}\\
\chi^2_\nu=\frac{N\chi^2}{(np-2)\sum_i w_i} \\
\chi^2=\sum_i w_i \frac{\left(y_i-cx_i-b\right)^2}{s^2+c^2}
\end{eqnarray}

The code also compute 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(i) \ne wts(i)$
We then used the best-fit relation to bring the weight to an absolute scale:
\begin{eqnarray}
 \sb&\approx&\frac{\chi^2_\nu}{\sum_i w_i/np}=\frac{np}{\nu \sum_i w_i} \sum_i w_i\frac{(y_i-cx_i-b)^2}{s^2+c^2}\\
 wtp(i)&=&\frac{w_i}{\sum_k w_k /np}\frac{1}{\sb} \\
 wts(i)&=&wtp_i/s^2
\end{eqnarray}
so that $\sb_{v_p(i)}=1/wtp(i)$ and $\sb_{v_s(i)}=1/wts(i)$

\subsubsection{Initialization of the fake SB1 system}
Let us write $n = 2np$. For $i=1...n/2$
\begin{eqnarray}
v^*(i)&=&\left(v_1(i)-\frac{b}{1-c}\right)\sqrt{c} \\
\left(w^*(i)\right)^{-1}&=&-c\sb_{v_p(i)}-\frac{c}{1-c^2}\sb_b+\left(\frac{v^*(i)}{2c}+\frac{b\sqrt{-c}}{\left(1-c\right)^2}\right)^2\sb_c
\end{eqnarray}
 For $i=(n/2+1)...n$
\begin{eqnarray}
v^*(i)&=&\left(v_2(k)-\frac{b}{1-c}\right)\frac{1}{\sqrt{c}} \hspace*{0.5cm}\mathrm{with }\hspace*{0.5cm} k=i-n/2\\
\left(w^*(i)\right)^{-1}&=&-\frac{\sb_{v_p(i)}}{c}-\frac{\sb_b}{c(1-c^2)}+\left(\frac{v^*(i)}{2c}+\frac{b}{\sqrt{-c}\left(1-c\right)^2}\right)^2\sb_c
\end{eqnarray}
QUESTION: Should be account for $\s_c$ and $\s_b$ while converting $V_1$ and $v_2$ into $v^*$. If not, then the previous equations can be rewritten:
\begin{eqnarray}
\left(w^*(i)\right)^{-1}&=&-c\sb_{v_p(i)}-\frac{c}{1-c^2}\sb_b
\end{eqnarray}
 For $i=(n/2+1)...n$
\begin{eqnarray}
\left(w^*(i)\right)^{-1}&=&-\frac{\sb_{v_p(i)}}{c}
\end{eqnarray}
and no absolute scale (Sect.\ref{ssect: absscale}) need to be defined. 

\subsection{Initializing the Wilsing-Russell method: INIT}
The weights are renormalized, accounting for all the (transformed) RVs. Then, from the $\{v^*(i),w^*(i)\}$, $i=1...n$, one constructs the matrix of the linear system 

\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)
\end{eqnarray} 
for the {\it fake} SB1 system were the separation between the upper and lower part of the matrix happen between n=np and n=np+1. For an SB1 binary, the system is:

\begin{eqnarray} 
\left( \begin{array}{c} | \\ v \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{ccc}
| &   |      & | \\
1 & \sin kM' &\cos kM' \\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \Gamma \\ S_k \\ C_k \\\end{array} \right)
\end{eqnarray} 

In either cases, the system  is solved using the Numerical Recipes routine GAUSSJ. This yields the bets-fit parameters $\Gamma$, $S^*_k$ and $C^*_k$ for the {\it fake} SB1 system or  $v_0$, $S_k$ and $C_k$ for a real SB1 binary.

Then computes the QUALITY of the least-square fit :
\begin{equation}\chi^2=\sum_{i=1}^{np} wt(i) \left( v^*_i - \Gamma - \sum_{k=1}^{k=m} \sin kM' S^*_k -\sum_{k=1}^{k=m} \cos kM' C^*_k \right)^2 + \sum_{i=np+1}^{2np} wt(i) \left( v^*_i + \Gamma + \sum_{k=1}^{k=m} \sin kM' S^*_k + \sum_{k=1}^{k=m} \cos kM' C^*_k \right)^2  
\end{equation}
for the {\it fake} SB1 system, while
\begin{equation}\chi^2=\sum_{i=1}^{N} wt(i) \left( v_i - \Gamma - \sum_{k=1}^{k=m} \sin kM' S_k -\sum_{k=1}^{k=m} \cos kM' C_k \right)^2 
\end{equation}
for a real SB1 system.
Hence 
\begin{eqnarray}
\chi^2_\nu&=&\frac{\chi^2}{N-(2m+1)} \label{eq: quality}  \\
r1^2=\frac{\chi^2}{(\sum_i w_i) (1-\frac{2m+1}{N}}
\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<5E-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 < 2.5E-13$) following exactly WHS67. Finally compute 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$ of the fit (Eq. \ref{eq: quality}):

\begin{equation}
A^2_j=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 \\
                                   \arctan(-S_j/C_j)+\pi  & if & C_j < 0 \\
                                   \frac{\pi}{2}  & if & C_j = 0 \& S_j \le 0 \\
                                   \frac{3\pi}{2}  & if & C_j = 0 \& 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}
\sb_{\a_j}=\frac{C_j^2*\sb_{S_j}+S_j^2*\sb_{C_j}}{A^4_j}
a \sin i &=& \frac{A_1}{b_1\mu}\\
T_0 &=& hjdl-\frac{\a_2-\a_1+\b_1-\b_2}{\mu}  \\
\mathcal{K} = \mu*(a \sin i)/dsqrt{1-e^2} \\
\sb_e=\frac{e^2*\sb_{A_1}+\sb_{A_2}}{A^2_1} \\
\sb_\o = 4\sb_{\a_1}+\sb_{\a_2} \\
\sb_\mathcal{K}=\sb_{A_1} \\
\sb_{a \sin i}=\sb_\mathcal{K}/\mu^2 \\
\sb_{T_0} = \frac{\sb_{\a_1}+\sb_{\a_2}}{\mu}
\end{eqnarray}

\subsection{The differential correction: DC}
When the following criterion is matched, the program enters apply a differential correction (D.C) phase as explained is Sect.~ 
\begin{eqnarray}e>2\s_e &\mathrm{and} & e> 0.03
\end{eqnarray}
If the option 'forced=1' has been chosen, then the D.C. is applied whatever the result of the above criterion.
If chosen the D.C. also allows to improve the adopted value for the period.
The routine DC thus builds the appropriate linear equation system.
First the true anomaly $\t_i$ corresponding to the approximate solution is computed for each observation and the $(O-C)_i=v^*_{\mathrm{obs},i}-v_\mathrm{th}(\t_i)$. Then the system 



\begin{eqnarray} 
\left( \begin{array}{c} | \\ (o_i-c_i)_1^* \\ | \\ | \\ (o_i-c_i)_2^* \\ | \\ \end{array} \right) & = & 
\left( \begin{array}{cccccc}
| &   |      & | \\
1 & e\cos\o+cos(\o+\t_i) & \mathcal{K}(\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} & -\mathcal{K} (e \sin\o+ \sin(\o+\t_i) & \mathcal{K}\mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} & \mathcal{K} (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}}\\ 
| &   |      & | \\
| &   |      & | \\
-1 & -e\cos\o+cos(\o+\t_i) & -\mathcal{K}(\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} & \mathcal{K} (e \sin\o+ \sin(\o+\t_i) & -\mathcal{K}\mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} & -\mathcal{K} (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}}\\ 
| &   |      & | \\
\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} 
for the {\it fake} SB1 system were the separation between the upper and lower part of the matrix happen between n=np and n=np+1. For an SB1 binary, the system is:

\begin{eqnarray} 
\left( \begin{array}{c} | \\ (o_i-c_i) \\ | \\  \end{array} \right) & = & 
\left( \begin{array}{cccccc}
| &   |      & | \\
1 & e\cos\o+cos(\o+\t_i) & K (\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} & -K (e \sin\o+ \sin(\o+\t_i) & K \mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} & K (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}}\\ 
| &   |      & | \\
\end{array} \right)
\left(\begin{array}{c} \Delta\gamma \\ \Delta K \\ \Delta e  \\ \Delta \o  \\ \Delta T_0  \\ \Delta \mu  \\  \end{array} \right)
\end{eqnarray} 
The matrix is inverted through the GAUSSJ routines and the correction are found.  The quality of the fit is estimated by the QUALITY2 routine which computed the (O-C) as well as 
\begin{eqnarray}
\chi^2&=&\sum_{i=1}^{i=np} w_i \left((O-C)_i- \Delta\Gamma -\Delta\mathcal{K} e\cos\o+cos(\o+\t_i) - \Delta e\mathcal{K}(\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} + \Delta \o \mathcal{K} (e \sin\o+ \sin(\o+\t_i) - \Delta T_0 \mathcal{K}\mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} - \Delta \mu \mathcal{K} (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} \right)^2   \\
&=&+\sum_{i=np+1}^{i=2np} w_i \left((O-C)_i+ \Delta\Gamma +\Delta\mathcal{K} e\cos\o+cos(\o+\t_i) + \Delta e\mathcal{K}(\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} - \Delta \o \mathcal{K} (e \sin\o+ \sin(\o+\t_i) + \Delta T_0 \mathcal{K}\mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} + \Delta \mu \mathcal{K} (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} \right)^2  
\end{eqnarray}
for the fake SB1 system, and\begin{eqnarray}
\chi^2&=&\sum_{i=1}^{i=N } w_i \left((O-C)_i- \Delta\gamma -\Delta         K  e\cos\o+cos(\o+\t_i) - \Delta e K(\cos\o-\sin\t_i \sin(\o+\t_i)\frac{2+e\cos\t_i}{1-e^2} + \Delta \o K  (e \sin\o+ \sin(\o+\t_i) - \Delta T_0 K \mu (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} - \Delta \mu K  (1+e\cos\t_i)^2\frac{\sin(\o+\t_i}{(1-e^2)^{3/2}} \right)^2   \end{eqnarray} 
for a real SB1 system. Then :
\begin{eqnarray}
\chi^2_\nu=\frac{\chi^2}{\nu}\\
r_1^2=\frac{N\chi^2}{\sum_i wt_i \nu}
\end{eqnarray}
where $\nu=N-6$ if the period (actually, $\mu=2\pi/P$) is included in the DC; $\nu=N-5$ otherwise. 

The improved parameters $\Gamma$, $\mathcal{K}$ (resp. $\gamma$ and $K$ for a real SB1 system), $e$, $\o$, $T_0$ and eventually $\mu=2\pi/P$. 
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 allowed iteration $mlabel$ has not been reached and as long as the convergence criterion is not matched. The latter requires that
\begin{equation}
\Delta^2_\mathcal{K}>\sb_\mathcal{K}/3
\end{equation}
where $\Delta_\mathcal{K}$ is the D.C. step along the $\mathcal{K}$ direction while $\s_\mathcal{K}$ is the uncertainty on the latter parameter.


\subsection{Computing the final elements: FINALEL}

 The FINALEL routine then 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. 


\subsection{Estimating the quality of the final SB2 system: MODELE}
The routine MODELE compute the O-C and the $\chi^2$ with respect to the final solution. Hence, for an SB2 system,
\begin{equation}
sv2=\sum_{i=1}^{np}wtp(i)\left(vp(i)-\gamma_1-K_1\left(e\cos\o+\cos(\o+\t_i)\right)\right)^2 + \sum_{i=1}^{np}wts(i)\left(vs(i)-\gamma_2+K_2\left(e\cos\o+\cos(\o+\t_i)\right)\right)^2 
\end{equation}
\begin{equation}
r_1^2=\frac{sv2*N}{\nu \sum_i(wtp_i+wts_i) }
\end{equation}
where $\nu=N-6-isw4$ if the period has been adjusted; $\nu=N-5-isw4$ otherwise.
For a real SB1 system, one obtains: 
\begin{eqnarray}
sv2=\sum_{i=1}{N}wt(i)\left(v(i)-\gamma-K\left(e\cos\o+\cos(\o+\t_i)\right)\right)^2 \\
r_1^2=\frac{sv2*N}{\nu \sum_i w_i }
\end{eqnarray}
where $\nu$ has the same definition as previously.

\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}
