Differential Correction of Orbits

In this section the correction of an orbit for a visual double star by means of a least squares fit using the method of differential correction will be explained in detail. The process of calculation follows the formulae and advises given in ?Double Stars? (Heintz, 1967)4.

First of all, a misunderstanding should be avoided. Many people think, using ?differential corrections? means that only small corrections can be applied and that the start orbit must be very close to the final solution. This is not correct at all.
When handling with care, quite large corrections can result or are possible, and the method is very powerful because only a limited number of orbits has to be calculated in order to get the final solution. During the approximation, the calculator can look at residuals and control the process of calculation. The procedure cannot be automized in many cases, especially if there are few measures or the measures are distributed over a part of the ellipse only.

In addition, the final solution of an orbit calculation does not depend so much on the method of calculation. It depends on the weights assigned to the observations. There was and is a lot of discussion concerning the weight to be given to visual measurements compared with photographic and speckle measurements.
Only experience and comparison of measurements with definitive or very well determined orbits or with pairs in very slow motion can allow to determine reliable weights.
The mathematical rule for the calculation of a weight W is:



where ? is the error of the measurement.

Now, the formulae for the corrections have to be set up. We follow the principle, that a change d?j of some coordinate ? is composed of elements corrections d?i given by the formula:



where the d?j are the residuals and the d?i are the improvements to be found. The quotients are calculated from initial elements found e.g. by the graphical method described in the previous chapter.

Heintz has given the following formulae when using polar coordinates. So, all differentials are expressed in units of degrees:

e = sin??
d? = 57.296 sec ? de
d?  = ? sec2 ? dT
dař = 57.296 da/a
dm = N sec2 ? d?

The assumption is that all observations have been collected (corrected for precession) and a start orbit is available.

For each observational position (it may consist of the mean of several observations because of economy of calculation and is called a ?normal position?) one has to calculate 2 Equations of Condition:

d? + B di + C d? + F d? + G dm + H d? = d?
d?? + b di + c d? + f d? + g dm + h d? = 57.296 d???

where d? and d? are the residuals for this observational position, i.e. d??= ? (observed) ? ? (calculated, using the start orbit) and d??= ? (observed) ? ? (calculated, using the start orbit).

The auxiliary functions are defined as follows:

? = sec ? sin ? (2 + e cos ?)
? = sec ? (1 + e cos ?)2
? = - sec ? cos ? (1 + e cos ?)
? = tan ? sin ? (1 + e cos ?)

and the coefficients B, C, F, G, H and b, c, f, g and h become:

B = - cos2 (? - ?) tan (? + ?) sin i
C = + cos2 (? - ?) sec2 (? + ?) cos i
F = - ? C
G = + ? C t
H = + ? C
b = - sin2 (? - ?) tan i
c = - sin (? - ?) cos (? - ?) sin i tan i
f = - (? c + ?)
g = + (? c + ?) ?
h = + ? c + ?

Each equation of condition has to be multiplied with the square root of the weight which is assigned to this measurement (left side and right side!). For determination of the weights: see the suggestion by Heintz in ?Double Stars?4, page 46 or Hartkopf et al: The weighting game in ?Binary star Orbits from Speckle Interferometry II?.7 For example, visual angle measurements obtain more weight than visual distance measurements.
Let?s assume that there are 50 normal places, i.e. 50 times 2 Equations of Condition and let us continue with the equations for d?.
There is a system of 50 linear equations with the 6 unknowns d?, di, d?, ?:
On the left side there is a 50 x 6 ? Matrix (let?s call it A), on the right side a 
50 x 1 ? vector (let?s call it R1) containing the residuals.

Now the system of Normal Equations has to be set up by Matrix ? multiplication:

(AT *  A) * ca = AT * R1 

where  AT  is the transposed Matrix of A and ca is the 6 x 1 - vector containing the unknown corrections:
	
	d?
 	di
ca =	d?
	d?
	dm
	d?

On the left there is a 6 x 6 ? Matrix  AT *  A = a, on the right a 6 x 1 ? Vector 
AT * R1 = ra. 
This is the system of Normal Equations.

Solving such a system of equations is pre-programmed in most mathematical programs and gives the solutions for d?, di, d?, d?? dm, d?. But the procedure works only if the observed arc represents the orbit fairly well and when there are many good measurements. Otherwise, large and unreliable corrections will result.

Now the corrections have to be applied to the elements of the start orbit; this yields the new elements.

The same can be done for the system of linear equations for d? to get the correction for ař and a, respectively. This gives us a 6 x 6 ? Matrix  BT *  B = b and a 6 x 1 ? Vector rb. Differential corrections with the Matrix b work well only if there are many good measurements for the distance like for example in the case of a binary such as McA14Aa. (Indeed,  McA14Aa can be corrected in rectangular coordinates also very well).

Is the corrected orbit better? This has to be checked for every new set of elements.
For every new orbit one calculates the 2 sums of squared errors. We call it SUM:

Wi = weight assigned to the ith ? equation of condition, n = number of normal places.
2 SUMs have to be calculated, one for the angles, the second one for the distances.

Important note: 
It is very much recommended to add up the two 6 x 6 normal equations and residual vectors for ? and ??in order to get a 7 x 7 matrix c and a 7 x 1 residual vector rc for correction of all 7 elements simultaneously. For adding up the matrixes and the vectors to get the 7 x 7 matrix c and the 7 x 1 vector rc, add zeros in this way:



and



Operating with c and rc, the algorithm is much more robust and unreliable excessive corrections happen less frequently. 

In case the observed arc is short one can also delete individual columns from the normal equations in order to reduce the number of corrected orbital elements in one step. Some calculators for example vary the period P by a step-width of a certain amount and look for the corrections of the remaining elements. Normally, the whole procedure is an iterative one, and it may take about 5 to up to about 20 to 50 steps, depending on the observational material and the observed arc.

For checking the reliability of the programs, first generate two artificial orbits X and Y with somewhat different orbital elements. From orbit X ?observation points? are calculated and orbit Y serves as start orbit. If the program works properly the corrected orbit must become orbit X. The observation points (about 15 to 20) should be well distributed around the complete ellipse.
For example:

X ? orbit: P=360, T=1903, a=0.70, e=0.30, i=131.5, omega=149.0, node = 120.0
Y ? orbit: P=340, T=1900, a=0.90, e=0.35, i=140.0, omega=140.0, node = 117.0

Although the two sets of elements differ considerably, the method of differential corrections will find the correct solution in 2 steps only.
In practice, for real observations one will see: the less complete the observational material and the less accurate the observations, the more difficult the procedure.
In case of a too undetermined arc (little curvature, periastron not observed or too far in the future etc.), the corrections may become very large and unreliable as mentioned above.

Possible ways to continue in such a case:

1. Correct only 5 or 4 or even only 3 elements simultaneously.
2. Do not try to correct at all, wait for new measurements.

Some calculators do it in the following way: They vary P, T and e in a 3 dimensional grid (fixed values) and correct the remaining elements thus calculating hundreds or even  thousands of orbits and as many values for the 2 SUMs of residuals. 

Calculations of errors of orbital elements can be performed if the observed arc contains the 2 ansae and hence defines the orbit already well:
Calculate the inverse of the 7 x 7 matrix c, this gives c-1. 
The diagonal elements are c-111, c-122, ?c-177. 
Calculate the sum of the residuals:



Now, the errors for the orbital elements can be calculated:



























1