For an interesting description of how a real satellite was placed into a geostationary orbit after a launch from India look at this article and video.
Some tasks that must be completed analytically prior to working with the code include:
- Determine the speed \(v_1\) required for the satellite to orbit in a circle in low-Earth orbit [LEO]. To determine this speed you should apply the momentum principle to the satellite under the assumption it is orbiting in a circle of radius \(R_1\).
This speed will be used to set the initial velocity of the satellite in the program. You will start the satellite on the \(-x\) axis so the initial position will be \(\langle -R_1, 0, 0\rangle\) and the initial velocity will be \(\langle 0, -v_1, 0\rangle\).
- Determine the radius \(R_2\) and speed \(v_2\) for a satellite in geostationary orbit [GSO]. Again, use the momentum principle to determine these values. You will need to require that the orbit of the satellite be 1 day.You will use this speed to place the satellite into GSO following the elliptical transfer orbit trajectory (described in 3 below).
This insertion will occur when the satellite crosses the positive \(x\) axis. The position and velocity of the satellite immediately following the GSO insertion “burn” will be \(\langle R_2, 0, 0\rangle\) and \(\langle 0, v_2, 0\rangle\).
- Determine the speed required at perigee \(v_p\) for the geostationary transfer orbit (GTO). To determine this speed you will need to apply the principles of conservation of energy and conservation of angular momentum to the Earth-satellite system. These two equations will have two unknowns, the speed at perigee \(v_p\), and the speed at apogee \(v_a\).
The GTO insertion burn will occur when the satellite passes the negative \(x\) axis, after completing the fourth circular low-Earth orbit. See below for more instructions on how to determine when to “burn”, i.e. when to increase the velocity of the satellite.
Use the following starter code and instructions from class to complete the project. The program will draw the Earth in a polar view, looking down on the North pole and will draw some semi-transparent rings at the LEO and GSO radii of \(R_1\) and \(R_2\). These rings can be used to help ensure that your program is operating properly.
scene = canvas(width=800, height=800)
########################################################################
# #
# Geostationary Transfer Orbit #
# #
# The purpose of this program is to show how to use a Hohmann transfer #
# to deliver a satellite from Low Earth Orbit (LEO) to a geostationary #
# orbit over the equator. #
# #
########################################################################
# Define required constant values
G = 6.674e-11 # Newton's Gravitational Constant [N m**2/kg**2]
Me = 5.974e24 # Mass of Earth [kg]
Re = 6.387e6 # Radius of Earth [m]
day = 24*60*60 # Number of seconds in one day
# Draw the Earth as a sphere with the proper radius and rotate it so that we
# have a polar view, looking down at the North pole. The center of the Earth
# is at the origin, the North pole is in the +z direction and the equator is
# in the x-y plane.
earth = sphere(radius=Re, texture=textures.earth)
earth.rotate(angle=pi/2, axis=vector(1,0,0))
earth.m = Me
# Define an altitude for Low Earth Orbit (LEO), an
h = 400e3 # Altitude of LEO above the Earth's surface [m]
R1 = Re + h # Radius of LEO.
########################################################################
# #
# 1. Compute the speed of the satellite required for a circular orbit #
# of radius R1. This is the speed in Low Earth Orbit. Store the #
# result in a variable V1 that is a vector of the form #
# #
# V1 = <0, -(LEO speed), 0> #
# #
# so that the satellite will move initially in the negative y #
# direction. #
# #
########################################################################
V1 = vector(0, 0, 0)
# Define an object for the satellite. Make its radius ridiculously large so
# that we can see it. Start the satellite at <-R1, 0, 0> with a velocity of
# V1 = <0, -(LEO Speed), 0>. If V1 is correct the satellite should move in a circle.
satellite = sphere(radius=5e5,make_trail=True, color = color.red, trail_radius=1e5, interval=10, retain=200)
satellite.m = 1500
satellite.pos = vector(-R1,0,0)
satellite.v = V1
satellite.p = satellite.m*satellite.v
# Draw a ring showing a corridor where the satellite should go when in LEO
LEO = ring(axis = vector(0,0,1), radius = R1, thickness = h/2, opacity = 0.5)
########################################################################
# #
# 2. Compute the radius and speed for a satellite that is in a #
# geostationary orbit. Store the radius in the variable R2 and the #
# speed in a vector of the form #
# #
# V2 = <0, +(GSO speed), 0> #
# #
# so that the satellite will initially move in the +y-direction at #
# the apogee kick following the elliptical geostationary transfer. #
# #
########################################################################
R2 = 3e7
V2 = vector(0, 0, 0)
# Draw a ring showing a corridor where the satellie should go when in GSO.
GSO = ring(axis = vector(0,0,1), radius = R2, thickness = h/2, opacity = 0.5)
########################################################################
# #
# 3. Compute the initial velocity of the satellite in the elliptical #
# geostationary TRANSFER orbit (GTO). This is the velocity of #
# perigee (the closest distance to Earth). #
# #
########################################################################
Vp = vector(0, 0, 0)
t = 0
dt = 60 # Try a time step of one minute = 60 seconds for now.
while(True):
rate(100)
# Show the Earth's rotation
earth.rotate(angle=2*pi*dt/day, axis=vector(0,0,1))
# Save the satellite position for computations and orbit determination
r = satellite.pos
# Compute the gravitational force on the satellite
# Update the momentum, velocity, and position of the satellite.
t = t + dt
Finding When/Where to “Burn”
Here are some hints to help determine how to transfer from one orbit to the other. It is easier to decide when to burn (when to increase the speed) by considering the position of the satellite instead of an amount of elapsed time. You will need to do two transfers:
- LEO -> GTO: You will need to transfer from a circular low-Earth orbit at radius \(R_1\) to an elliptical orbit that has it’s perigee at \(R_1\) and apogee at \(R_2\). Here you will need to change the speed of the satellite from \(v_1\) to \(v_p\). We will do this burn when the satellite crosses the negative \(x\) axis from positive \(y\) to negative \(y\).
- GTO -> GSO: This is a transfer from the elliptical orbit to a circular geostationary orbit. You will need to increase the speed from \(v_a\) to \(v_2\). You should make this “burn” first pass of the satellite across the positive \(x\) axis after it is in the elliptical orbit. In our geometry this will be a crossing of the positive \(x\) axis.
Notice that in the computational loop the position of the satellite before the update is stored in the variable r. After the position update the satellite position is in the variable satellite.pos So, it should be possible to test to see if the satellite has crossed the negative x-axis using code similar to what is shown below where we test for the previous \(y\) coordinate being positive and the current \(y\) coordinate is negative.
# Search for crossings of the y-axis by comparing our previous y-coordinate
# (r.y) to our latest y-coordinate (satellite.pos.y).
# Crossings of the y-axis from the positive side to the negative side will
# represent the perigee of our Geostationary Transfer Orbit (GTO).
#
# This is the point where we need to give a delta-V to the satellite that will
# change its speed to the speed at perigee for our GTO.
if (r.y>0 and satellite.pos.y<0):
perigee = perigee + 1
# Give the "kick" on the fourth "perigee" crossing.
if (perigee == 4):
satellite.p = satellite.m*Vp
Notice that this code, which should be in your while loop after the momentum and position of the satellite has been updated, also counts how many perigee crossings have occurred. It will perform the burn, that is change the velocity to \(\vec{V}_p\), on the fourth perigee crossing.
For this code to work you will need to initialize the variable perigee to zero prior to the start of the while loop.
perigee = 0
You will need to also add similar code to perform the burn at apogee to insert into the geostationary orbit.
