Celestial Mechanics

Seventh Annual Physics Workshop

Victor Irby

Physics Department

University of South Alabama

Celestial Mechanics

(The following Maple command "resets" or clears all variables).

> restart:

The Ellipse

> plot([4*(1-.5)*(1+.5)*cos(theta)/(1+.5*cos(theta)) + 2, 4*(1-.5)*(1+.5)*sin(theta)/(1+.5*cos(theta)),theta=0..2*Pi], scaling=constrained);

The above drawing is that of an ellipse. The semi-major axes (denoted by "a") has a length of 4 and the semi-minor axes (denoted by "b") a length of about 3.5. The eccentricity of an ellipse is related to "a" and "b" by the relation:

> eccentricity := sqrt(1-(b/a)^2);

Since a=4 and the eccentricity is 0.5, the exact length of the semi-minor axes is given by;

> b := (4)*sqrt((1-((.5)^2)));

We can re-draw our ellipse such that one of the focus points is at our origin (note each focus is at distance of "e times a" (e is eccentricity));

> plot([4*(1-.5)*(1+.5)*cos(theta)/(1+.5*cos(theta)), 4*(1-.5)*(1+.5)*sin(theta)/(1+.5*cos(theta)),theta=0..2*Pi], scaling=constrained);

With our origin defined as above, the "perihelion" distance is equal to 2. What is the "aphelion" distance? (answer = 6).

What happens when we change the eccentricity?

> epsilon := .9;

> plot([4*(1-epsilon)*(1+epsilon)*cos(theta)/(1+epsilon*cos(theta)), 4*(1-epsilon)*(1+epsilon)*sin(theta)/(1+epsilon*cos(theta)),theta=0..2*Pi], scaling=constrained);

> epsilon := 0;

> plot([4*(1-epsilon)*(1+epsilon)*cos(theta)/(1+epsilon*cos(theta)), 4*(1-epsilon)*(1+epsilon)*sin(theta)/(1+epsilon*cos(theta)),theta=0..2*Pi], scaling=constrained);

> epsilon := 'epsilon';

(Note: the last Maple command clears the value of epsilon).

>

Orbital Equation

A planet's orbital distance from the Sun can be expressed using the equation for an ellipse. In polar coordinates:

> r := A/(1+epsilon*cos(theta));

> A := a*(1-epsilon)*(1+epsilon);

> A := ra*(1-epsilon);

>

where "ra" is the aphelion distance from the Sun and " e " is the orbital eccentricity. The first equation, listed above, can also be used when given the semi-major axes "a".

As an example, the orbital eccentricity of the planet Mercury is 0.206, with an aphelion distance of 0.467 A.U. (1 A.U. is equal to 150 billion meters). Thus, the orbital equation for Mercury is:

> rmercury := .371/(1+.206*cos(theta));

A plot of Mercury's orbit about the Sun, (with the Sun situated at the origin), is given below.

> plot([rmercury*cos(theta), rmercury*sin(theta), theta=0..2*Pi], scaling=constrained);

Compare Mercury's orbit to a "circular" orbit centered at the Sun.

> plot({[rmercury*cos(theta), rmercury*sin(theta), theta=0..2*Pi], [0.3076*cos(theta), 0.3076*sin(theta), theta=0..2*Pi]}, scaling=constrained);

The following sequence of commands will plot the orbits of the inner planets.

> rvenus := .723/(1+.007*cos(theta));

>

> rearth := 1/(1+.017*cos(theta));

>

> rmars := 1.51/(1+.093*cos(theta));

>

While we are at it, let us plot the orbits of ALL the planets.

> rjupiter := 5.2/(1+.048*cos(theta));

>

> rsaturn := 9.5/(1+.056*cos(theta));

>

> ruranis := 19.2/(1+.047*cos(theta));

>

> rneptune := 30.1/(1+.009*cos(theta));

>

> rpluto := 39.5/(1+.247*cos(theta));

>

> plot({[rmercury*cos(theta), rmercury*sin(theta), theta=0..2*Pi], [rvenus*cos(theta), rvenus*sin(theta), theta=0..2*Pi], [rearth*cos(theta), rearth*sin(theta), theta=0..2*Pi], [rmars*cos(theta), rmars*sin(theta), theta=0..2*Pi], [rjupiter*cos(theta), rjupiter*sin(theta), theta=0..2*Pi], [rsaturn*cos(theta), rsaturn*sin(theta), theta=0..2*Pi], [ruranis*cos(theta), ruranis*sin(theta), theta=0..2*Pi], [rneptune*cos(theta), rneptune*sin(theta), theta=0..2*Pi], [rpluto*cos(theta), rpluto*sin(theta), theta=0..2*Pi]}, scaling=constrained);

>

Binary Star System

The following commands will create an "animation" of two stars orbiting each other. The masses of the stars are "m1" and "m2" respectively. Using Maple's programming capabilities, we can change individual masses of the stars.

> m1 := 50;

> m2 := 50;

> r1 := 10;

> r2 := m2*r1/m1;

> x1 := -r1*cos(theta);

> x2 := r2*cos(theta);

> y1 := -r1*sin(theta);

> y2 := r2*sin(theta);

> size := 1;

> with(plots):

Warning, the name changecoords has been redefined

> animate({[size*cos(z) + x1, size*sin(z) + y1, z=-Pi..Pi], [size*cos(z)+x2, size*sin(z)+y2, z=-Pi..Pi]}, theta=0..2*Pi, frames=75, scaling=constrained);

> m1 := 500;

> m2 := 50;

> r1 := 10;

> r2 := m2*r1/m1;

> x1 := -r1*cos(theta);

> x2 := r2*cos(theta);

> y1 := -r1*sin(theta);

> y2 := r2*sin(theta);

> size := 1;

> with(plots):

> animate({[size*cos(z) + x1, size*sin(z) + y1, z=-Pi..Pi], [size*cos(z)+x2, size*sin(z)+y2, z=-Pi..Pi]}, theta=0..2*Pi, frames=75, scaling=constrained);

>

Orbital Transfer

A 6000 kg planetary probe is in a circular solar orbit with a radius of 108 million kilometers (orbital radius of Venus). The space agency wants to put the probe in the same orbit as Earth, a circular orbit with a radius of 150 million kilometers. Can you think of a way to do this? Here are the two orbits:

>

> plot({[150*cos(theta), 150*sin(theta), theta=0..2*Pi], [108*cos(theta), 108*sin(theta), theta=0..2*Pi]},scaling=constrained);

Maybe the following animation can give you an idea.

>

> rp := 108;

> ra := 150;

>

> r := A/(1+epsilon*cos(theta*t/10));

> epsilon := (1-rp/ra)/(1+rp/ra);

> A := rp*(1+epsilon);

>

> animate({[150*cos(theta), 150*sin(theta), theta=0..2*Pi], [108*cos(theta), 108*sin(theta), theta=0..2*Pi], [r*cos(theta*t/10), r*sin(theta*t/10), theta=0..2*Pi]}, t=0..10, frames=75, color=blue, scaling=constrained);

The first step is to increase the probe's speed so that the probe's orbit is changed into an elliptical orbit, with a "perihelion" distance equal to the radius of its initial orbit and an "aphelion" distance equal to the radius of the final orbit. What is the necessary change in speed to place the probe into an elliptical orbit?

The probe's initial speed can be obtained by setting the gravitational force equal to the probe's mass times its centripetal acceleration and solving for "v" (R is the probe's distance to the Sun)

> vi := sqrt(G*M/R);

>

The total energy of the probe is always the sum of its kinetic energy and potential energy

> E := m*v^2/2 -GMm/R;

For a circular orbit the total energy is

> Ec := -G*M*m/(2*R);

For elliptical orbits, the total energy is (a is the semi-major axes length)

> Ee := -G*M*m/(2*a);

If we can change the probe's speed quickly , the potential energy of the probe will be the same after increasing its speed. The change in the probe's kinetic energy can then be obtained from:

> KEchange := Ee-Ec;

( I don't know why Maple likes to write it this way ).

Thus, in order to change the probe's orbit from the initial circular one, into an elliptical one, the final kinetic energy is obtained by adding the "change" in kinetic energy to its "initial" kinetic energy:

> KEfinal := m*(vi^2)/2 + KEchange;

We can now assign appropriate values.

> G := 6.673*10^(-11); M := 1.99*10^(30); m := 6000; R := 108*10^9;

> KEfinal;

We need one more quantity, "a", the semimajor axes. One can easily show that "a" is equal to the sum of the perihelion and aphelion distances divided by 2..

> a := (108*(10^9) + 150*(10^9))/2;

Thus, the final kinetic energy is

> KEfinal;

and the final speed of the probe is

> vfinal := sqrt(2*KEfinal/m);

Therefore, the probe's initial orbital speed needs to be increased by (delta v)

> delv := vfinal - vi;

This corresponds to an increase in orbital speed of 6144 mph! (a 4.67 minute rocket burn for an acceleration equal to "g").

However, one must not forget the next step . Once the probe reaches Earth's orbit at its "aphelion", what must then be done? I will let you figure it out. Here is a final animation.

>