We recently posted a blog entry celebrating the anniversary of the Apollo 11 landing on the Moon. Now, just a couple weeks later, we are preparing for another first: the European Space Agency’s attempt to orbit and then land on a comet. The Rosetta spacecraft was launched in 2004 with the ultimate goal of orbiting and landing on comet 67P/Churyumov–Gerasimenko. Since the launch, Rosetta has already flown by asteroid Steins, in 2008, and asteroid 21 Lutetia, in 2010.

NASA and the European Space Agency (ESA) have a long history of sending probes to other solar system bodies that then orbit those bodies. The bodies have usually been nice, well-behaved, and spherical, making orbital calculations a fairly standard thing. But, as Rosetta recently started to approach comet 67P, we began to get our first views of this alien world. And it is far from spherical.

Credits: ESA/Rosetta/MPS for OSIRIS Team MPS/UPD/LAM/IAA/SSO/INTA/UPM/DASP/IDA

The comet’s nucleus can be vaguely described like a flat platform with a spheroidal rock attached to one end, with the entire system rotating in a complex manner. How do you orbit such an object? The exact dynamics to model are beyond the scope of this post, but we can make use of Wolfram|Alpha and the Wolfram Language, including new features in *Mathematica* 10, to do some thought experiments. Can we generate a simple model of the comet, and then try to simulate the probe’s approach and orbit of this irregular body? Keep in mind that comet nuclei are relatively small bodies with small masses. This means that there’s not a lot of gravity to hold you in orbit or to help you land. Initial estimates suggest this comet has an escape velocity of about .5 m/s. If you move any faster than this, you will fly off into space and never come back, so the approach must be very slow, on the order of 10 cm/s, or risk never being captured.

Wolfram|Alpha recently added support for physical systems that we can use to model our comet. First, we need the shape. Let’s model the comet as a massive and flat cuboid with a sphere sitting on one end. We can access the new Wolfram|Alpha functionality using `EntityValue` in the Wolfram Language.

Now we can render the resulting model.

The next step is to obtain the gravitational potential of this system. Let’s assume that 2/3 of the total mass is in the cuboid and 1/3 of the total mass is in the sphere. Let’s also assume that the densities of the sphere and the cuboid are the same.

Now, we can plot a slice of the potential using `ContourPlot`. A preliminary estimate of the comet’s total mass is also included.

We can also visualize the gravitational potential in 3D using `ContourPlot3D`. Here, we’ve chosen to look at a specific iso-potential near the surface of the model.

The force of gravity can be found by differentiating the gravitational potential.

One of the more difficult problems is that in order to orbit and land on the comet, you have to approach the comet and then slow down so that the gravity of the comet can capture you, the orbiter. This typically requires a rocket to be fired to slow the orbiter down. So here we define a new force, created by this rocket, that opposes the forward velocity of the probe.

For visualization purposes, we want to tag the times when the rocket burns and shuts off, so we create a couple of state variables. `NDSolve` can be used to handle all of the acceleration, velocity, and position changes. In addition, we can use `WhenEvent`‘s functionality to handle several possible conditions to watch out for, and take specific actions if those conditions are met.

We can start with one possible scenario where we choose an arbitrary starting point and velocity, but the velocity and burn times result in the orbiter hitting the comet instead of obtaining a capture orbit.

We can see that a single rocket burn takes place.

Red segments show where the rocket was fired to slow down.

With adjustments to the initial velocity and burn distances, we can achieve a more slowly decaying orbit that involves multiple burns.

In this case, there are nine separate rocket burns resulting in orbits that start out highly eccentric, but then settle down into closer, more circular orbits.

We can plot the velocity over the course of the orbital insertion and orbit phases. Velocities are measured in meters/second, and time is in seconds.

The ESA will have a much more difficult problem, likely requiring many separate burns over several months, and this is even before it attempts to drop a lander onto the comet in November. As the rocket approaches the comet, more information becomes available, such as refined mass estimates and the location of potentially hazardous debris from jets issuing from the icy comet. All of this must be carefully considered before refining the approach and orbital maneuvers. The next few months should be really exciting!

To experiment with the parameters of this simulation yourself, try using the deployed version created using `CloudDeploy` in the Wolfram Language.

Download this post as a Computable Document Format (CDF) file.

Awesome! I have been using Mathematica for many years and I will admit that this is well beyond my ability but glad someone is able to not only “do the math” but share it with the world. I will surely show this post to my students this year.

It’s times like these that I miss maths. I had a chance to study it as part of my graduate degree but I ended up dumping it for some other degree. It’s simply amazing how mathematics is involved in our everyday life and in our quest for knowledge of the outside world.

I almost shed a nostalgic mathematic tear here, my friend. I’m also an amateur fan of astronomy so, needless to say, I’m happy to have found this neat treasure of a website as is this one!

Romeo