#### Useful Links

## 1. Introduction

Lots of engineering problems do not require advanced math or physics. However, they may be very complicated since they contain many different objects and links between them. This article is devoted to one such problem. All code of this article is compatible with Visual C# 2010 Express.

## 2. Background

Determination of orbits of satellites, in fact, is a statistical processing of measurements. As well as education of former centuries required knowledge of Greek and Latin languages, any specialist of determination of orbits should know the General linear method and Kalman filter. This article contains a profound description of the application of the General linear method and a brief description of the technology of Kalman filter.

## 3. External Libraries

Any advanced software requires interoperability with third party software. Two additional libraries are used in this article. The first one calculates Gravity of Earth. The second one calculates parameters of Atmosphere of Earth.The adapter pattern is used for software interoperability. Let us show usage of this pattern. The following code contains declaration of function which calculates Gravity of Earth:

private void Forces(int N0, int NK, double X, double Y, double Z, out double FX,
out double FY, out double FZ)

This function is not known by the framework. However framework contains `IObjectTransformer`

interface.

public interface IObjectTransformer
{
string[] Input
{
get;
}
string[] Output
{
get;
}
object GetInputType(int i);
object GetOutputType(int i);
void Calculate(object[] input, object[] output);
}

This interface is an abstract transformation operation. Members of this interface have the following meaning:

| Member name | Comment |

1 | `Input` | Names of input variables |

2 | `Output` | Names of output variables |

3 | `GetInputType(int i)` | Type of i-th input variable |

4 | `GetOutputType(int i)` | Type of i-th output variable |

5 | `Calculate(object[] input, object[] output)` | Calculation |

The following code snippet shows implementation of this interface:

public class Gravity : IObjectTransformer, ICategoryObject, IPropertiesEditor,
IObjectFactory
{
static private readonly string[] inps = new string[] { "x", "y", "z" };
static private readonly string[] outs = new string[] { "Gx", "Gy", "Gz" };
string[] IObjectTransformer.Input
{
get { return inps; }
}
string[] IObjectTransformer.Output
{
get { return outs; }
}
object IObjectTransformer.GetInputType(int i)
{
return ret;
}
object IObjectTransformer.GetOutputType(int i)
{
return ret;
}
void IObjectTransformer.Calculate(object[] input, object[] output)
{
double x = (double)input[0];
double y = (double)input[1];
double z = (double)input[2];
double fx, fy, fz;
Forces(n0, nk, x, y, z, out fx, out fy, out fz);
output[0] = fx;
output[1] = fy;
output[2] = fz;
}
}

Let us explain this code. Input variables are named by "`x`

", "`y`

" and "`z`

" respectively. These names are used for interoperability. Similarly output variables are named "`Gx`

", "`Gy`

","`Gz`

". So we have three input and three output variables. All these variables have the "`Double`

" type. System of types used by the framework does not coincide with the .NET one, because the framework should distinguish arrays with different length. The `ArrayReturnType`

is used for arrays. However for `double `

variable type, the `const Double ret = 0;`

is used. The `Calculate`

calls `Forces`

function. So `Forces`

function becomes implicitly visible by the framework since it is called by method of `IObjectTransformer`

interface. Now let us explain how the framework can find necessary objects. The following interface is used for this purpose:

public interface IObjectFactory
{
string[] Names
{
get;
}
ICategoryObject this[string name]
{
get;
}
}

The `Names`

property returns names of objects which provides the library. The `ICategoryObject this[string name]`

returns object by name. The following code presents implementation of this interface:

string[] IObjectFactory.Names
{
get { return new string[] { "Gravity 36 * 36" }; }
}
ICategoryObject IObjectFactory.this[string name]
{
get { return this; }
}

The "Gravity 36 * 36" is indicated in user interface by the following way. First of all, we will set external object.

If we open its editor of properties and select library (*Grav36.dll*), then the list of accessible objects will appear in the `combobox `

below:

Then the necessary object ("Gravity 36 * 36") should be selected. After this operation, the object shall have a new icon:

The following code represents searching of factory of objects:

static public IObjectFactory GetFactory(Assembly ass)
{
Type[] types = ass.GetTypes(); foreach (Type type in types)
{
if (type.GetInterface("IObjectFactory") != null)
{
FieldInfo fi = type.GetField("Object");
if (fi != null)
{
return fi.GetValue("Object") as IObjectFactory;
}
else
{
return type.GetConstructor(new Type[] { }).Invoke(new object[]
{ }) as IObjectFactory;
}
}
}
return null;
}

This code defines such type of assembly which implements `IObjectFactory`

interface. If this type contains singleton (named "`Object`

"), then this singleton is used. Otherwise default constructor of object is used.

## 4. Containers

Complicated applications require a lot of squares of objects and arrows. It is very difficult to understand such a picture. However, you can merge a set of squares of objects and arrows. You can also combine sets of squares into containers. Then, the result picture shall be clear and easy to understand. Determination of orbits requires equations of motions that require gravity and atmosphere models.

In this picture, **Gravity** and **Atmosphere** use dynamic link libraries as it had been described in 3. It is convenient to merge these components to one container. For the creation of the container, the user should create a scenario and select the *Tools/Container designer* option of the main menu. After selecting the *Tools/Container designer* option, we have:

The right part of this form enables us to select public components by checking checkboxes. The left one enables us to set geometrical positions of public components. When these positions are set, you can save the container to the **.cont* file. In this example, the public components are **Atmosphere**, **Gravity** and **Motion Equations**. Then, we should create a palette button. To do this, you should place the **.cont* to the *Application path/cont* directory, and create an icon of the button and place it in the application path. Then, you should create or edit the *Containers.xml* file. This file contains information about additional tabs of toolbar additional buttons and their icons and filenames of containers. The contents of this file are clear and omitted here.

Creation or edition of the *Containers.xml* file results in the appearing of new tabs and buttons, like it is presented in the following picture:

We have an additional button, and can put new components on the desktop.

## 5. Motion Model

The motion model of an artificial satellite is a system of ordinary differential equations that describe the motion in a Greenwich reference frame:

where are the coordinates and components of velocity respectively, is angular velocity of Earth, and are components of the summary external acceleration. The current model contains gravitational force and the aerodynamic one. These accelerations can be defined by the following formulae:

The following table contains meaning of these parameters:

Component of ordinary differential equations is used for above equations. This component does not have a rich system of symbols and uses lower case variables only. However this component supports rich system of comments which is presented below:

These comments can contain mappings between variables and physical parameters. Independent variables of equations (*x*, *y*, *z*, *u*, *v*, *w*) are checked in right part. Mapping between these variables and

The component represents differential equations in the following way:

Besides independent variables equations contain additional parameters. The following table contains mapping of these parameters:

Some of these parameters are constants. Other parameters should be exported. Constants are checked and exported parameters are not checked below:

Exported parameters concern with gravity and atmosphere models. These models are contained in external libraries. The following picture contains components of motion model.

The **Motion equations** component is considered above component of ordinary differential equations. These equations have feedback. Right parts of equations depend on gravitational forces. Otherwise gravitational forces depend on coordinates of the satellite. The **Vector** component is used as feedback of **Motion equations**.

This picture has the following meaning. Variables *x*, *y*, *z*, *u*, *v*, *w* of **Vector** component are assigned to *x*, *y*, *z*, *u*, *v*, *w* variables of **Motion equations**. Input and output parameters of **Vector** are presented below:

Components **Gravity**, **Atmosphere** are objects obtained from external libraries which calculate Gravity of Earth and Earth's atmosphere respectively. Both these objects implement `IObjectTransformer`

interface. The **G-transformation** and **A-transformation** objects are helper objects which are used for application of **Gravity**, **Atmosphere** objects respectively. Properties of **G-transformation** are presented below:

The above picture has the following meaning. The **G-transformation** is connected to **Gravity** object. This object corresponds to exported `Gravity`

class (code of `Gravity`

class is above this text). The `Gravity`

class implements `IObjectTransformer`

interface. Input parameters of `Gravity`

are coordinates which are labeled by "*x*", "*y*", "*z*" labels.So above comboboxes are labeled by *x*", "*y*", "*z* labels. The above picture means that parameters *x*", "*y*", "*z*" are assigned to *Formula_1*, *Formula_2*, *Formula_3* of **Vector** component. Otherwise *Formula_1*, *Formula_2*, *Formula_3* of **Vector** are assigned to *x*, *y*, *z* of **Motion equations**. However **Motion equations** component uses parameters of **G-transformation**:

This picture means that parameters *a*, *b*, *c* of *a*, of **Motion equations** are assigned to parameters *Gx*, *Gy*, *Gz* of **G-transformation**. So of Earth's `Gravity `

function is exported to **G-transformation**. The Earth's atmosphere is exported in a similar way.

## 6. Description of the Problem

The following picture illustrates the problem of determination of the orbit.

We have an artificial Earth's satellite. There are a set of ground stations which perform measurements. So we have measurements of distance and relative velocity between stations and the satellite. Processing of these measurements provides determination of orbit.

### 6.1 Input Data

Input data contains measurements of distances and relative velocities. Special components are used for this purpose. These components are presented below:

Ordinates of **Distance** component are measured distances and ordinates of **Velocity** one are measured relative velocities. Abscises of these components are indifferent. Other input data are represented below:

The full list of input parameters is represented in the following table:

| Component name | Comment |

1 | Distance | Array of distance measurements |

2 | Velocity | Array of velocity measurements |

3 | Distance Time | Array of times of distance measurements |

4 | Velocity Time | Array of times of velocity measurements |

5 | Distance X | Array of *X* - coordinates of ground station which correspond to distance measurements |

6 | Distance Y | Array of *Y* - coordinates of ground station which correspond to distance measurements |

7 | Distance Z | Array of *Z* - coordinates of ground station which correspond to distance measurements |

8 | Velocity X | Array of *X* - coordinates of ground station which correspond to velocity measurements |

9 | Velocity Y | Array of *Y* - coordinates of ground station which correspond to velocity measurements |

10 | Velocity Z | Array of *Z* - coordinates of ground station which correspond to velocity measurements |

11 | Distance Sigma | Array of standard deviations of distance measurements |

12 | Velocity Sigma | Array of standard deviations of velocity measurements |

Let us comment on this table. We have two arrays of measurements. First array *d*_{1}, *d*_{2}, ... contains measurements of distance. Second one *V*_{1}, *V*_{2} ... contains measurements of velocity. First two rows of the above table correspond to these two arrays. Measurements *d*_{1}, *d*_{2}, ... are performed in different times. Third row of above table corresponds to array of times of these measurements. Similarly 4 - th row corresponds to array of times of *V*_{1}, *V*_{2} ... measurements. A set of ground stations is used for measurements. Suppose that *d*_{i} is obtained by ground station which has coordinates *X*_{i}, *Y*_{i}, *Z*_{i}. The 5, 6, 7 th rows of table correspond to *X*_{i}, *Y*_{i}, *Z*_{i} arrays. Similarly 8 - 10 rows correspond to coordinates of ground stations and *V*_{i} measurements. Note that the different measurements have no equal accuracy. The General linear method requires to use different weights for such measurements. A new component was developed for this purpose. It enables us to create weighted selections. In case of orbit determination, measurements have different variances. So we have weighted least squares situation. So variances of measurements are needed. Rows 11 and 12 contain arrays variances of *d*_{i} and *V*_{i} respectively. Components **Distance Selection** and **Velocity Selection** are "weighted selections" which contain measurements and their variances. Properties of **Distance Selection** are presented below:

This picture has the following meaning:

- Array of measurements is contained in ordinates of
**Distance** component.
- Array of variances is contained in ordinates of
**Distance Sigma** component.

The **Velocity Selection** component is constructed in a similar way. I've used the graph component for the storage of the selection. It has been done since the user would be able to try this example without a database. Really, this scenario should use the connection to the database graph .

### 6.2 Nonlinear Regression

**Nonlinear regression** in statistics is the problem of fitting a model.

to multidimensional *x*,*y* data, where *f* is a nonlinear function of *x*, with regression parameter θ.

Nonlinear regression operates with selections. This software implements two methods of manipulation with selections. The first method loads selection iteratively, i.e., step by step. The second one loads selection at once.

The architecture of nonlinear regression software that uses iterative regression is presented in the following scheme:

Let us describe the components of this scheme. The iterator provides data-in selections *x* and *y*. The *y* is the **Left part** of fitting the equations. The **Transformation** corresponds to the nonlinear function *f*, and generates the **Left part** of the fitting model. The **Processor** coordinates all the actions and corrects the **Regression parameters**.

The architecture of the software with loading selection at once is presented in the following scheme:

The processor compares **Calculated parameters** with **Selection**, calculates residuals, and then corrects **Regression parameters**. In these two schemas, the **Iterator**, **Selection**, and others are not concrete components. They are components that implement interfaces of **Iterator**, **Selection** etc. For example, a **Camera** component may implement the **Selection** interface. Here, I'll consider examples in brief. A more profound description is available in the documentation of this project. Second scheme is used for determination of orbits.

### 6.3 Ingredients

The above picture contains general ingredients of nonlinear regression. Let us consider ingredients related to orbit determination of orbits. The **Selection** is already considered in 6.1. **Regression parameters** are initial values of coordinates and components of velocity of the satellite. These parameters are contained in **Motion equations** component. These parameters are presented in the following picture:

As it is already noted, parameters *x*, *y*, *x*, *u*, *v*, *w* are coordinates and components of velocity. Top part of above picture contains values of these parameters before orbit determination. Bottom one contains values after orbit determination. These parameters are being changed during orbit determination. Transformation stage contains the following steps:

- Calculation of trajectory of satellite
- Calculation of time shifts
- Calculation of distances and velocities

The **Accumulator** component performs calculation of trajectory parameters as actual functions. Let us explain what is the actual function. There are two points of definitions of function. According to potential point function *f*(*t*) is calculation of *f* by known argument *t*. Actual point means that function *f* is single (actual) object. Actual function can be obtained from potential by calculation of values of *f*(*t*) at different values of *t*. The **Motion equations** component performs potential calculation of motion parameters and **Accumulator** object transforms these potential functions to actual ones. Properties of **Accumulator** are presented below:

According to these properties, **Accumulator** calculates *f*(*t*_{0}), *f*(*t*_{1}), *f*(*t*_{2}), ... where *t*_{0} = 1770457504 (Start), *t*_{i+1} - *t*_{i} = 10 (Step). Number of values of argument is equal to 2500 (Step count). So **Accumulator** stores *f*(*t*_{0}), *f*(*t*_{1}), *f*(*t*_{2}) and for any *t*_{0} < *t* < *t*_{2500}. Then *f*(*t*) is obtained by interpolation by polynom of third degree. It should be noted that measurements of distances and velocities are not instant since speed of light is finite. Ground station sends signal at *t*_{1} time, then satellite sends reflected signal at *t*_{2} time and then ground station receives signal at *t*_{3} time. The time of measurement is *t*_{3}. However measurement of distance is not distance between satellite and station at *t*_{3}. It is equal (*d*_{1} - *d*_{2}) / 2 where *d*_{1} is distance between station at *t*_{1} time and satellite at *t*_{2} time and similarly *d*_{2} is distance between station at *t*_{3} time and satellite at *t*_{2} time. The **Difference** and **Time + Shift** components enable us take into account this circumstance. Properties of **Difference** component are presented in the following picture:

Output of this component is used by **Time + Shift**. Properties of **Time + Shift** are presented below:

This component calculates arrays of times of reflection of signals by the satellites. Component **Delta** is helper. It has the following properties:

The **Measurements** component calculates measurements of distance and velocity. This component has the following properties:

This parameters are compared with selections. The **Processor** component performs such computation or (using another term) regression. Regression formula is presented below:

**Processor** component reflects this formula. Properties of this component are presented below:

Left part of above window corresponds to *x* of regression formula. In considered case components of *x* are coordinates and components of velocity. The following table contains components of vector *x*.

| Physical parameter | Identifier |

1 | *X* | *Motion Equations/Motion Equations.x* |

2 | *Y* | *Motion Equations/Motion Equations.y* |

3 | *Z* | *Motion Equations/Motion Equations.z* |

4 | *V*_{x} | *Motion Equations/Motion Equations.u* |

5 | *V*_{y} | *Motion Equations/Motion Equations.v* |

6 | *V*_{z} | *Motion Equations/Motion Equations.w* |

Middle part of above window contains *f*(*x*) of regression formula. The *f*(*x*) vector contains two subvectors *Formula_1* and *Formula_2* which belong to **Measurements** components. These subvectors are calculated values of distances and velocities. Right part of above window corresponds to *y* of regression formula. Vector *y* contains two selections which correspond to distance and velocity measurements respectively. Click to *Iterate* button initiate step of nonlinear regression.

## 7. Kalman Filter

This section contains the technology of Kalman filter without any determination of orbits. Common theory and simple example is considered only.

### 7.1 Theory

The Kalman filter is a mathematical method named after Rudolf E. Kalman. Its purpose is to use measurements that are observed over time that contain noise (random variations) and other inaccuracies, and produce values that tend to be closer to the true values of the measurements and their associated calculated values. Kalman filter theory does not depend on physical phenomenon. So Kalman filter can be used for estimation of state different types of objects (mechanical, electrical, pneumatic, etc.). Moreover Kalman filter theory does not depend on types of measurements. Any type of measurements are acceptable (inertial, electrical, temperature, optical, nuclear). So Kalman Filter is engineering object which is invariant to physical background. We have universal engineering object and this idea is very close to idea of universal engineering software. How can Kalman Filter theory exist without any physical background? It is based on abstract equations. Some of these equations are presented below:

where *x*_{k- 1}, *x*_{k} are state vectors at *k* - 1and *k*-th steps respectively. The *z*_{k} is vector of measurements. Vector functions *f* and *h* are respectively called transition function and function of measurements. Vectors *w*_{k} and *v*_{k}. All above notions are abstract. These notions do not directly related to any special engineering problem. Kalman Filter uses matrixes of partial derivations:

### 7.2 Implementation

Kalman filter implementation contains the following set of matrix formulae:

Implementation of these formulae is presented below:

Kalman filter uses operations with vectors and matrixes. The vector as matrix components has been developed for this purpose. A typical example of the usage of the Kalman filter is presented below:

This picture has all the necessary Kalman filter ingredients:

- Covariance matrix
- Gain derivation
- Relationship to recursive Bayesian estimation

The editor for the properties of the matrix is presented below:

It links the elements of the matrix with external data. Operations with matrixes is being performed by a formula editor. As it had been noted in part 4, the editor is case sensitive. For example, the meaning of the formula *ax* may be different. If *a* and *x* are real numbers, then this formula means an ordinary product. If *a* and *x* are matrixes, then we have a matrix product. If *a* is a matrix and *x* is a vector, then *ax* is a matrix - vector product etc. In the following formulas of this situation:

*a* and *h* are matrixes and *y* and *z* are vectors. So having vector and matrix components, we can construct any modification of the Kalman filter.

## 8. Forecast of Artificial Earth's Satellite

Besides calculation, any software needs generation of documents. Facilities generation of documents are exhibited on forecast problem. Forecast document contains parameters of satellite motion in times of equator intersections. The following picture explains the forecast problem:

Satellite moves along its orbit. Parameters of motion correspond to the time when it moves from South to North and intersects equator. Satellite intersects equator when its *Z* - coordinate is equal to 0. If we have coordinates and time as function of *Z*, then we can see necessary parameters are values of these functions where *Z* = 0. The following picture represents a solution of this task.

The **Delay** component has `FunctionAccumulator`

type. This type is very similar to Transport delay of Simulink. However `FunctionAccumulator`

is more advanced. Properties of **Delay** are presented below:

The **Delay** component is connected to **Data** to data one. It stores latest output values of **Data**. The *Step count* = 4 means that **Delay** stores 4 last values of each output parameters. The *Formula_4* of data is *Z* coordinate of satellite. So all parameters of **Delay** are considered as functions of *Z*. Let us consider all components of forecast. The **Satellite** component is a container which contains **Motion equations** of the satellite. The **Data** component is connected to **Motion equations** and adds time parameters to motion parameters. Properties of **Data** are presented below:

.

The **Recursive** component is similar to Memory component of Simulink. This component calculates time delay of *Z* coordinate. Properties of the **Recursive** element are presented below:

The **Condition** element is used for determination of time of equator intersection. Properties of **Condition** are presented below:

Meaning of parameters of this component is contained in comments:

And at last, **Output** component contains calculation of forecast parameters. It has the following formulae:

The 86 400 number is number of seconds in a day. The *time* function converts `double`

variable to `DateTime`

. The *b* is condition of equator intersection. Input parameters are motion parameters as functions of . Full list of parameters is contained in the following comments:

And at last, **Chart** component has the following properties:

These properties have the following meaning. Output is performed when condition variable is `true`

. In this case, condition is *Formula_1* of **Condition**

object. It is condition of equator intersection. Output parameters are *Formula_1* - *Formula_7* of **Output**

object. Labels *X*, *Y*, *Z*, *Vx*, *Vy*, *Vz* are labels which are reflected in XML document. The XML document is presented below:

<Root>
<Parameters>
<Parameter Name="T" Value="Monday, June 28, 1999 11:17:05 PM.457" />
<Parameter Name="X" Value="-6595.47050815637" />
<Parameter Name="Y" Value="-2472.05799683873" />
<Parameter Name="Z" Value="0" />
<Parameter Name="Vx" Value="-0.541295154046453" />
<Parameter Name="Vy" Value="1.46945520971716" />
<Parameter Name="Vz" Value="7.45051367374592" />
</Parameters>
<Parameters>
<Parameter Name="T" Value="Tuesday, June 29, 1999 12:55:08 AM.068" />
<Parameter Name="X" Value="-7026.76544757067" />
<Parameter Name="Y" Value="487.053173194772" />
<Parameter Name="Z" Value="0" />
<Parameter Name="Vx" Value="0.117122882887633" />
<Parameter Name="Vy" Value="1.56173832654579" />
<Parameter Name="Vz" Value="7.45055871464937" />
</Parameters>
<Root>

The first parameter has `DateTime`

type and the other ones have `double`

type. The XML document reflects this circumstance. The XSLT (Extensible Stylesheet Language Transformations) is used for generation of documents. The following code contains XSLT file.

<xsl:stylesheet version="1.0" xmlns:xsl=http://www.w3.org/1999/XSL/Transform
xmlns:msxsl="urn:schemas-microsoft-com:xslt"
xmlns:user="urn:my-scripts">
<xsl:template match="/">
<html>
<body>
<p>
<h1>Forecast</h1>
</p>
<xsl:apply-templates />
</body>
</html>
</xsl:template>
<xsl:template match="Parameters">
<p></p>
<p>
<table border="1">
<xsl:apply-templates />
</table>
</p>
</xsl:template>
<xsl:template match="Parameter">
<xsl:apply-templates />
<tr>
<td>
<xsl:value-of select="@Name" />
</td>
<td>
<xsl:value-of select="@Value" />
</td>
</tr>
</xsl:template>
</xsl:stylesheet>

This XSLT code provides the following HTML document.

## Forecast

T | Monday, June 28, 1999 11:17:05 PM.457 |

X | -6595.47050815637 |

Y | -2472.05799683873 |

Z | 0 |

Vx | -0.541295154046453 |

Vy | 1.46945520971716 |

Vz | 7.45051367374592 |

T | Tuesday, June 29, 1999 12:55:08 AM.068 |

X | -7026.76544757067 |

Y | 487.053173194772 |

Z | 0 |

Vx | 0.117122882887633 |

Vy | 1.56173832654579 |

Vz | 7.45055871464937 |

## Points of Interest

There exists an opinion that deep research requires particular specialty. However, sir Isaac Newton was a great scientist and a master of the Mint. Indeed, the occupation of abstract algebra (part 5) is not an obstacle for everyday engineering problems.