Three-Dimensional Modeling with Application
to Long-Range Artillery -Satellite
and Planet Orbits
The path of a projectile that returns
to Earth is thought of as a portion of an orbit, a sub-orbit. At a sufficiently
high altitude and given appropriate initial velocity, the projectile path can have
the appearance of a spiraling near elliptic shape that may circumnavigate, orbit,
the globe many times before returning to Earth.
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:
A plot of the object's altitude versus time is seen
The Earth is Not Spherical
The World Geodetic System,
WGS, describes Earth in terms
of a Coordinate Frame, a Reference Ellipsoid for its surface and a Geoid, (gravitational
equipotential surface), that defines nominal sea level.
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
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
Point-to-point distances in a city are adequately represented
as straight-line distances as measured by a surveyor's chain, a flat Earth is presumed.
Over the horizon distances generally employ
great circle distance.
The usual assumption is that the two points are on the surface of a spherical earth
at the same altitude. Intercontinental air flights will usually employ
a great circle route. Given the coordinates of two points, latitude, longitude
and their common altitude it is a relatively simple matter to compute the length
of the great circle arc that joins the points. This is the distance model that is employed
Distances between points are determined in a more complicated manner on the
A Rotating Spherical Earth
Previous chapters have shown that when dealing with the trajectory in a short-range
situation, such as with a launched grenade or an ordinary artillery shell, the assumption of a flat Earth with uniform gravity field produces acceptable results. For long-range situations, a more precise treatment for the trajectory is achieved by accounting for a spherical Earth that rotates on its axis.
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
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
The following figure reproduces some of the essential elements of the reference.
This figure shows a vector r
in a coordinate system with centre
with Cartesian components x, y, z and spherical coordinates r, θ, φ. The r in the figure represents the distance from the origin, which happens to be
the magnitude of the vector (i.e. the scalar distance from O), while θ and φ are
angles measured from the z and x axes as shown.
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
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 V
at position r
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 u
φ based at position r
. They show the direction
for each of the velocity vector's spherical components: u
radially outward, u
θ and u
φ lie in a plane which
is perpendicular to r
θ pointing in the
θ direction and u
φ pointing in the φ direction.
If we imagine a
sphere with centre at O passing through the point r
, we can see that
the plane formed by u
θ and u
φ is tangent to the
sphere. We will refer to this plane as the tangent plane or, from the viewpoint
of an observer at r
in orbit around O, the horizon plane.
velocity components Vr, Vθ, Vφ are the projection of the velocity V
onto the directions defined by u
It is interesting to note that for a given V
the components Vx,
Vy, Vz are invariant with position whereas the components Vr, Vθ, Vφ vary with position
because the unit vectors u
change direction depending on the position of r
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.
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 u
r, can be projected
onto the x,y plane as Vr * sin(θ) and this result projected onto the x direction
as Vr * sin(θ) * cos(φ). So the contribution to Vx by Vr is Vr * sin(θ) * cos(φ).
Similarly, for Vθ we note that it is parallel to u
θ which is contained
in the z,
r plane and points downward by an amount θ. Hence, its
projection onto the x,y plane is Vθ * cos(θ) and this result projected onto the
x direction is Vθ * cos(θ) * cos(φ). For Vφ we note that it is parallel to u
which is contained in the x,y plane and points in the negative x direction by an
amount φ. Hence, its projection onto the x direction is -Vφ * sin(φ). Summing the contributions
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(θ)
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 r
= (x, y, z). From basic trigonometry,
the definitions of sine and cosine allow us to write:
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
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
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
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.
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)
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
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.
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
φ direction. The velocity components
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 /
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
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 = - ∆t * [ G * M * (1 - ph) * x / r
^ 3 - Fdx ] recall that Fd is the
∆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
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.
provides the fundamental relationship (see the section "Keplerian Orbits",
paragraphs and sentences are relative to that section):
r = p / (1 + e * cos(θ)) paragraph
= 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(θ))
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:
IntersectArcAngle: 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
semiMajor: For elliptic orbits is (R0 +R1) / 2
OrbitPeriod: Only for elliptic orbits:
= (4 * π^2 * semi-major^3 / G * M)^0.5
escapeVelocity: = (2 * G * M / (earthRadius + launchAlt))^0.5
siderialRadius: Geostationary orbit radius, calculated for the given G * M and siderialDay:
= (siderialDay^2 * G * M / (4 * π^2))^ (1 / 3)
siderialAlt: Geostationary orbit altitude:
= siderialRadius - earthRadius
siderialVelocity: Geostationary tangent velocity:
= siderialRadius * 2 * π / siderialDay
Keplerian Orbital Elements
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
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
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
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
= acos(hz / h) in the range 0 to 180 degrees, acos() is the inverse
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:
= 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.
A spreadsheet macro for the three-dimensional model that takes into account rotation,
gravitational attraction, atmospheric resistance, and buoyancy will be described.