Calculating orbital velocity with Python

Earlier we calculated escape velocity for a rocket using Python, and today, we are going to look at calculating orbital velocity using Python. For a background on some of the terms I use, you may want to refer to the above linked article.

Orbital velocity is what the Space
Shuttle needs to reach to get to orbit.

Orbital velocity is very similar to escape velocity, except that it is instead what is necessary to achieve orbit, rather than totally escape a planet's gravity well. That is, if we're trying to get into orbit from the ground (like the Space Shuttle), we calculate orbital velocity. If we're trying to get to Mars (i.e. escape Earth's gravity well), we calculate escape velocity (as well as Mars' escape velocity in that case).

The formula for orbital velocity is as follows:

Vo ≈ sqrt[ GM / r ]

where

Vo = Orbital velocity (m/s)
G = Gravitational constant
M = Mass of planet (kg)
r = Distance from planet's center of gravity the rocket is launching at (m)

We can easily write a Python function to solve for this. In fact, it's almost identical to the one we wrote for calculating escape velocity:

from decimal import *

_GRAVITATIONAL_CONSTANT = Decimal('6.67384E-11')

def orbital_velocity(mass, radius):
    '''Calculates the delta-V orbital velocity.

    Returns delta-V in meters per second for the orbital
    velocity for a given object given its mass in kilograms
    and its radius in meters.'''

    if mass == Decimal('0') or radius == Decimal('0'):
        # User gave a mass or radius which is unrealistic,
        # abort.
        return None
    
    result = Decimal((_GRAVITATIONAL_CONSTANT * mass) / radius)
    result = result.sqrt()

    return result

Let's run this using our figures for Earth (mean radius of 6,371 km, mass of 5.9736E24 kg), assuming the above code is in sci.py:

>>> from decimal import *
>>> import sci
>>> sci.orbital_velocity(Decimal('5.9736E24'), Decimal('6371000'))
Decimal('7910.467706232366495042659538')

So, 7.9 km/s is approximately the orbital velocity according to our function. And, if we check the Wikipedia article for orbital velocity, it mentions that for low Earth orbit, the orbital velocity ranges from 6.9 km/s to 7.8 km/s! That's what I call close enough for government work.

Calculating escape velocity with Python

As a kid I was fascinated by rockets, the ideals of space exploration, and the promise of a brighter tomorrow. While it looks like our dreams are on hold for now with the evisceration of NASA, we can still have some fun by running through some of the necessary calculations ourselves using Python, a free scripting language.

For this, I'm assuming you are already familiar with Python. If not, you will probably want to look at the Python 2.7.x tutorial.

Today we'll be calculating escape velocity. Escape velocity is the speed which is necessary to escape the influence of gravity from a planet or other heavenly body. That is, if a rocket wants to launch from Earth and then leave Earth's orbit, it has to have enough fuel to reach escape velocity.

The equation for escape velocity is as follows:

Ve = sqrt[ (2 * G * M) / R ]

Ve = Escape velocity (m/s)
G = Gravitational constant (N)
M = Mass of the planet/heavenly body (kg)
R = Distance from the planet/heavenly body's center of gravity we're interested in solving for (m)

Note that R is a little funky. If we're launching at sea level, we would just use the radius of the planet as R, but if we're higher up, like on a mountain, the escape velocity will be less.

How can we implement this in Python? To increase re-usability, I decided to make a function to solve for escape velocity, given a planet's mass and radius:

from decimal import *

GRAVITATIONAL_CONSTANT = Decimal('6.67384E-11')

def escape_velocity(mass, radius):
    '''Calculates the delta-V escape velocity.

    Returns delta-V in meters per second for the escape
    velocity for a given object given its mass in kilograms
    and its radius in meters.'''

    if mass == Decimal('0') or radius == Decimal('0'):
        # User gave a mass or radius which is unrealistic,
        # abort.
        return None
    
    result = Decimal((2 * GRAVITATIONAL_CONSTANT * mass) / radius)
    result = result.sqrt()

    return result

Let's work through this piece by piece.

I decided to use the decimal module, as I wanted greater precision than the built-in float was going to give me. Binary floating point numbers can't represent decimals exactly, so I decided that since we're dealing with such large and small numbers here, we had better use a unit type that's designed for that type of work.

On the next line, I define the gravitational constant, using the Decimal type.

Now we get to the actual function. For best results, the user should pass in Decimals for both mass and radius.

Python then checks to see if either mass or radius are zero, which is highly unrealistic, so it aborts. Otherwise it continues.

The function then calculates the result, and returns it as a Decimal.

How does this look in practice? Let's say we were in the Python interpreter, and the previous code was in sci.py in the local directory. We'll calculate the escape velocity for Earth, which according to Wikipedia has a mass of 5.9736E24 kg, a mean radius of 6,371 km (or 6,371,000 m), and an escape velocity of 11.186 km/s (or 11,186 m/s):

>>> import sci
>>> from decimal import *
>>> sci.escape_velocity(Decimal('5.9736E24'), Decimal('6371000'))
Decimal('11187.09071486820095554468233')

As we can see, it is calculating an escape velocity of 11,187 m/s. Variation could be introduced depending on what radius we use for Earth; I used mean radius as opposed to the equatorial radius or polar radius.

Let's try another celestial object, the Moon:

>>> import sci
>>> from decimal import *
>>> sci.escape_velocity(Decimal('7.3477E22'), Decimal('1737100'))
Decimal('2376.108079493511717080823539')

This is based off a mass of 7.3477E22 kg, and a mean radius of 1,737.10 km. Wikipedia gives an answer of 2.38 km/s (or 2,380 m/s), which is again close to our calculation.