Astrodynamics library

Description

Astrodynamics is a C++ based, Python wrapped, high performance, low dependency software library. It includes algorithms for math operations, numeric techniques, and for solving problems of astrodynamics, astronomy, geodesy, geomatics and other linked disciplines.

Purpose

This is my private educational and research project. It serves 4 main goals:

This is my way to study astrodynamics and linked disciplines.
This is my very own and very independent library of tools for solving astrodynamics problems.
This project allows me to improve my coding skills.
This project also serves as my portfolio.

Features

Astrodynamics library contains tools for astrodynamics calculations, space flight simulation and mission planning. The toolset now includes more than 230 tools.

List of implemented functions (click the arrow to expand)

Set of math functions
- Extended set of basic math functions
- Set of classes for linear algebra
  - vectors
  - matrices
  - quaternions
- Set of functions for numeric methods
  - basic numeric methods (Newton-Raphson, power series)
  - extremum search methods (split-half, golden section, parabolic interpolation)
  - root-finding methods (Brent, solvers for quadratic, cubic & quartic equations)
  - Integrating routines (including solutions for elliptic integrals)
  - Legendre polynomial, ALF & fnALF solutions
  - Numeric integration routines
    - generic ODE numeric solvers (Verlet, RK4, RKF, DP, CK)
    - generic 2-nd order solvers (RKN54)
    - specialized integrators (ABM, SC, GJ, and other multistep solutions with orders 3-16)
- Spline smoothing, blending and interpolation routines
Set of application astrodynamic functions
- Basic Kepler problem solutions
  - COE <–> RV elements conversion
  - Time of flight
  - Kepler propagator
- Coordinate and time conversion functions
  - IAU2000 CIO ITRF<–>GCRF conversion functions
  - IAU 76 FK5 ITRF<–>GCRF conversion functions
  - Carthesian <–> Spherical (geocentric & geodetic) conversion functions
  - Universal UTC <–> JD conversion functions
  - UTC -> UT1, TAI, GPS, TT, TDB, etc. time conversions
- Observation processing functions
  - IJK -> RA/Dec
  - ECEF ->RA/Dec
  - SEZ<–> R/Az/El
- Celestial phenomena functions
  - Sun position, rise and set
  - Moon position, rise and set
  - Planet positions
  - Umbra/Penumbra, LoS detection
- Set of functions basic space maneuvers
  - Coplanar maneuvers
  - Non-colpanar maneuvers
  - Relative motion
    - Hills equations
    - Hill EQCM <–> ECI conversion
- Lambert problem solutions
  - Lambert problem solvers
  - Targeting algorithm
- Gravity models
  - J2-J6 model
  - EGM 2008 model (degree/order 2190 x 2159)
- Atmospheric models
  - Exponential atmosphere model
- Launch & aiming solutions
  - Launch window calculator
  - Launch azimuth calculator
Astrodynamics constant library
Set of auxiliary functions

This list is not final – it will be extended with new functional on my further studies.

Algorithms

Original algorithms are taken from works (including books, scientific articles, public documents and periodical papers) of:

David A. Vallado & Wayne D. McClain
Richard H. Battin
James R. Wertz
Oliver Montenbtuck & Eberhard Gill
Theresa & Granino Korn
Boris Demidowich
Roger R. Bate, Donald D. Mueller & Jerry E. White
Howard Curtis
and many other (will be explicitly enumerated later)

Part of algorithms are re-worked by me to apply for practical purposes. Some of algorithms I developed myself.

Architecture

The library designed with strict requirements to reliability, long-term maintainability and practical applicability.

Architecture of the library represents 3 separated components:

C++ calculation core
C ABI
Python wrapper / API

Most of heavy calculations executed inside high performance C++ core. It is made as highly-modular, easily modified set of functions, utilizes multi-threading and advanced calculation algorithms.

Lightweight C ABI provides standard unified access to the core for any third-party software.

Python wrapper provides simple and agile interface for the astrodynamics tools, allows to use them with any modules of modern Python ecosystem.

Dependencies

C++ core of Astrodynamics library has 0 (zero) external dependencies, utilizing only standard library of C++ 11. See Fig. 1:

Python based wrapper uses Python 3.9.13, with only 1 dependency (NumPy). It seemed to be impractical to avoid NumPy, because most of other Python software widely uses it.

Code

The key code features that I tried to maintain in this project were simplicity, modularity, and uniformity.

Here is an example of the code – C++ Vector class.

Vector.h:

#pragma once
#include <cmath>
#include <stdexcept>

#include "mathematics.h"

namespace Astrodynamics
{
	namespace Mathematics
	{
		class Vector
		{
		public:
			double x = 0;
			double y = 0;
			double z = 0;

			Vector(const double x = 0.0, const double y = 0.0, const double z = 0.0);

			static bool AreEqual(const Vector v1, const Vector v2, const double tolerance = 1e-9);

			bool operator == (const Vector other) const;
			bool operator != (const Vector other) const;
			Vector operator + (const Vector other) const;
			Vector operator - (const Vector other) const;
			Vector operator * (const double scalar) const;
			Vector operator / (const double scalar) const;
			double Magnitude();
			Vector Normalize();

			static double DotProduct(const Vector v1, const Vector v2);
			static Vector CrossProduct(const Vector v1, const Vector v2);
			static double AngleBetween(Vector v1, Vector v2);
			static double SignedAngleBetween(Vector v1, Vector v2, Vector a);
			static Vector ProjectOnPlane(Vector v, Vector n);

			double operator * (const Vector other) const;
			Vector operator - () const;

			double& operator[](const int axis)
			{
				switch (axis)
				{
				case 0:
					return x;
				case 1:
					return y;
				case 2:
					return z;
				default:
					throw std::out_of_range("Out of range");
				}
			}
			const double& operator[](const int axis) const //required if vector is const
			{
				switch (axis)
				{
				case 0:
					return x;
				case 1:
					return y;
				case 2:
					return z;
				default:
					throw std::out_of_range("Out of range");
				}
			}
		};

		Vector UnitVector(const int axis);

		double mag(Vector v);
		Vector norm(Vector v);
		double dot(const Vector v1, const Vector v2);
		Vector cross(const Vector v1, const Vector v2);
		double ang(Vector v1, Vector v2);
		double sgdang(Vector v1, Vector v2, Vector a);
		Vector proj(Vector v, Vector n);

		Vector operator*(const double s, const Vector v);
	}
}

#pragma once
#include <cmath>
#include <stdexcept>

#include "mathematics.h"

namespace Astrodynamics
{
	namespace Mathematics
	{
		class Vector
		{
		public:
			double x = 0;
			double y = 0;
			double z = 0;

			Vector(const double x = 0.0, const double y = 0.0, const double z = 0.0);

			static bool AreEqual(const Vector v1, const Vector v2, const double tolerance = 1e-9);

			bool operator == (const Vector other) const;
			bool operator != (const Vector other) const;
			Vector operator + (const Vector other) const;
			Vector operator - (const Vector other) const;
			Vector operator * (const double scalar) const;
			Vector operator / (const double scalar) const;
			double Magnitude();
			Vector Normalize();

			static double DotProduct(const Vector v1, const Vector v2);
			static Vector CrossProduct(const Vector v1, const Vector v2);
			static double AngleBetween(Vector v1, Vector v2);
			static double SignedAngleBetween(Vector v1, Vector v2, Vector a);
			static Vector ProjectOnPlane(Vector v, Vector n);

			double operator * (const Vector other) const;
			Vector operator - () const;

			double& operator[](const int axis)
			{
				switch (axis)
				{
				case 0:
					return x;
				case 1:
					return y;
				case 2:
					return z;
				default:
					throw std::out_of_range("Out of range");
				}
			}
			const double& operator[](const int axis) const //required if vector is const
			{
				switch (axis)
				{
				case 0:
					return x;
				case 1:
					return y;
				case 2:
					return z;
				default:
					throw std::out_of_range("Out of range");
				}
			}
		};

		Vector UnitVector(const int axis);

		double mag(Vector v);
		Vector norm(Vector v);
		double dot(const Vector v1, const Vector v2);
		Vector cross(const Vector v1, const Vector v2);
		double ang(Vector v1, Vector v2);
		double sgdang(Vector v1, Vector v2, Vector a);
		Vector proj(Vector v, Vector n);

		Vector operator*(const double s, const Vector v);
	}
}

Vector.cpp:

#include "mathematics_Vector.h"

using namespace Astrodynamics::Mathematics;

Vector::Vector(const double x, const double y, const double z)
{
	this->x = x;
	this->y = y;
	this->z = z;
}

bool Vector::AreEqual(const Vector v1, const Vector v2, const double tolerance)
{
	return (Equal(v1.x, v2.x, tolerance) && Equal(v1.y, v2.y, tolerance) && Equal(v1.z, v2.z, tolerance));
}

bool Vector::operator==(const Vector other) const
{
	return AreEqual(*this, other);
}

bool Vector::operator!=(const Vector other) const
{
	return !(*this == other);
}

Vector Vector::operator+(const Vector other) const
{
	return Vector(x+other.x, y + other.y, z + other.z);
}

Vector Vector::operator-(const Vector other) const
{
	return Vector(x - other.x, y - other.y, z - other.z);
}

Vector Vector::operator*(const double scalar) const
{
	return Vector(x * scalar, y * scalar, z * scalar);
}

Vector Vector::operator/(const double scalar) const
{
	if (scalar != 0.0)
	{
		return Vector(x / scalar, y / scalar, z / scalar);
	}
	else
	{
		return Vector();
	}
}

double Vector::Magnitude()
{
	return sqrt(pow(x, 2.0) + pow(y, 2.0) + pow(z, 2.0));
}

Vector Vector::Normalize()
{
	double magnitude = Magnitude();
	if (magnitude > 0.0)
	{
		return *this / magnitude;
	}
	else
	{
		return Vector();
	}
}

double Vector::DotProduct(const Vector v1, const Vector v2)
{
	return (v1.x * v2.x) + (v1.y * v2.y) + (v1.z * v2.z);
}

Vector Vector::CrossProduct(const Vector v1, const Vector v2)
{
	return Vector
	(
		(v1.y * v2.z) - (v1.z * v2.y),
		(v1.z * v2.x) - (v1.x * v2.z),
		(v1.x * v2.y) - (v1.y * v2.x)
	);
}

double Vector::AngleBetween(Vector v1, Vector v2)
{
	double dot = DotProduct(v1, v2);
	double magnitudeProduct = v1.Magnitude() * v2.Magnitude();
	if (magnitudeProduct > 0.0)
	{
		return degrees(acos(std::max(-1.0, std::min(1.0, dot / magnitudeProduct))));
	}
	else
	{
		return 0.0;
	}
}

double Vector::SignedAngleBetween(Vector v1, Vector v2, Vector a)
{
	Vector aNormalized = a.Normalize(); //always must be normalized for correct results
	double dot = DotProduct(v1, v2);
	Vector cross = CrossProduct(v1, v2);
	return degrees(atan2(DotProduct(cross, aNormalized), dot));
}

Vector Vector::ProjectOnPlane(Vector v, Vector n)
{
	Vector nNormalized = n.Normalize();
	Vector p = nNormalized * DotProduct(v, nNormalized);
	return v - p;
}

double Vector::operator*(const Vector other) const
{
	return Vector::DotProduct(*this, other);
}

Vector Vector::operator-() const
{
	return Vector(-x, -y, -z);
}

Vector Astrodynamics::Mathematics::UnitVector(const int axis)
{
	switch (axis)
	{
	case 0:
		return Vector(1, 0, 0);
	case 1:
		return Vector(0, 1, 0);
	case 2:
		return Vector(0, 0, 1);
	default:
		throw std::out_of_range("Out of range");
	}
}

double Astrodynamics::Mathematics::mag(Vector v)
{
	return v.Magnitude();
}

Vector Astrodynamics::Mathematics::norm(Vector v)
{
	return v.Normalize();
}

double Astrodynamics::Mathematics::dot(const Vector v1, const Vector v2)
{
	return Vector::DotProduct(v1, v2);
}

Vector Astrodynamics::Mathematics::cross(const Vector v1, const Vector v2)
{
	return Vector::CrossProduct(v1, v2);
}

double Astrodynamics::Mathematics::ang(Vector v1, Vector v2)
{
	return Vector::AngleBetween(v1, v2);
}

double Astrodynamics::Mathematics::sgdang(Vector v1, Vector v2, Vector a)
{
	return Vector::SignedAngleBetween(v1, v2, a);
}

Vector Astrodynamics::Mathematics::proj(Vector v, Vector n)
{
	return Vector::ProjectOnPlane(v, n);
}

Vector Astrodynamics::Mathematics::operator*(const double s, const Vector v)
{
	return v * s;
}

#include "mathematics_Vector.h"

using namespace Astrodynamics::Mathematics;

Vector::Vector(const double x, const double y, const double z)
{
	this->x = x;
	this->y = y;
	this->z = z;
}

bool Vector::AreEqual(const Vector v1, const Vector v2, const double tolerance)
{
	return (Equal(v1.x, v2.x, tolerance) && Equal(v1.y, v2.y, tolerance) && Equal(v1.z, v2.z, tolerance));
}

bool Vector::operator==(const Vector other) const
{
	return AreEqual(*this, other);
}

bool Vector::operator!=(const Vector other) const
{
	return !(*this == other);
}

Vector Vector::operator+(const Vector other) const
{
	return Vector(x+other.x, y + other.y, z + other.z);
}

Vector Vector::operator-(const Vector other) const
{
	return Vector(x - other.x, y - other.y, z - other.z);
}

Vector Vector::operator*(const double scalar) const
{
	return Vector(x * scalar, y * scalar, z * scalar);
}

Vector Vector::operator/(const double scalar) const
{
	if (scalar != 0.0)
	{
		return Vector(x / scalar, y / scalar, z / scalar);
	}
	else
	{
		return Vector();
	}
}

double Vector::Magnitude()
{
	return sqrt(pow(x, 2.0) + pow(y, 2.0) + pow(z, 2.0));
}

Vector Vector::Normalize()
{
	double magnitude = Magnitude();
	if (magnitude > 0.0)
	{
		return *this / magnitude;
	}
	else
	{
		return Vector();
	}
}

double Vector::DotProduct(const Vector v1, const Vector v2)
{
	return (v1.x * v2.x) + (v1.y * v2.y) + (v1.z * v2.z);
}

Vector Vector::CrossProduct(const Vector v1, const Vector v2)
{
	return Vector
	(
		(v1.y * v2.z) - (v1.z * v2.y),
		(v1.z * v2.x) - (v1.x * v2.z),
		(v1.x * v2.y) - (v1.y * v2.x)
	);
}

double Vector::AngleBetween(Vector v1, Vector v2)
{
	double dot = DotProduct(v1, v2);
	double magnitudeProduct = v1.Magnitude() * v2.Magnitude();
	if (magnitudeProduct > 0.0)
	{
		return degrees(acos(std::max(-1.0, std::min(1.0, dot / magnitudeProduct))));
	}
	else
	{
		return 0.0;
	}
}

double Vector::SignedAngleBetween(Vector v1, Vector v2, Vector a)
{
	Vector aNormalized = a.Normalize(); //always must be normalized for correct results
	double dot = DotProduct(v1, v2);
	Vector cross = CrossProduct(v1, v2);
	return degrees(atan2(DotProduct(cross, aNormalized), dot));
}

Vector Vector::ProjectOnPlane(Vector v, Vector n)
{
	Vector nNormalized = n.Normalize();
	Vector p = nNormalized * DotProduct(v, nNormalized);
	return v - p;
}

double Vector::operator*(const Vector other) const
{
	return Vector::DotProduct(*this, other);
}

Vector Vector::operator-() const
{
	return Vector(-x, -y, -z);
}

Vector Astrodynamics::Mathematics::UnitVector(const int axis)
{
	switch (axis)
	{
	case 0:
		return Vector(1, 0, 0);
	case 1:
		return Vector(0, 1, 0);
	case 2:
		return Vector(0, 0, 1);
	default:
		throw std::out_of_range("Out of range");
	}
}

double Astrodynamics::Mathematics::mag(Vector v)
{
	return v.Magnitude();
}

Vector Astrodynamics::Mathematics::norm(Vector v)
{
	return v.Normalize();
}

double Astrodynamics::Mathematics::dot(const Vector v1, const Vector v2)
{
	return Vector::DotProduct(v1, v2);
}

Vector Astrodynamics::Mathematics::cross(const Vector v1, const Vector v2)
{
	return Vector::CrossProduct(v1, v2);
}

double Astrodynamics::Mathematics::ang(Vector v1, Vector v2)
{
	return Vector::AngleBetween(v1, v2);
}

double Astrodynamics::Mathematics::sgdang(Vector v1, Vector v2, Vector a)
{
	return Vector::SignedAngleBetween(v1, v2, a);
}

Vector Astrodynamics::Mathematics::proj(Vector v, Vector n)
{
	return Vector::ProjectOnPlane(v, n);
}

Vector Astrodynamics::Mathematics::operator*(const double s, const Vector v)
{
	return v * s;
}

I’m well aware that my code is primitive and crude. I know it looks like a student’s code. But this code works, and it will always work. And most importantly, anyone who comes after me will understand it.

Reliability

The library has 99% test coverage.

Here is Python coverage output:

Name                                    Stmts   Miss  Cover
-----------------------------------------------------------
setup.py                                    6      1    83%
test_astrodynamics.py                     925      1    99%
test_blend.py                              92      1    99%
test_elliptic.py                           12      1    92%
test_frame.py                              22      1    95%
test_integrate.py                          17      1    94%
test_interface.py                        1016      1    99%
test_interpolate.py                        20      1    95%
test_linear.py                            130      1    99%
test_mathematics.py                        96      1    99%
test_numeric.py                            99      1    99%
test_numeric_integration.py               185      1    99%
test_numeric_integration_multistep.py     183     17    91%
test_quaternion.py                        192      1    99%
test_sanity.py                              7      1    86%
test_solve.py                              34      1    97%
test_spline.py                             75      1    99%
-----------------------------------------------------------
TOTAL                                    3111     33    99%

Name                                    Stmts   Miss  Cover
-----------------------------------------------------------
setup.py                                    6      1    83%
test_astrodynamics.py                     925      1    99%
test_blend.py                              92      1    99%
test_elliptic.py                           12      1    92%
test_frame.py                              22      1    95%
test_integrate.py                          17      1    94%
test_interface.py                        1016      1    99%
test_interpolate.py                        20      1    95%
test_linear.py                            130      1    99%
test_mathematics.py                        96      1    99%
test_numeric.py                            99      1    99%
test_numeric_integration.py               185      1    99%
test_numeric_integration_multistep.py     183     17    91%
test_quaternion.py                        192      1    99%
test_sanity.py                              7      1    86%
test_solve.py                              34      1    97%
test_spline.py                             75      1    99%
-----------------------------------------------------------
TOTAL                                    3111     33    99%

C++ test coverage can not be provided, because Visual Studio Community 2022 does not provide such functional. But tests of C++ core are 100% identical to those of Python wrapper.

As an example, here is Python implementation of test for Algorithm 08: Kepler from “Fundamentals of Astrodynamics and Applications”, Fourth Edition, by David A. Vallado:

def test_Kepler(self):
        vector_r0 = vector(1131.340, -2282.343, 6672.423)
        vector_v0 = vector(-5.64305, 4.30333, 2.42879)
        dt = 40.0 * 60.0
        vector_r, vector_v = Kepler(vector_r0, vector_v0, dt, mu)
        self.assertTrue(eq(vector_r, vector(-4219.7527, 4363.0292, -3958.7666), 1e-3))
        self.assertTrue(eq(vector_v, vector(3.689866, -1.916735, -6.112511), 1e-5))

def test_Kepler(self):
        vector_r0 = vector(1131.340, -2282.343, 6672.423)
        vector_v0 = vector(-5.64305, 4.30333, 2.42879)
        dt = 40.0 * 60.0
        vector_r, vector_v = Kepler(vector_r0, vector_v0, dt, mu)
        self.assertTrue(eq(vector_r, vector(-4219.7527, 4363.0292, -3958.7666), 1e-3))
        self.assertTrue(eq(vector_v, vector(3.689866, -1.916735, -6.112511), 1e-5))

And here is its C++ analog:

		TEST_METHOD(Test08_Kepler)
		{
			Vector vector_r0 = Vector(1131.340, -2282.343, 6672.423);
			Vector vector_v0 = Vector(-5.64305, 4.30333, 2.42879);
			double dt = 40.0 * 60.0;
			Vector vector_r, vector_v;
			Kepler(vector_r, vector_v, vector_r0, vector_v0, dt, mu);
			Assert::IsTrue(Vector::AreEqual(vector_r, Vector(-4219.7527, 4363.0292, -3958.7666), 1e-3));
			Assert::IsTrue(Vector::AreEqual(vector_v, Vector(3.689866, -1.916735, -6.112511), 1e-5));
		}

		TEST_METHOD(Test08_Kepler)
		{
			Vector vector_r0 = Vector(1131.340, -2282.343, 6672.423);
			Vector vector_v0 = Vector(-5.64305, 4.30333, 2.42879);
			double dt = 40.0 * 60.0;
			Vector vector_r, vector_v;
			Kepler(vector_r, vector_v, vector_r0, vector_v0, dt, mu);
			Assert::IsTrue(Vector::AreEqual(vector_r, Vector(-4219.7527, 4363.0292, -3958.7666), 1e-3));
			Assert::IsTrue(Vector::AreEqual(vector_v, Vector(3.689866, -1.916735, -6.112511), 1e-5));
		}

As you can see, the tests are similar.

All functions of the library are checked using tests from textbooks which multiple times proved their correctness. I wrote a separate article dedicated to validation of the library algorithms utilize exercises from textbook “Satellite Orbits. Models, Methods and Applications” © 2000, by Oliver Montenbruck & Eberhard Gill.

In addition, different versions of algorithms from Astrodynamics library were practically tested number of times during my previous activities.

Extendability

Fixed C ABI and simple interface to Python provides to Astrodynamics library high extendability. I used this for make number of effective extensions, including Python-based simple Astrodynamics plotter, based on NumPy & MatplotLib, which I widely use in my studies. On Fig. 2 – 5 presented examples of plots made by this app.

Further development

As I said in the beginning of this post, Astrodynamics library is my private project I am using for study astrodynamics and navigation in space. So I am planning to improve its functionality and accuracy. You can find more information about this project on other pages of my site.

Legal information

Source code examples were published for educational and evaluation purposes only. The source code may be viewed and studied, but may not be copied, modified, redistributed, or used in any commercial, governmental, or revenue-generating projects without explicit written permission from the author.

Presented code is not intended for use in commercial or operational systems.

More information can be obtained in ULTRADYNE.SPACE LEGAL NOTICE AND TERMS OF USE.

ultradyne.space

Jacchia 71 atmosphere model

EGM 2008 2190 x 2159

Multi-step integration – practical test