Math 142 Project 3: Celestial Mechanics
2025-05-02
Chapter 1 Introduction
In this project we will derive ODEs that govern the motion of a small satellite as it orbits the earth, as it orbits the moon, and as it returns to earth. We will consider the influence of both the gravitational force from the Earth and the moon, as well as the effect of drag once the satellite enters the Earth’s atmosphere. Though our discussion is limited to the satellite, earth, and moon, the exact same ideas apply to all objects in the galaxy (so long as you don’t try to account for relativistic effects). We will restrict ourselves to a two-dimensional model of space, but a three-dimensional model would not be hard to include. It would, however, include more notation and more typing.
1.0.1 Notation
In this project we will consider the position of the satellite as well as the velocity of the satellite. Denote the \(x\)- and \(y\)-coordinates of the satellite as \(x\) and \(y\), respectively. We will use \(vx\) and \(vy\) to denote the velocity of the satellite in the \(x\)- and \(y\)- directions.
In this project we will assume that the center of the earth is at position \((0,0)\), and that it does not move. The center of the moon is approximately \(3.844 \times 10^{8}\) meters from the center of the earth. We will assume that the center of the moon is in position \((3.844 \times 10^{8},0)\) and that it does not move.
Figure 1.1: The configuration of the Earth and moon for this project.
Note: In many programming languages, it is common to write 3.844e8 rather than 3.844 * 10^8. Either will work, but R defaults to the former notation, and it is less typing.
Note: In this project we will often make plots with vastly differing scales on the \(x\) and \(y\) axis. For example, on the figure above, the \(x\) axis covers \(4 \times 10^{8}\) units of space, while the \(y\) axis covers only \(2 \times 10^{7}\). The \(x\)-axis represents 20 times as much space as the \(y\)-axis. If we drew those axes to be the same length on the page, then everything would look stretched out in the \(y\)-direction, and our circular planets would look like ellipses. See the example below. To avoid this, we’ll often set the aspect ratio of the figure manually, using asp=. It is totally fine to not set the aspect ratio manually in your figures, you’ll just have to deal with circular planets that look like ellipses and elliptical orbits that look circular.
Figure 1.2: Same figure as 1.1 without the aspect ratio set correctly
If you want, you can use the following function to calculate a proper aspect ratio for you. Copy and paste this line of code, and execute it. Then you will have the function get_asp available to you for later use.
Once get_asp is available, you can set the aspect ratio with the option asp=get_asp(soln). Again, it is totally fine to skip setting the aspect ratio. You’ll just end up with some squished looking plots.
Note: In order to aid visualization, it is nice to include a representation of the moon and the Earth in figures. This is not required, it is just nice. If you wish, copy and run the following code.
add.earth=function(scale=1){
rad_m=6378137*scale
x0=0
f=makeFun(sqrt(rad_m^2-(x-x0)^2)~x)
p=USAFACalc::plotFunFill(f(x)~x,bottom=-f(x)~x,xlim=c(x0-rad_m,x0+rad_m),add=TRUE,col="darkslategray1")
return(p)
}
add.moon=function(scale=1){
rad_m=1737400*scale
x0=384399000
f=makeFun(sqrt(rad_m^2-(x-x0)^2)~x)
p=USAFACalc::plotFunFill(f(x)~x,bottom=-f(x)~x,xlim=c(x0-rad_m,x0+rad_m),add=TRUE,col="gray42")
return(p)
}To add the moon to any static figure, run add.earth(). If you wish to make it appear larger or smaller than its actual size, use scale=0.5 to make the Earth appear half-size, or scale=2 to make it double size. Similarly, use add.moon() to add the moon to any static figure.
1.1 Gravitational Force
The formula for the gravitational force exerted by one object on another is \[F=\frac{m_1m_2 G}{r^2}\] where \(m_1\) and \(m_2\) are the masses of the two objects, \(G\) is the universal gravitational constant and \(r\) is the distance between the two objects.
The formula is above is correct, but it only tells us how strong the force is, it does not tell us in which direction the force acts. The gravitational force in the \(x\) direction acting on object 1 induced by object 2 is \[F_x = \frac{m_1 m_2 G \cdot (x_2 - x_1)}{\left((x_2-x_1)^2+(y_2-y_1)^2\right)^{3/2}}.\] In this formula, \((x_1, y_1)\) is the position of object 1, \((x_2,y_2)\) is the position of object 2. If \(F_x\) is negative, then the force is pushing object 1 to the left, if it is positive then the force is acting to the right.
Similarly, the force acting on object 1 in the \(y\) direction is \[F_y = \frac{m_1 m_2 G \cdot (y_2 - y_1)}{\left((x_2-x_1)^2+(y_2-y_1)^2\right)^{3/2}}.\] We can create an R function that calculates the force in the \(x\) direction experienced by the satellite as a result of gravitational interaction with earth if the satellite is in position \((x,y)\). We just need to fill in the formula for \(F_x\) above with the mass of the earth, the mass of the satellite, and the universal gravitational constant.
The universal gravitational constant is approximately \[ G=6.674 \times 10^{-11} \text{N m}^2\text{kg}^{-2} \] The mass of the earth is approximately \[m_{earth} = 5.972 \times 10^{24}\text{ kg}.\] The mass of the satellite is approximately \[m_{sat} = 6 \times 10^{3}\text{ kg},\] which is typical of a heavy weather satellite. We know that the earth is in position \((0,0)\).
G=6.6743e-11 #Universal gravitational constant in N m^2/kg^2
mass_e = 5.9722e24 #Mass of the earth, in kg
mass_sat=6000 #Mass of satellite, in kg
Fearth_x=makeFun(mass_sat*mass_e*G*(0-x)/((0-x)^2+(0-y)^2)^(3/2)~x&y)1.1.1 Exercises
Write function that calculates the gravitational force acting on the satellite in the \(y\) direction.
Calculate the force in the \(x\) and \(y\) directions acting on the satellite if the satellite is in position \((1 \times 10^{7},0)\). Your results are in units of Newtons.
Calculate the force experienced by the satellite as a result of gravitational interaction with earth if the satellite is in position \((-4 \times 10^{6},8.5 \times 10^{6})\)
This force is acting to the right and down, which both make sense given the position of the satellite relative to earth.
1.2 Newtonian Mechanics
Newton’s second law says that the acceleration experienced by an object is related to the net force acting on the object via \[ma=F.\] If we apply this formula to motion in the \(x\) direction, remember that acceleration is the derivative of velocity, and assume that gravitational interaction with earth is the only force acting on the satellite, then we arrive at the differential equation \[m_{sat} vx'=F_{earth,x},\] where \(m_{sat}\) is the mass of the satellite.
Equivalently, \[vx'=\frac{1}{m_{sat}}F_{earth,x}.\]
Similarly, the acceleration experienced by the satellite in the \(y\) direction is \[vy'=\frac{1}{m_{sat}}F_{earth,y}.\]
1.2.1 Exercises
Calculate the acceleration in the \(x\) direction experienced by the satellite if it is in position \((1 \times 10^{7},0)\).
Calculate the acceleration in the \(y\) direction experienced by the satellite if it is in position \((1 \times 10^{7},0)\).
Calculate the acceleration experienced by the satellite if it is in position \((-4 \times 10^{6},8.5 \times 10^{6})\).
1.3 Systems of Differential Equations
We are now in a position to write a system of differential equations governing the motion of the satellite. First, remember the relationship between position, \(x\), and velocity in the \(x\) direction, \(vx\). \[x'=vx.\] Similarly, \[y'=vy.\] Now, recall from the section above that \[vx'=\frac{1}{m_{sat}}F_{earth,x} = \frac{1}{m_{sat}}\cdot \frac{G\cdot m_{sat}\cdot m_{earth} \cdot (x_{earth}-x)}{((x_{earth}-x)^2+(y_{earth}-y)^2)^{3/2}} = \frac{-G\cdot m_{earth}\cdot x}{(x^2+y^2)^{3/2}},\] and \[vy'=\frac{1}{m_{sat}}F_{earth,y}=\frac{1}{m_{sat}}\cdot \frac{G\cdot m_{sat}\cdot m_{earth} \cdot (y_{earth}-y)}{((x_{earth}-x)^2+(y_{earth}-y)^2)^{3/2}} = \frac{-G\cdot m_{earth}\cdot y}{(x^2+y^2)^{3/2}}.\] This is a system of four ordinary differential equations with which we can model the position and velocity of our satellite, assuming that only gravitational force from earth acts on the satellite.
1.3.1 Example
The system of differential equations governing the motion of our satellite is \[\begin{aligned} &x'=vx\\ &y'=vy\\ &vx'=\frac{1}{m_{sat}}F_{earth,x}\\ &vy'=\frac{1}{m_{sat}}F_{earth,y}. \end{aligned}\]
We are not going to be able to find an analytical solution for this system of equations (one does exist, but if the model becomes any more complicated, then no analytical solution exists). Instead, we’ll use Euler’s method to predict the future motion of the satellite.
This is a system of four differential equations, so we need four initial conditions. Suppose that the initial conditions are
\[x(0)=4.218 \times 10^{7} \text{ m},\quad y(0)=0 \text{ m},\quad vx(0)=0\text{ m/s},\quad vy(0)= 3.074 \times 10^{3}\text{m/s}.\]
We can use Euler’s method to approximate a solution to the system of ODEs using these initial conditions. Euler allows you to name the equations, so that it is easy to keep track of which equations are which, and which initial conditions are which. When naming equations, use _t to indicate derivative.
The units we have chosen are kilograms, meters, and seconds. To simulate the position of the satellite for one day, we will need to go from time \(t=0\) seconds to \(t=24\cdot60\cdot60 = 8.64 \times 10^{4}\) seconds. We will use a small step size, 20 seconds, to help Euler produce accurate results.
soln=Euler(c(x_t=vx,
y_t=vy,
vx_t=1/mass_sat*Fearth_x(x,y),
vy_t=1/mass_sat*Fearth_y(x,y))~x&y&vx&vy,
ic=c(x=4.218e7,y=0,vx=0,vy=3074),
tlim=c(0,24*60*60),stepSize=20)We can view the results of our simulation.
| t | x | y | vx | vy |
|---|---|---|---|---|
| 0 | 42180000 | 0.0 | 0.000000 | 3074.000 |
| 20 | 42180000 | 61480.0 | -4.480814 | 3074.000 |
| 40 | 42179910 | 122960.0 | -8.961614 | 3073.993 |
| 60 | 42179731 | 184439.9 | -13.442390 | 3073.980 |
| 80 | 42179462 | 245919.5 | -17.923132 | 3073.961 |
| 100 | 42179104 | 307398.7 | -22.403832 | 3073.935 |
We can plot the \(x\) and \(y\) coordinates of our satellite against time.
Finally, we can plot the trajectory of the satellite without a time axis.
We can see that our differential equations, combined with Euler’s method, predict that with these initial conditions, the satellite will orbit the earth. It further predicts that the satellite will make nearly one complete orbit every 24 hours. Indeed, the parameters from this satellite are approximately those of the GOES-16 satellite operated by NOAA and NASA. That satellite does in fact move with a speed in such a way as to keep it geostationary, which means they keep pace with the rotation of the earth, so one cycle every 24 hours.
Notice that it looks like our satellite won’t quite make an orbit. This is due to error in Euler’s method. We could make the step size smaller to help Euler’s method be more accurate, but no step size will result in perfect results. We will always have to accept that the results of Euler’s methods are estimates of where the satellite will be at the given time and how fast it will be moving.
If we want, we can use animateODESolution to animate our results, so that we can see our satellite orbit earth. The argument duration specifies how long it should take to get through the animation. By default this number is 5 seconds, but for this animation 10 seconds looks better.
1.3.2 Exercises
Use the same initial conditions as above, \[x(0)=4.218 \times 10^{7} \text{ m},\quad y(0)=0 \text{ m},\quad vx(0)=0\text{ m/s},\quad vy(0)= 3.074 \times 10^{3}\text{m/s},\] but simulate the position and velocity of the satellite over the next week.
Suppose that somehow the satellite got in a state where \[x(0)=4.218 \times 10^{7} \text{ m},\quad y(0)=0 \text{ m},\quad vx(0)=-2000\text{ m/s},\quad vy(0)= 3500\text{m/s}.\]
Simulate the motion of the satellite over the next eight days.
1.3.3 Example: Decomissioning the Satellite
After a long and fruitful mission, the operators of the satellite decide it is time for decommissioning. To do so, the operators fire thrusters that slow the satellite’s velocity. At the end of thruster firing, the state of the satellite is \[x=4.218 \times 10^{7} \text{ m},\quad y=0\text{ m},\quad vx=0\text{ m/s},\quad vy=1.5 \times 10^{3}\text{m/s}.\] We can simulate the position of the satellite over the next 10 hours using Euler’s method, with the same system of ODEs, and initial conditions that reflect the state above. We need a very small step size to get accurate results for this example.
soln=Euler(c(x_t=vx,
y_t=vy,
vx_t=1/mass_sat*Fearth_x(x,y),
vy_t=1/mass_sat*Fearth_y(x,y))~x&y&vx&vy,
ic=c(x=4.218e7,y=0,vx=0,vy=1500),
tlim=c(0,10*60*60),stepSize=1)If we plot the trajectory of the satellite, we find that, as expected, the satellite no longer makes a circular orbit of the earth. Instead, the trajectory is now an ellipse.
plotPoints(y~x,data=soln,col='firebrick',xlim=c(-1e7,5e7),ylim=c(-2e7,2e7),asp=get_asp(soln))
add.earth()
It looks like the satellite might hit the earth. To investigate, remember that the radius of the earth is approximately \(6.378 \times 10^{6}\) meters. At time \(t\), the distance from the satellite to the center of the earth is \(\sqrt{x^2(t)+y^2(t)}\). We can plot the distance from the satellite to the center of the earth against time, and see if the distance ever slips below \(6.378 \times 10^{6}\).
It looks like there is a possibility that the satellite collides with the earth sometime around time \(t=19000\) and/or sometime around time \(t=20000\). We can zoom in on these parts of the graph to confirm or reject our suspicions.
plotPoints(sqrt(x^2+y^2)~t,data=soln,xlim=c(16000,20000),ylim=c(0,1e7))
plotFun(0*x+6.378e6~x,add=TRUE,col="black",lty=2,lwd=2)
It looks like around time \(t=18000\), or \(5\) hours from the initial time, the satellite will impact the earth.
This graphical analysis gives us an idea when the satellite impacts the earth, but it does not tell us at what position this happens, nor with what velocity the satellite will impact the earth. Both of these are important questions that the output of Euler’s method can provide us. If we look at the data in soln we find the answers we seek.
| t | x | y | vx | vy | |
|---|---|---|---|---|---|
| 18039 | 18038 | -4822511 | 4182963 | -4133.123 | -9550.566 |
| 18040 | 18039 | -4826644 | 4173413 | -4125.735 | -9556.975 |
| 18041 | 18040 | -4830770 | 4163856 | -4118.329 | -9563.378 |
| 18042 | 18041 | -4834888 | 4154292 | -4110.906 | -9569.776 |
| 18043 | 18042 | -4838999 | 4144723 | -4103.466 | -9576.169 |
Unfortunately, the time we are interested in does not occur until row 18040 in the solution data. It isn’t reasonable to scroll through 18040 rows of data and calculate distance to the earth in each row. Advanced data manipulation is beyond the scope of this class. However, the best way to solve this problem is to use the filter command, which filters tables of data to find rows that meet criteria we specify. To find times at which the distance from the satellite to the center of the earth is less than \(6.378 \times 10^{6}\) we could use the following command.
| t | x | y | vx | vy |
|---|---|---|---|---|
| 18040 | -4830770 | 4163856 | -4118.329 | -9563.378 |
| 18041 | -4834888 | 4154292 | -4110.906 | -9569.776 |
| 18042 | -4838999 | 4144723 | -4103.466 | -9576.169 |
| 18043 | -4843103 | 4135146 | -4096.008 | -9582.557 |
| 18044 | -4847199 | 4125564 | -4088.534 | -9588.939 |
| 18045 | -4851287 | 4115975 | -4081.042 | -9595.316 |
From this table we can see that at time \(t=18040\) the satellite will have distance to center of the earth \[\sqrt{(-4830770)^2+(4163856)^2}=6.378 \times 10^{6} \text{ meters}.\] This is less than the radius of the earth, and so at this time the satellite impacts the earth. It impacts the earth at position \(x=-4830770\) m, \(y=4163856\). The speed with which the satellite impacts the earth is \[\sqrt{(-4118.329)^2+(-9563.378)^2} =1.041 \times 10^{4} \text{ m/s}.\]
1.4 The Moon
1.4.1 Gravitational Force
We have considered the gravitational impact of the earth on the satellite, but what about the moon? The center of the moon is about \(3.844 \times 10^{8}\) meters from the center of the earth. Assume that the center of the moon is at \((3.844 \times 10^{8},0)\). The mass of the moon is about \(7.348 \times 10^{22}\) kg.
We can use the general equations for gravitational force to calculate the gravitational force acting on the satellite generated by the moon.
\[F_{moon,x} = \frac{G\cdot mass_{moon}\cdot mass_{sat} (x_{moon}-x)}{\left((x_{moon}-x)^2+(y_{moon}-y)^2\right)^{3/2}}\]
\[F_{moon,y} = \frac{G\cdot mass_{moon}\cdot mass_{sat} (y_{moon}-y)}{\left((x_{moon}-x)^2+(y_{moon}-y)^2\right)^{3/2}}.\]
We can create R functions that compute the gravitational force exerted by the moon on the satellite when the satellite is in position \((x,y)\) by filling in the above formula.
mass_m=7.347673e+22 #Mass of the moon, in kg.
Fmoon_x=makeFun(G*mass_m * mass_sat*(3.844e8-x)/((3.844e8-x)^2+(0-y)^2)^(3/2)~x&y)1.4.1.1 Exercise
Write a function that computes the gravitational force in the \(y\) direction of the moon upon the satellite when the satellite is at position \((x,y)\).
Calculate the gravitational force in the \(x\) and \(y\) directions exerted by the moon on the satellite when the satellite is in position \(x=4.218 \times 10^{7}\) m, \(y=0\) m.
Calculate the gravitational force in the \(x\) and \(y\) directions exerted by the moon on the satellite when the satellite is in position \(x=4.218 \times 10^{7}\) m, \(y=0\) m. Compare it to the gravitational force exerted by the earth.
1.4.2 Systems of Differential Equations
Incorporating the effect of the moon on our satellite is not difficult. Remember that Newton’s second law says that \[ma=F,\] where \(F\) is the net force acting on the object. If we are now incorporating gravity from the moon in our model, then our force has two parts: gravity from the earth and gravity from the moon. Therefore, \[\begin{aligned} &m_{sat}vx'=F_{earth,x} + F_{moon,x}\\ &m_{sat}vy'=F_{earth,y} + F_{moon,y}. \end{aligned}\] We can rewrite those equations to get \(vx'\) and \(vy'\) alone, and we end up with the system of ODEs \[\begin{aligned} &x'=vx\\ &y'=vy\\ &vx'=\frac{1}{mass_{sat}}\left(F_{earth,x} + F_{moon,x}\right)\\ &vy'=\frac{1}{mass_{sat}}\left(F_{earth,y} + F_{moon,y}\right). \end{aligned}\]
1.4.2.1 Exercise
Predict the motion of the satellite over the next four hours if it is initially in position \[x(0)=3.844 \times 10^{8}, \quad y(0)=1.739 \times 10^{6},\quad vx(0)=-1.679 \times 10^{3},\quad vy(0)=0.\] Include the effect of the earth. Use a small step size, 2 seconds.
Predict the motion of the satellite over the next five days if it is initially in position \[x(0)=3.844 \times 10^{8}, \quad y(0)=4.428 \times 10^{6},\quad vx(0)=-1.263 \times 10^{3},\quad vy(0)=631.\] Include the effect of the earth. Use step size 10 seconds.
1.5 Lagrange Points
Including gravitational effects from the Earth and the moon, the differential equations guiding the motion of the satellite are \[\begin{aligned} &x'=vx\\ &y'=vy\\ &vx'=\frac{1}{mass_{sat}}\left(F_{earth,x} + F_{moon,x}\right)\\ &vy'=\frac{1}{mass_{sat}}\left(F_{earth,y} + F_{moon,y}\right). \end{aligned}\]
Let us calculate the equilibrium solutions for this system of ODEs. To calculate equilibrium solutions, we will need to solve
\[\begin{aligned}
&0=vx\\
&0=vy\\
&0=\frac{1}{mass_{sat}}\left(F_{earth,x} + F_{moon,x}\right)\\
&0=\frac{1}{mass_{sat}}\left(F_{earth,y} + F_{moon,y}\right).
\end{aligned}\]
findZeros can do this for us. We will provide the option roundDigits=-3 so that we round all potential solutions to the thousands (not thousandths) place. If we don’t include this option, findZeros will find many solutions that are all similar to, but different from, each other. We also give findZeros the option within=3.844e8. By default, findZeros will look for solutions that have coordinates in the range -100 to 100. The within option tells findZeros to look for solutions with coordinates in the range \(-3.844 \times 10^{8}\) to \(3.844 \times 10^{8}\).
findZeros(c(vx,
vy,
1/mass_sat*(Fearth_x(x,y)+Fmoon_x(x,y)),
1/mass_sat*(Fearth_y(x,y)+Fmoon_y(x,y)))~x&y&vx&vy,
within=3.844e8,roundDigits=-3)| x | y | vx | vy |
|---|---|---|---|
| 346020000 | 0 | 0 | 0 |
We can see that there is one equilibrium solution. This result tells us that if we can get the satellite in to position \(x=346020000\), \(y=0\) with velocity \(vx=0\), \(vy=0\), then the satellite will never move. This is interesting, if we wanted to park a satellite in a fixed position and let the earth rotate below the satellite, an equilibrium position would be a good place to park the satellite, as it would take no energy to keep the satellite in this position. Positions in which a satellite can be “parked” and never move, while expending no fuel, are called Lagrange Points. Our model of the earth and the moon does not account for the motion of earth and the moon, and so it does not predict the existence of other Lagrange points (there are four more of them). The other Lagrange points are actually used for parking satellites, while the one our model found is not used due to its long distance from earth.
It turns out that the equilibrium solution we found is an unstable equilibrium solution. That means that the satellite will only stay in position if it is perfectly in position. That is not possible, so our satellite won’t really stay in position with no energy spent. Instead, it will have to occasionally adjust position to guide itself back to the Lagrange point. It turns out that two of the other Lagrange points are stable.
1.6 Incorporating Wind Resistance
In our analysis of a satellite returning to earth in 1.3.3, we neglected any resistance the satellite would encounter from the atmosphere (drag). This is an important factor to consider, as it will slow the satellite considerably. For our purposes in this section, we will assume that the satellite is not being recovered for future use or analysis, and therefore there is no parachute trying to land the satellite gently.
The drag force is generally depends on a) the shape, size, and material of the object under consideration, b) the speed with which the object is travelling, and c) the density of the air.
Again, we will consider the drag force in the \(x\) and \(y\) directions. One common model is
\[F_{drag,x} = -\frac12 \alpha \rho |vx| vx \] \[F_{drag,y} = -\frac12 \alpha \rho |vy| vy \] Drag force, in this example, is measured in Newtons.
The coefficient \(\alpha\) is the drag coefficient. It takes into account the size, shape, and material of the object. For a spherical satellite roughly the size of GOES-16, a good estimate is \(\alpha \approx 43.2\).
The coefficient \(\rho\) is the density of the atmosphere. This changes with height above the surface of the earth. One simple model for the density of the atmosphere \(h\) meters above the surface of the earth is \[\rho = 1.2e^{-h/8400}.\] The height of the satellite above the earth depends on the satellite’s position, and can be calculated as \[h=\sqrt{x^2+y^2}-6.378 \times 10^{6}.\] If we put this all together, we find that the drag forces acting on the satellite when it is in position \((x,y)\) and moving with velocities \(vx\), \(vy\) are
\[\begin{aligned} &F_{drag,x}=-\frac{1}{2}\cdot 43.2\cdot 1.2\cdot e^{-\frac{\left(x^2+y^2\right)^{1/2}-6.378 \times 10^{6}}{8400}}|vx|vx = -25.2e^{-\frac{\left(x^2+y^2\right)^{1/2}-6.378 \times 10^{6}}{8400}}|vx|vx\\ &F_{drag,y}=-\frac{1}{2}\cdot 43.2\cdot 1.2\cdot e^{-\frac{\left(x^2+y^2\right)^{1/2}-6.378 \times 10^{6}}{8400}}|vy|vy = -25.2e^{-\frac{\left(x^2+y^2\right)^{1/2}-6.378 \times 10^{6}}{8400}}|vy|vy.\\ \end{aligned}\]
Though these formulas are complicated, it is easy enough to convert them into R code. We can create a function to calculate the drag force in the \(x\) direction if the satellite is in position \(x\), \(y\) and moving with velocity \(vx\), \(vy\) as below.
1.6.1 Exercise
Write a function that calculates the drag force in the \(y\) direction acting on the satellite if it is in position \((x,y)\) and has velocities \(vx\), \(vy\).
Calculate the drag forces encountered by the satellite if it is in position \(x=-5.037 \times 10^{6}\), \(y=4.131 \times 10^{6}\) and moving with velocities \(vx=-4.097 \times 10^{3}\), \(vy=-9.504 \times 10^{3}\).
1.7 Incorporating Drag into ODEs
Adding drag into the model is not difficult. Once again, since \[m\cdot a=F\] where \(F\) is the net force acting on the satellite, we just need to add in one more term for the drag force. We end up with \[\begin{aligned} &x'=vx\\ &y'=vy\\ &vx'=\frac{1}{mass_{sat}}\left(F_{earth,x} + F_{moon,x}+F_{drag,x}\right)\\ &vy'=\frac{1}{mass_{sat}}\left(F_{earth,y} + F_{moon,y}+F_{drag,y}\right). \end{aligned}\]
1.7.1 Example: Decomissioning the Satellite with Drag
Let us reconsider the scenario in 1.3.3. The operators of the satellite have decided to decommission, and have fired thrusters to slow down the satellite so that it can return to earth. The first time the satellite will encounter any measurable drag from the Earth’s atmosphere is when it is position \[x=-4.775 \times 10^{6},\quad y=4.687 \times 10^{6}\] and is moving with velocity \[vx=-4.496 \times 10^{3},\quad vy=-9.134 \times 10^{3}.\]
We can use Euler’s method to compare the trajectory of the satellite over the next 10 minutes of motion where we only account for Earth’s gravitational force and when we account for the moon (negligible here) and drag (not negligible here). We will use very small step sizes to help ensure accuracy.
First, Euler’s method accounting only for Earth’s gravity. This will reproduce the results of 1.3.3.
soln_nodrag=Euler(c(x_t=vx,
y_t=vy,
vx_t=1/mass_sat*Fearth_x(x,y),
vy_t=1/mass_sat*Fearth_y(x,y))~x&y&vx&vy,
tlim=c(0,10*60),
ic=c(x=-4.775e6,y=4.687e6,vx=-4.496e3,vy=-9.134e3),
stepSize = 0.1)Next, we account for the moon’s gravity and drag.
soln_drag=Euler(c(x_t=vx,
y_t=vy,
vx_t=1/mass_sat*(Fearth_x(x,y)+Fmoon_x(x,y)+Fdrag_x(x,y,vx,vy)),
vy_t=1/mass_sat*(Fearth_y(x,y)+Fmoon_y(x,y)+Fdrag_y(x,y,vx,vy)))~x&y&vx&vy,
tlim=c(0,10*60),
ic=c(x=-4.775e6,y=4.687e6,vx=-4.496e3,vy=-9.134e3),
stepSize = 0.1)Now, we can visually compare the trajectory of the satellite under the two different models.
plotPoints(y~x,data=soln_nodrag,col="cadetblue4",asp=get_asp(soln_nodrag))
plotPoints(y~x,data=soln_drag,col="firebrick4",add=TRUE)
add.earth()
This figure is a little hard to read, but we can see that the two trajectories are certainly quite different. We can zoom in on the with-drag trajectory to examine further. The window to look at isn’t at all obvious, and is carefully selected so that a) the figure looks good, and b) the figure is not misleading.
plotPoints(y~x,data=soln_drag,col="cadetblue4",xlim=c(-5.5e6,-4.5e6),ylim=c(3.5e6,4.5e6),asp=1)
plotPoints(y~x,data=soln_nodrag,col="firebrick4",add=TRUE)
add.earth()
Now we can see that the two trajectories are very different. The trajectory without drag continues on a familiar elliptical path as it nears the earth. The trajectory with drag, however, starts to “skip” off the atmosphere as it nears earth. It eventually slows and then returns to earth nearly straight down (from earth’s perspective).
To compare further, let us compare impact information under the two models. For the no-drag model, we can find impact information by finding the time at which the satellite is less than one earth-radius from the center of the earth.
| t | x | y | vx | vy |
|---|---|---|---|---|
| 108.2 | -5221086 | 3663150 | -3718.726 | -9780.792 |
| 108.3 | -5221458 | 3662172 | -3717.924 | -9781.355 |
| 108.4 | -5221829 | 3661194 | -3717.121 | -9781.918 |
| 108.5 | -5222201 | 3660216 | -3716.319 | -9782.480 |
| 108.6 | -5222573 | 3659238 | -3715.516 | -9783.043 |
| 108.7 | -5222944 | 3658259 | -3714.713 | -9783.605 |
We can see that under this model the satellite impacts the earth at time \(t=108.2\) seconds. It impacts the earth at position \(x=-5221086\), \(y=3663150\), with speed \(\sqrt{(-3718.726)^2+(-9780.792)^2}=1.046 \times 10^{4}\) m/s.
To find information about impact under the with-drag model, we can filter the data similarly.
| t | x | y | vx | vy |
|---|---|---|---|---|
| 442.7 | -5209438 | 3679756 | 44.06178 | -37.10016 |
| 442.8 | -5209434 | 3679752 | 44.04648 | -37.08722 |
| 442.9 | -5209430 | 3679748 | 44.03118 | -37.07429 |
| 443.0 | -5209425 | 3679745 | 44.01590 | -37.06138 |
| 443.1 | -5209421 | 3679741 | 44.00063 | -37.04847 |
| 443.2 | -5209416 | 3679737 | 43.98537 | -37.03556 |
We can see that under this model, the satellite impacts the earth much later, at time \(t=442.7\) seconds. It impacts the earth at position \(x=-5209439\), \(y=3679756\), with speed \(\sqrt{(44.06181)^2+(-37.10014)^2}=57.6\) m/s. Quite different from the no-drag model indeed!