Multibody Gravitational Interaction Problem
Student Kristian Dolghier
Professor Pavan Pillalamarri
Barrett Honors College Honors Contract
May 8th 2020
One earth, Escape velocity
Setting, M=5.972*10^24 kilograms and R=3.67*10^6 meters(radius of the earth), with a velocity in the positive y
At a speed of 13700 m/s, the object begins to escape the earth's orbit and no longer is in a stable orbit
At the higher speed of 14700 m/s the object is truly in escape velocity
Note: In real life the escape velocity is lower, but it is very difficult to see on the graph so I included the first escape velocity where it is easily visually distinguishable that it is escaping.
Geosynchronous orbit
To find the radius required we must first do this calculation
GM/r=v^2/r where v=sqrt(G M / r)=2 pi r/T from these two equations we find the radius is equal to r=(GMT^2/4pi^2)^⅓=((6.67*10^-11*5.972*10^24*(86400)^2/(4pi^2)))^⅓= 4.22*10^7 meters
Then replace the current radius with this.so we find the velocity is v=2 pi r/T= 2pi4.22*10^7/86400= 3068.9m/s
Inserting this speed and radius we get this geostationary orbit
2 planets interacting with 1 object
Case 1:
Earth and moon and the ISS
The ISS orbits the earth at a minimum height of 205 mi (329915.52m) and a maximum of 255 mi (410382.72m). Getting an average height we find that the mean is 370149m. Therefore with the radius of the earth the total R is 4040149m. The speed of the ISS is 7.66 km/sec or 7660m/s. Finally the mass of the earth is 5.972*10^24kg.
The moon has a mass of 7.348 × 10²² kg, and has an average distance from the earth of 384,403 km or 3.48*10^8m. The for loop is run for 1*10^7 iterations( it always is ) and the following graph is the result.
Case 2:
Two earth-like planets with different distances from the object.
With two planets with the same mass of 5.972*10^24kg, and the same initial velocity of 7660m/s, the only difference is the radius. The planet to the left of the object is 3.67*10^7meters away while the object on the right is 3.67*10^6 meters away.
Case 3:
A moon-like object and an earth-like object with the moon being closer to the object.
Set one of the masses to be the mass of the earth(5.972*10^24 kg) and the radius of the moon, where the radius is the radius of the earth or 3.67*10^6 m. The second planet has a mass of 5.972*10^22 kg and a radius of 3.67*10^5 m. The velocity is still 7660 m/s.
The following is the graph.
Case 4
The effects of the sun on the ISS orbiting earth.
Setting the mass of the sun to be M=1.989*10^30 kg, the R is found by the distance from the sun subtracting the height of the ISS and adding the radius R=1.5*10^11-329915.52+6.95*10^8 m = 1.5069*10^11 m, the velocity is still 7660 m/s. The Mass of the earth is the same at 5.972*10^24 kg with a radius of 3.67*10^6+329915.52. The interaction is as follows.
Note: this is done for 600000/.1 iterations
Reflection
This was a very hard project to finish and basically complete in around a week. It took a little bit to create the code for the gravitational interaction of a single object, but that code seemed to be wasted. Instead of thinking of new ways to make my code work for a two body system I instead was trying to retrofit the old code to the new one as I have done similarly before. Instead of looking at the core concepts of my code and trying to find where the problem is with the idea of my code, I was tracking back where I was getting errors from and then became confused as the code was the exact same for object one and yet I was getting different values. In the future, I must look at the flowchart behind my code, and not just try to put patches on bigger problems but look at the core of my code and figure out exactly where my conceptualization of the code went wrong.
Matlab Code
%Gravitational interaction
clc, clear
close all
M=5.972*10^24; %kilograms;
R=3.67*10^6; %meters;
v=7660; %meters per second;
vx(1)=.707*v;
vy(1)=.707*v;
x(1)=0;
Dist(1)=R;
y(1)=0;
M2=5.972*10^23; %m0on in kg
R2=-3.67*10^6; %meters; %cant set radius to 0, must be very small or error will occur
G=6.67*10^-11;
delta_t=.01;
run=10000/delta_t; %10 is the number of seconds it runs for
for i=1:run
Dist(i)=sqrt((x(i)-R).^2+y(i).^2);
a(i)=-G*M/(Dist(i).^2);
ax(i)=a(i)*(x(i)-R)/Dist(i);
ay(i)=a(i)*y(i)/Dist(i);
Dist2(i)=sqrt((x(i)-R2).^2+y(i).^2);
a2(i)=-G*M2/(Dist2(i).^2);
ax2(i)=a2(i)*(x(i)-R2)/Dist2(i);
ay2(i)=a2(i)*y(i)/Dist2(i);
axn(i)=ax(i)+ax2(i);
ayn(i)=ay(i)+ay2(i);
vx(i+1)=vx(i)+axn(i)*delta_t;
vy(i+1)=vy(i)+ayn(i)*delta_t;
x(i+1)=x(i)+vx(i)*delta_t+.5*axn(i)*delta_t.^2;
y(i+1)=y(i)+vy(i)*delta_t+.5*ayn(i)*delta_t.^2;
%plot(x(1:i),y(1:i),'color','r')
%drawnow()
end
plot(x,y) %must be no (i)