Hands-On Math

Last Revision March, 2011

Last Revision March, 2011

In the absence of an atmosphere the path of an orbiting body is governed by gravitational attractions and can closely approach the shape of a ellipse or circle and may be treated analytically as such.

Analytic determination of the velocity and position of a body subject to gravitational attraction, which is presumed to be moving along an elliptical path, turns out to be a tough problem that has required either the use of an infinite series or iterative numerical methods for solution.

The approach of numerical modeling with a spreadsheet that is taken in this topic, does not explicitly encounter that problem, does not require the path to be elliptic and the presence of Earth's atmosphere can be included in the model.

For an orbit that is decaying in the atmosphere see next:

The World Geodetic System,

The dawn of long range missiles and satellites in the fifties sparked intense efforts
to refine the WGS that continue today. The now quite familiar GPS, Global Positioning
System, employs WGS 84. There was a major revision to WGS 84 in 1996, *Earth Gravitational
Model 96, or *
*EGM 96**. *

*The precision that the WGS provides does not add to the level of understanding of motion
in orbits and sub orbits that HOM provides and its implementation requirements
are beyond HOM's resources.*

HOM models the Earth as Spherical.

**Distances on a Spherical Earth**

Over the horizon distances generally employ

Distances between points are determined in a more complicated manner on the

For artillery and similar situations, the basic idea is to consider the trajectory of the object as being a three-dimensional sub-orbital path governed by the inverse square law of gravitational attraction. This approach can be easily generalized to cover the motion of a satellite or similar object, such as a planet, revolving around a central mass. A simplifying assumption is made that the central mass is very large in comparison to the satellite. This assumption yields quite accurate results in the practical cases that we investigate.

In this topic we concentrate on projectiles that are in sub-orbit or in orbit around the Earth. The coordinate system will therefore be Earth-centred. For Earth-centred coordinates, the initial conditions and description for the orbit are more easily described and documented in spherical coordinates (r, θ, φ). However, performing numerical calculations of orbits in spherical coordinates requires a division by sin(θ), where θ is the polar angle. This division produces uncertain results when θ approaches 0 or 180 degrees (near the poles). Using Cartesian coordinates (x, y, z) avoids this problem.

The design of the spreadsheet is such that the user may input orbit and other parameters in spherical coordinates as necessary. The results of calculations are also displayed in spherical coordinates where necessary. Internally, a part of the macro that performs orbit calculations must transform between spherical and Cartesian coordinates.

Transforming between spherical and Cartesian coordinates is straightforward but must be done carefully since inverse trigonometric functions can place a result in an incorrect quadrant. This is particularly true for calculation of bearing angle on the surface of the Earth.

A description of spherical coordinates, in .PDF format, can be found here. The following figure reproduces some of the essential elements of the reference.

This figure shows a vector

In our use, the vector could represent the position of a satellite in orbit, with O being the Earth's centre.

The process for transforming between spherical and Cartesian coordinates is based on the following:

From the figure, using basic geometry, it is evident that the projection of the vector onto the x,y plane produces a component of length r * cos(θ-90). Likewise, we may project this result onto the x and y axes and find that the x component is r * cos (θ-90) * cos (φ) and the y component is r * cos (θ-90) * sin(φ). In the same manner the z component is r * cos (θ). Using the fact that cos(θ-90) = sin(θ) we have the fundamental transformation of coordinates r, θ, φ to x, y, z as:

x = r * sin(θ) * cos(φ) (1)

y = r * sin(θ) * sin(φ) (2)

z = r * cos(θ) (3)

To obtain r, θ, φ from x, y, z: sum the square each of (1) (2) and (3) and use the identity sin(φ)^2 + cos(φ)^2 = 1 to obtain:

r^2 = x^2 + y^2 + z^2 (4)

Looking at the figure and using the definition of tangent, we can see that

tan(φ) = y / x or

φ = atan(y / x) (5) where atan is the inverse tan function

Also, tan(θ) = (x^2 + y^2)^0.5 / z or

θ = atan((x^2 + y^2)^0.5 / z) (6)

In physics, it is frequently useful to deal with a type of vector that defines the state or condition of a body at a point in space. In our case such vectors are the forces acting on an orbiting body as well as its velocity and acceleration. These are useful to treat as vectors because they have both magnitude and direction. As an example, consider the following figure where a satellite has velocity

As with position, the velocity can be resolved into 3 Cartesian components Vx, Vy, Vz or spherical components Vr, Vθ, Vφ. Note in the figure the three unit vectors

If we imagine a sphere with centre at O passing through the point

The spherical velocity components Vr, Vθ, Vφ are the projection of the velocity

In many or our studies of projectiles and satellites, we employ the initial velocity in components closely related to Vr, Vθ, Vφ and need to transform these to Vx, Vy, Vz to perform further calculations.

Consider the derivation of Vx given Vr, Vθ, Vφ. Each of Vr, Vθ, Vφ will have a component in the x direction. If each of these components is found, then their sum will yield Vx.

Vr, for example, because it lies along

Vx = Vr * sin(θ) * cos(φ) + Vθ * cos(θ) * cos(φ) - Vφ * sin(φ) (7)

A similar procedure gives the other component transformations as:

Vy = Vr * sin(θ) * sin(φ) + Vθ * cos(θ) * sin(φ) + Vφ * cos(φ) (8)

Vz = Vr * cos(θ) - Vθ * sin(θ) (9)

When we study the orbital motion of projectiles and satellites around the Earth, our model assumes a spherical Earth so massive that the influence of the body in motion on the Earth is negligible. Thus the force of gravity is radial acting from the centre of the Earth. If the origin of the coordinate system is placed at the centre of the Earth, then the force of gravity on the body, based on Newton's theory of gravity, is totally radial and of magnitude:

F = G * M * m / r^2

As has been noted, equations (7) - (9) apply to any vector quantity related to the body and thus they can be re-cast for the force of gravity. However, we can take advantage of some simplifications, since the force is strictly radial; that is, if the force has components Fr, Fθ, Fφ then Fr = F, Fθ = Fφ = 0, then the equations corresponding to (7) - (9) are:

Fx = F * sin(θ) * cos(φ)

Fy = F * sin(θ) * sin(φ)

Fz = F * cos(θ)

The body is at location

cos(θ) = z / r

= z / (x^2 + y^2 + z^2)^0.5

sin(θ) = (1 - cos(θ)^2)^0.5

= 1 - z^2 / (x^2 + y^2 + z^2))^0.5

= (x^2 + y^2) ^0.5 / ( x^2 + y^2 + z^2)^0.5

= (x^2 + y^2)^0.5 / r

Note that the quantity (x^2 + y^2)^0.5 amounts to the projection of r onto the x-y plane. In that case, we also have, by definition:

cos(φ) = x / projection of r on x-y plane = x / (x^2 + y^2)^0.5

sin(φ) = y / projection of r on x-y plane = y / (x^2 + y^2)^0.5

Taking all expressions for sine and cosine as well as the gravitational expression for F into the formulas for Fx, Fy, Fz gives:

Fx = F * sin(θ) * cos(φ)

= F * (x^2 + y^2)^0.5 / r * x / (x^2 + y^2)^0.5

= F * x / r thus:

= G * M * m / r^2 * x / r

= G * M * m * x / r^3

Fy = F * sin(θ) * sin(φ)

= G * M * m / r^2 * (x^2 + y^2)^0.5 / r * y / (x^2 + y^2)^0.5

= F * y / r

= G * M * m / r^2 * y / r

= G * M * m * y / r^3

Fz = F * cos(θ)

= F * z / r

= G * M * m / r^2 * z / r

= G * M * m * z / r^3

Buoyancy is also a radial vector force and can be resolved into components in a like manner.

From these equations we see:

Fi = F * cos(θ) = F * i/r (i = x, y, z).

The factor i/r acts to resolve the force along the direction i and is generally referred to as the direction cosine of the force for the direction i. A vector quantity that is non-radial requires slightly more involved expressions for direction cosine. For example, the force or air resistance is directed along a line that is common to the velocity vector but is opposite in direction and has a magnitude which can be calculated from the velocity. If this force is Fd, and the velocity has magnitude Vmag, then using the fact that it is opposite in direction to the velocity vector allows us to find its components as:

Fdi = Fd * (-Vi) / Vmag where (-Vi) / Vmag acts as direction cosine for direction i, i = x, y, z

To show how the transformations between coordinate systems can be used, we look at the case of simulating the flight of a projectile fired from a cannon. The requirement is to describe the initial conditions of flight in terms of parameters that are easily grasped.

Thus the cannon location is configured with input parameters Latitude, Longitude and Altitude. By convention, the direction of aim is configured by barrel Elevation (the angle above the plane of the horizon) and Azimuth (the angle of the barrel measured from the direction of true north in the plane of the horizon).

The plane of the horizon is tangent to the Earth at the point of the cannon's location and thus the perpendicular to the plane passes through the Earth centre and the plane contains the θ and φ directions. Elevation and Azimuth specify the direction of initial muzzle velocity in terms of the local horizon plane. We specify the initial velocity magnitude as Vm, a characteristic of the cannon and its charge.

The spreadsheet macro must first obtain the initial conditions of location (x, y, z) and velocity (Vx, Vy, Vz) from the above user-configured parameters. The following steps are used to do this.

We assume that the coordinate system origin is at the Earth's centre. The coordinate system is fixed in space (does not rotate with the Earth). For convenience, but without loss of generality, one can assume that the equatorial plane is coincident with the x,y plane and the x,z plane is coincident with the cannon location at the instant the shell leaves the muzzle (at time t = 0).

The initial position in spherical coordinates is obtained using the fact that r = EarthRadius + Altitude. The angle θ is 90 - Latitude and φ = 0 from the alignment of the cannon on the x,z plane. With this initial alignment, the Longitude does not bear on the initial conditions. Its configuration is provided to allow calculation of distance to a target, using spherical trigonometry, which is also configured with a Latitude and Longitude.

Assuming:

Lat = Latitude, the latitude of the location on Earth of the launch device (cannon).

Lon = Longitude, the longitude of the location on Earth of the launch device.

EarthRadius = radius of the Earth, which we most often take as 6371010 meters, and is a configurable parameter.

Alt = Altitude, height of cannon above Earth's surface

The equations (1), (2), (3) above give the initial position of the projectile as:

x = (EarthRadius + Alt) * sin(90 - Lat) * cos(0)

= (EarthRadius + Alt) * cos(Lat)

y = (EarthRadius + Alt) * sin(90 - Lat) * sin(0)

= 0

z = (EarthRadius + Alt) * cos(90 - Lat)

= (EarthRadius + Alt) * sin(Lat)

The initial velocity in Cartesian coordinates is obtained by transforming from the velocity in spherical coordinates. Since the cannon (or other launch device) is fixed to the Earth and rotates with the Earth, the Earth's rotational speed for the given Latitude must be added to the φ component of initial velocity. This augmentation of initial velocity does not apply to the case where a satellite is inserted into an orbit since in that case the launch vehicle is not attached to the Earth.

Note that a northerly Azimuth (from -90 to 90 degrees) corresponds to the negative sense of the θ component.

Assuming:

Elev = Elevation, the angle of the barrel above its local horizon plane

Az = Azimuth, the angle of the barrel from true North.

Vm = muzzle velocity, the initial velocity at t = 0

siderialDay = the time the Earth takes to complete one revolution with respect to the fixed coordinate system

The initial spherical velocity components are obtained by projecting the initial velocity, Vm, onto the spherical coordinates using the Elev and Az with procedures similar to those already described. The following figure is provided to show some of the details in obtaining the spherical components of the velocity:

The augmentation of Vφ due to Earth's rotation is 2 * π * EarthRadius * cos(Lat) / siderialDay. This is the tangent velocity of the launch point and is always in the

Vr = Vm * sin(Elev)

Vθ = - Vm * cos(Elev) * sin(Az)

Vφ = Vm *cos(Elev) * cos(Az) + 2 * π * EarthRadius * cos(Lat) / siderialDay

These can be transformed to Cartesian coordinates using equations (4), (5), (6) .

As the Earth continues to rotate it carries the cannon with it so the spherical angular coordinates of the cannon at any time t are:

θcannon = 90 - Lat

φcannon = 2 * π * EarthRadius * (t / siderialDay)

Atmospheric air resistance to the motion of the projectile is modeled as proportional to the square of the speed of the projectile with respect to the atmosphere. One can assume that the atmosphere rotates with the Earth and its velocity components (relative to the fixed coordinate system) are:

Wr = 0

Wθ = 0

Wφ = cos(Lat) * 2 * π * r / siderialDay

r = radial distance from Earth centre

= earthRadius + altitude

These can be transformed to Cartesian coordinates in the same manner discussed earlier.

Once the initial conditions are determined, integration can begin to calculate the trajectory of the projectile. The motion of the projectile is governed by the forces of gravity, atmospheric buoyancy and atmospheric resistance. The following definitions and equations are used:

r = radial distance of projectile from Earth centre

Vmag = magnitude of the projectile velocity relative to air

= ((Vx - Wx) ^ 2 + (Vy - Wy) ^ 2 + (Vz - Wz) ^ 2 )^0.5

ph = air density at altitude r - earthRadius

= 1.225 at sea level

Fd = air resistance force at velocity Vmag due to drag

= Cd * ph * Vmag^2 / 1.225

Fdi = air resistance force in direction i

= - Fd * (Vi - Wi) / Vmag where i = x, y, z

In the last definition above, the factor (Vi - Wi) / Vmag represents the direction cosine of the force and is used to resolve the air resistance force along Cartesian coordinates (as has been discussed).

The total force on the projectile is obtained by summing the force of gravity, buoyancy and air resistance. Note that gravity and buoyancy act in radial inward (negative) and outward directions, respectively. They are resolved into x,y,z directions by corresponding factors x/r, y/r and z/r.

The effects of buoyancy and air resistance are covered in Chapter 5 and are used here in a three-dimensional setting. Using that reference we see that the components of force due to gravity and buoyancy = G * M * (1 - ph) * i /r^3 for i = x, y, z.

The total force is obtained by adding to these components the force due to air resistance.

From the relation acceleration = force / mass we may take acceleration as = ∆V / ∆t and thus for a small time interval (step size) ∆t the change in velocity in the x direction is ∆Vx :

∆Vx = - ∆t * [ G * M * (1 - ph) * x / r ^ 3 - Fdx ] recall that Fd is the drag force

more generally:

∆Vi = - ∆t * [ G * M * (1 - ph) * i / r ^ 3 - Fdi ] where i = x, y, z

The accumulation of small amounts of velocity, ∆Vi, to form the ongoing velocity and future positions of an object is the fundamental integration process used in HOM's calculators to solve the differential equations for the motion of a mass subject to gravitational acceleration in an atmosphere.

This integration process was first illustrated and demonstrated in the second topic of HOM's Chapter 4.

It is this expression that is implemented in our three-dimensional spreadsheet that is described in the next topic of this Chapter.

When an object in an orbit or sub orbit does not contend with an atmosphere it becomes
possible to form many accurate analytical conclusions about the object's trajectory
and motion. Such conclusions are useful references from which to gauge the
influence of atmosphere. Our 3D spreadsheet provides analytic results for user reference
prior to commencing numerical integration. The following two-dimensional
diagram shows the classic elliptic orbit and fundamental parameters:

See next for variable names that we employ on the spreadsheet.

**V0:** Is the velocity of the object at the maximum altitude. It is
totally tangential.

**V1:** Is the velocity of the object at the minimum altitude.
It is totally tangential.

**R0:** Is the distance from the centre of Earth at maximum altitude.
This point in the orbit is the Apogee.

**R1:** Is the distance from the centre of Earth at minimum altitude.
This point in the orbit is the Perigee.

**Semi-Major:** this distance is (R0 + R1)/2 and usually represented by
letter **a** in most references.

**Alt0**: Is the height above Earth surface at Apogee = R0 - Earth Radius

**Alt1:** Is the height above Earth surface at Perigee = R1 - Earth Radius

**r:** **r** is the vector from Earth centre to the
satellite. Its magnitude is r in the following equations.

**θ:** Is the angle from semi-major axis to the vector **r**.
This θ should
not be confused with the 3 dimensional θ used in polar coordinates.

**e:** Is the eccentricity of the orbit.

**p:** Is the semi-latus-rectum of the ellipse of the orbit.

**C:** is related to angular momentum and eccentricity as C = G * M
* e / H^2

**H, h:** is used for the per unit mass angular momentum of the orbiting object

Assume the Earth's centre is the origin of a two-dimensional coordinate system with
r the distance from the origin to the satellite and θ the angle from the Semi-Major
axis to the satellite.
Celestial Mechanics
provides the fundamental relationship (see the section "Keplerian Orbits",
paragraphs and sentences are relative to that section):

r = p / (1 + e * cos(θ)) paragraph
4/sentence 1

= 1 / (1/p + e / p * cos(θ))

The reference also provides the relationships:

p = **a** * ( 1 - e^2) paragraph 4/sentence 9

e^2 = 1 - H^2 / (**a *** G * M) paragraph 7/sentence 7

Note that the reference
uses lower case h for angular momentum. In our equations we will use upper case
H for momentum. Using the last two equations:

p = **a** * (1 - e^2)

= **a** * (1 - 1 + H^2 / (**a** * G * M)) = H^2 / (G * M) same as
paragraph 9/sentence 4

Thus, the equation for an elliptical orbit has the form

r = 1 / (G * M / H^2 + G * M * e / H^2 * cos(θ))

= 1 / (G * M / H^2 + C * cos(θ))
where

C = G * M * e / H^2

The analytic calculations produce a number of values that are displayed
on the spreadsheet before numeric integration is performed:

I**ntersectArcAngle:** An output for sub orbital trajectories
used to calculate the analytic distance along Earth's surface is:

2 * **ArcAngle
IntersectArcDistance:** For sub orbital trajectories, an analytic distance
along the surface of the Earth is:

earthRadius * IntersectArcAngle

= (4 * π^2 * semi-major^3 / G * M)^0.5

= (siderialDay^2 * G * M / (4 * π^2))^ (1 / 3)

= siderialRadius - earthRadius

= siderialRadius * 2 * π / siderialDay

It is common in the literature to characterize the orbit of an object in terms of its Keplerian orbital elements since they are easy to visualize. There are six Keplerian elements and they describe the orbit in terms of its parameters, of an ellipse or other type of conic section curve.

The first three of the set of elements are the semi major
axis (**a**),
eccentricity (e), and inclination (i). These three elements determine
the shape of the orbit.

The remaining three elements are: argument of periapse (ω), right ascension of ascending node (Ω), and periapse passage (T). These are used to determine the position of the orbit and the position of the object within the orbit at a specific date and time.

Of the six orbital elements, we will look at only the first three that deal with
basic geometry, **a**, e and i. Aspects of orbital position for a given date and time
are beyond our needs for the problems that we wish to explore. If we assume
that the object is acted on by only the central force of gravity, then the three
parameters are constants.

When other forces such as air resistance or atmospheric buoyancy act on the object the parameters are generally time varying. But in orbits, such as satellites in high orbit around the Earth, these forces are small and the variation of orbital elements due to these forces with time is small. These forces are ignored for this discussion.

As has been explained, the integration of the equations of motion occurs in Cartesian coordinates. The six values that determine the objects motion in Cartesian coordinates are x, y, z, Vx, Vy, Vz. These can be used to obtain the Keplerian parameters using the following steps:

Calculate the angular momentum. The angular momentum, h, is the distance multiplied by the tangential velocity of the object in orbit. It is a constant for an object in orbit and has both magnitude and direction. It has hx, hy, hz components and magnitude h:

hx = y * Vz - z * Vy

hy = z * Vx - x * Vz

hz = x * Vy - y * Vx

h = (hx^2 + hy^2 + hz^2)^0.5

Compute the following quantities r, V and Specific Energy (SE):

r = (x^2 + y^2 + z^2)^0.5

V = (Vx^2 + Vy^2 + Vz^2)^0.5

SE = V^2 / 2 - GM / r

The semi major axis (**a**), eccentricity (e), inclination (i), are then computed as follows:

**a** = - GM / ( 2 * SE )

e = (1 - h^2 / (a * GM))^0.5

i = acos(hz / h) in the range 0 to 180 degrees, acos() is the inverse cosine.

The formulae for converting Cartesian coordinates to orbital elements can be found as a free download as this .pdf file.

Consider an object in an elliptic orbit around the Earth. Although acted upon by the force of gravity, there is no sense of any force on the object and therefore no sense of acceleration. Astronauts in orbiting vehicles, for example, feel weightless. Whenever a non-gravitational force acts on the object, acceleration is sensed and the natural orbit is disturbed.

In our modeling of orbits, we consider a total of four forces: force of gravity, air resistance, buoyancy and thrust (due to propulsion rockets). The sum of these forces will induce an acceleration of the object, call it A. The acceleration will have three components Ax, Ay, Az.

At any position in the orbit there is a gravitational acceleration, g, which varies with radial distance, r, from the centre of the Earth and is directed toward the Earth centre. This acceleration will have components gx, gy, gz that can be calculated as:

g = GM / r ^ 2

gx = - g * x / r

gy = - g * y / r

gz = - z / r.

We can define a sensed acceleration, AccSense, as being the difference between the total acceleration less the acceleration due to gravity:

AccSensex = Ax - gx

AccSensey = Ay - gy

AccSensez = Az - gz

It is common to discuss the magnitude of this acceleration in terms of how many gs an object is subjected to when it is accelerated.

AccSenseMag = ((Ax - gx)^2 + (Ay - gy)^2 + (Az - gz)^2)^0.5

Note that in orbit, neglecting air resistance and buoyancy, Ax = gx, etc., and therefore the sensed acceleration is zero.

When there is thrust from propulsion rockets, the number of gs experienced by an object can become large. For NASA's Space Shuttle, the guideline is a limit of 3 g. The Main Engines and Solid Rocket Boosters, used for initial ascent are powerful enough to reach this limit and, at times, require moderation of thrust to avoid damaging the vehicle and its occupants.

After Main Engine Cutoff (MECO), smaller and less powerful engines collectively
called the Orbital Maneuvering System (OMS) are employed for orbital insertion as
well as for rendezvousing with the ISS. The OMS thrusts result in much smaller g
forces and thrust moderation is not required.

Top | Next Topic | Topics |