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.