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 ith input variable 
4 
GetOutputType(int i) 
Type of ith 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 Gtransformation and Atransformation objects are helper objects which are used for application of Gravity, Atmosphere objects respectively. Properties of Gtransformation are presented below:
The above picture has the following meaning. The Gtransformation 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 Gtransformation:
This picture means that parameters a, b, c of a, of Motion equations are assigned to parameters Gx, Gy, Gz of Gtransformation. So of Earth's Gravity
function is exported to Gtransformation. 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 datain 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 kth 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 <sum>Z. 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:schemasmicrosoftcom:xslt"
xmlns:user="urn:myscripts">
<xsl:template match="/">
<html>
<body>
<p>
<h1>Forecast</h1>
</p>
<xsl:applytemplates />
</body>
</html>
</xsl:template>
<xsl:template match="Parameters">
<p></p>
<p>
<table border="1">
<xsl:applytemplates />
</table>
</p>
</xsl:template>
<xsl:template match="Parameter">
<xsl:applytemplates />
<tr>
<td>
<xsl:valueof select="@Name" />
</td>
<td>
<xsl:valueof 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.