Universal Framework for Science and Engineering - Part 6: Determination of Orbits of Artificial Satellites

• Download source code - 3.21 MB
• Download gravity field project - 218.48 KB
• Download atmosphere density model project - 224.47 KB
• Download installer - 1.7 MB
• Download Astroframe project documentation - 1.73 MB
• Download sample of motion model of artificial satellite - 187.07 KB
• Download sample of General linear method - 230.29 KB
• Download sample of Kalman filter - 4.63 KB
• Download sample of forecast - 129.65 KB

Useful Links

• The Universal Framework for Science and Engineering Homepage
• Project Astroframe homepage
• Visual C# 2010 Express

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:

/// <summary>
/// Calculation of gravity accelerations
/// </summary>
/// <param name="N0">Degree of Legendre polinomials</param>
/// <param name="NK">Number of harmonics</param>
/// <param name="X">X - coordinate</param>
/// <param name="Y">Y - coordinate</param>
/// <param name="Z">Z - coordinate</param>
/// <param name="FX">X - component of gravitational acceleration</param>
/// <param name="FY">Y - component of gravitational acceleration</param>
/// <param name="FZ">Z - component of gravitational acceleration</param>
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.

/// <summary>
/// Transformer of objects
/// </summary>
public interface IObjectTransformer
{
    /// <summary>
    /// Input variables
    /// </summary>
    string[] Input
    {
        get;
    }

    /// <summary>
    /// Output variables
    /// </summary>
    string[] Output
    {
        get;
    }

    /// <summary>
    /// Gets type of i - th input variable
    /// </summary>
    /// <param name="i">Variable index</param>
    /// <returns>The type</returns>
    object GetInputType(int i);

    /// <summary>
    /// Gets type of i - th output variable
    /// </summary>
    /// <param name="i">Variable index</param>
    /// <returns>The type</returns>
    object GetOutputType(int i);

    /// <summary>
    /// Calculation
    /// </summary>
    /// <param name="input">Input</param>
    /// <param name="output">Output</param>
    void Calculate(object[] input, object[] output);
}

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

The following code snippet shows implementation of this interface:

/// <summary>
/// Gravity of Earth
/// </summary>
public class Gravity : IObjectTransformer, ICategoryObject, IPropertiesEditor,
    IObjectFactory
 {
    /// <summary>
    /// Inputs
    /// </summary>
    static private readonly string[] inps = new string[] { "x", "y", "z" };

    /// <summary>
    /// Outputs
    /// </summary>
    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;
    }

    /// <summary>
    /// Calculation
    /// </summary>
    /// <param name="input">Input</param>
    /// <param name="output">Output</param>
    void IObjectTransformer.Calculate(object[] input, object[] output)
    {
        // Input cast
        double x = (double)input[0];
        double y = (double)input[1];
        double z = (double)input[2];

        double fx, fy, fz;

        // Call of forces
        Forces(n0, nk, x, y, z, out fx, out fy, out fz);

        // filling of output
        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:

/// <summary>
/// The factory of category object
/// </summary>
public interface IObjectFactory
{
    /// <summary>
    /// Names of objects
    /// </summary>
    string[] Names
    {
        get;
    }

    /// <summary>
    /// Object by name
    /// </summary>
    /// <param name="name">The name</param>
    /// <returns>The object</returns>
    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:

/// <summary>
/// Gets factory from assembly
/// </summary>
/// <param name="ass">The assembly</param>
/// <returns>The factory</returns>
static public IObjectFactory GetFactory(Assembly ass)
{
    Type[] types = ass.GetTypes(); // Types of assembly
    foreach (Type type in types)
    {
        if (type.GetInterface("IObjectFactory") != null)
        {
            // Looking for singleton
            FieldInfo fi = type.GetField("Object");
            if (fi != null)
            {
                return fi.GetValue("Object") as IObjectFactory;
            }
            else
            {
                // Calling of default constructor
                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.

• Click here to download a sample of the container

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:

Let us comment on this table. We have two arrays of measurements. First array d₁, d₂, ... contains measurements of distance. Second one V₁, V₂ ... contains measurements of velocity. First two rows of the above table correspond to these two arrays. Measurements d₁, d₂, ... 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₁, V₂ ... 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₀), f(t₁), f(t₂), ... where t₀ = 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₀), f(t₁), f(t₂) and for any t₀ < t < t₂₅₀₀. 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₁ time, then satellite sends reflected signal at t₂ time and then ground station receives signal at t₃ time. The time of measurement is t₃. However measurement of distance is not distance between satellite and station at t₃. It is equal (d₁ - d₂) / 2 where d₁ is distance between station at t₁ time and satellite at t₂ time and similarly d₂ is distance between station at t₃ time and satellite at t₂ 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.

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.

• Download this sample - 229 Kb

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:

• Click here to download this sample.

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: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

• Download this sample - 132 KB

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.