Orbital Mechanics: Part 1

tl;dr. You can find the source code here.

As an exercise in learning more of Swift, and gaining an understanding of orbital dynamics, I have started building a Swift framework that can compute orbital parameters from input parameters. The first version, really version 0.1, takes a celestial body (like the Earth), and the initial altitude and velocity to compute the periapsis (closest approach) and apoapsis (furthest away). First assumption here is that the orbiting body is small in relation to the body being orbited, so the mass of the smaller body is insignificant.

With that, the first class is the Body. At this point it contains properties necessary to describe the Body for purposes of the orbital calculation.

class Body
{
    let mass: Double// in kilograms
    let radius: Double// Of body, in meters
    let gm: Double// Universal Gravitational constant * mass of Body, in meters cubed / seconds squared
    
    init(withMass mass: Double, radius: Double, gm: Double)
    {
        self.mass = mass
        self.radius = radius
        self.gm = gm
    }
}


There are three constants associated with a Body. It’s mass, it’s radius, and GM, which is the mass times the universal gravitational constant. This term is used a lot in orbital mechanics, and some sources claim that we know the combined term with more precision than either the mass or gravitational constant. So I use the term here as it’s own value.

 

I use the Earth as an example, and as such, create an Earth class as a subclass of Body:

class Earth: Body

{
    init() {
        let mass = 5.97237e24// kilograms
        let radius = 6378100.0// meters, equatorial radius, since it's the largest
        let gm = 3.986004418e14// gravitational constant * mass of Earth, in meters cubed / square root seconds
        
        super.init(withMass: mass, radius: radius, gm: gm)
    }
}


The satellite is located at an altitude above a specific point on the surface, defined by latitude and longitude. As a result, I have created a class that contains the three values to define the position of the satellite. I called it distance, but Position is probably better.


class Distance
{
    var latitude: Double// 0-90 degrees, + North
    var longitude: Double// 0-180 degrees, + East
    var altitude: Double// Height above mean body surface, in meters
    
    init(withLatitude: Double, longitude: Double, altitude: Double)
    {
        self.latitude = withLatitude
        self.longitude = longitude
        self.altitude = altitude
    }
}


Notice that latitude and longitude are in degrees and altitude in meters. More about units later.


Velocity also requires a number of parameters to describe, so it is also a class:

 

class Velocity
{
    var horizontal: Double// Speed normal to the radius, in meters/sec
    var vertical: Double// Speed tangential to the radius, in meters/sec
    var speed: Double// Total speed, in meters/sec
    
    init(withHorizontal: Double, vertical: Double)
    {
        self.horizontal = withHorizontal
        self.vertical = vertical
        self.speed = pow(self.horizontal * self.horizontal + self.vertical * self.vertical, 0.5) // Apply Pythagoreus
    }
}


Speed uses meters/second as the units.

 

An orbit requires a satellite:

 

class Satellite
{
    let initialRadius: Double// Distance from the center of the planet (Body), in meters
    let initialVelocity: Double// Speed of spacecraft at initial point, in meters/second
    let initialZenithAngle: Double// The angle between the position and velocity vectors, in radians
    
    init(withRadius: Double, velocity: Double, zenithAngle: Double)
    {
        self.initialRadius = withRadius
        self.initialVelocity = velocity
        self.initialZenithAngle = zenithAngle
    }
}

Finally, we need a class to actually do some calculations:

 

/**
 * http://www.braeunig.us/space/orbmech.htm
 */
class OrbitalParameters
{
    let planet: Body
    let satellite: Satellite
    
    init(withPlanet: Body, satellite: Satellite)
    {
        self.planet = withPlanet
        self.satellite = satellite
    }
    
    func computePeriapsisAndApoapsisWith(radius: Distance,
                                       velocity: Velocity,
                                         zenith: Double) -> (periapsis: Double, apoapsis: Double)
    {
        let r1 = planet.radius + radius.altitude
        let sinZenith = sin(zenith)
        let C = (2.0 * planet.gm) / (r1 * velocity.speed * velocity.speed)
        let squareRoot = pow(C * C - 4 * (1 - C) * -(sinZenith * sinZenith), 0.5)
        let root1 = ((-C - squareRoot) / (2 * (1.0 - C))) * r1
        let root2 = ((-C + squareRoot) / (2 * (1.0 - C))) * r1
        var periapsis = 0.0
        var apoapsis = 0.0
        if root1 < root2 {
            periapsis = root1
            apoapsis = root2
        } else {
            periapsis = root2
            apoapsis = root1
        }
        
        return (periapsis, apoapsis)
    }
}
 

We first initialize the object with a satellite and planet. Once initialized, computePeriapsisAndApoapsisWith(radius: velocity: zenith) return periapsis and apoapsis. This points out one of the advantages of Swift, in that it can return a tuple, not just a single parameter. 


The basic equation from Braeunig’s paper is a quadratic equation with two roots:

 



In the above method, I solve for Rp and solve for the two roots. First off, since the actual radius is the sum of the planet’s radius and the altitude of the satellite, I calculate the total radius:

 

        let r1 = planet.radius + radius.altitude

 


Since the sine of the zenith angle is used in multiple places, I compute it once and use the result in multiple places:

        let sinZenith = sin(zenith)


C is a term used in the equation that’s defined as shown below:


 

        let C = (2.0 * planet.gm) / (r1 * velocity.speed * velocity.speed)

 


I calculate the square root term separately for the sake of performance:

        let squareRoot = pow(C * C – 4 * (1 – C) * -(sinZenith * sinZenith), 0.5)


Then I calculate the two roots:
        let root1 = ((-C - squareRoot) / (2 * (1.0 - C))) * r1
        let root2 = ((-C + squareRoot) / (2 * (1.0 - C))) * r1

Once the roots are solved, I need to determine which is the periapsis and which is the apoapsis:
        if root1 < root2 {
            periapsis = root1
            apoapsis = root2
        } else {
            periapsis = root2
            apoapsis = root1
        }

And Swift allows returning the tuple: 

        return (periapsis, apoapsis)

 

Finally, I have a unit test. This is not a full set, but it shows that things work.

    func computeCWith(gm: Double, radius: Double, velocity: Double) -> Double
    {
        let c = (2 * gm) / (radius * velocity * velocity)
        return c
    }
 
    func testZeroVelocity() {
        let altitude = 200000.0// meters
        let speed = 30000000.0 / 3600.0// meters/second
        let initialDistance = Distance(withLatitude: 0.0, longitude: 0.0, altitude: altitude)
        let initialVelocity = Velocity(withHorizontal: speed, vertical: 0.0)
        
        let earth = Earth()
        let satellite = Satellite(withRadius: altitude, velocity: speed, zenithAngle: M_PI / 2.0)
        
        let orbit = OrbitalParameters(withPlanet: earth, satellite: satellite)
        
        var periapsis = 0.0
        var apoapsis = 0.0
        let zenithAngle = M_PI / 2.0
        (periapsis, apoapsis) = orbit.computePeriapsisAndApoapsisWith(initialDistance,
                                                                      velocity: initialVelocity,
                                                                        zenith: zenithAngle)
        
        let initialRadius = altitude + earth.radius
        let C = computeCWith(earth.gm, radius: initialRadius, velocity: speed)
        let oneMinusC = 1.0 - C
        let denominator = 2 * oneMinusC
        let sinZenith = sin(zenithAngle)
        let root = pow(C * C - 4 * oneMinusC * -(sinZenith * sinZenith), 0.5)
        let expectedPeriapsis = ((-C + root) / denominator) * initialRadius
        XCTAssertEqual(expectedPeriapsis, periapsis)
        let expectedApoapsis = ((-C - root) / denominator) * initialRadius
        XCTAssertEqual(expectedApoapsis, apoapsis)
    }

Where from here

Units are a problem. If the calling function and the called function don’t agree on units, bad things will happen. I have some ideas for classes that will handle the units, and enforce unit cohesion. I will have more on this later as I develop them. I have a basic design for C++ template classes, but I haven’t figured out the Swift classes, yet.
 
Before going to that step, I want to refactor the source code into the directory structure expected by the open source Swift compiler on Linux. I could go on to build a Continuous Integration (CI) tool in Swift, to run on OS X and Linux. I’d try to base it on convention and keep things simple. Initially everything would be hard coded. Then I need to build tools to demonstrate using the library. Something like an iOS app that uses the framework. Could also have a version that performs the calculation on the Linux server and publish results in a web service interface. 
 
After that, there are a lot of additions to the framework. This is just a simple beginning.



Posted in Swift | Tagged , , , | Leave a comment

What’s wrong with Comey’s response

Well worth listening to.

https://www.youtube.com/watch?v=bC1Mc6-RDyQ

Posted in Politics | Leave a comment

How we won our liberty

If it sounds familiar, then so be it.

Posted in Politics, Self Defense | Leave a comment

Ending the Apollo Cargo Cult

A good look at how we can settle the Solar System, including Mars, without the expense of SLS/Orion. Also looks at the real history of Apollo.

Posted in New Space | Leave a comment

Swarm technology

https://www.youtube.com/watch?v=AyguXoum3rk

Posted in Engineering | Leave a comment

Cute machines

https://www.youtube.com/watch?v=tf7IEVTDjng

Posted in Engineering | Leave a comment

Don’t get close to the Clintons

This list illustrates why.

Posted in Politics | Leave a comment

Armed Citizens stop mass murderers

A list of instances where armed citizens stopped mass murderers.

Posted in Self Defense | Leave a comment

New Rocket Facility

Blue Origin breaks ground.

Posted in Commercial Space | Leave a comment

17 Equations that changed the world.

Explained here.

Posted in Science | Leave a comment