Author Topic: Simulate orbit using Matlab  (Read 43611 times)

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Simulate orbit using Matlab
« on: 11/22/2008 06:37 PM »
Hi!! I'm trying to make an orbit simulation with Matlab, and I'm having some trouble making it simulate real scenarios such as the ISS, the Moon or a sattelite in Geosynchronous orbit.

My model is extremely simple, I avoided so many facts I'm kind of embarassed presenting this here, but perhaps you could help me with some advices.

Basically what I did was just consider the gravitational attraction between two spots in a cartesian coordinate system. Every body is just a spot in the system. The Earth does not move, it's not affected by the satellite's gravitational pull.

The key is the orbital speed, I believe my mistake is here. The sattelite is given an initial distance from the center of the earth, an initial orbital speed and a mass.
Position is calculated as:
   pos_sat_x = pos_sat_x + vel_sat_x*escala_tiempo;
   pos_sat_y = pos_sat_y + vel_sat_y*escala_tiempo;
Orbital speed as:
   vel_sat_x  = vel_sat_x + acc_tierra_x*escala_tiempo;
   vel_sat_y  = vel_sat_y + acc_tierra_y*escala_tiempo;
And acceleration as:
   distancia = sqrt(pos_sat_y^2+pos_sat_x^2);
   gravedad  = -((G*(masa_tierra*masa_satelite))/(distancia^2));
   angulo = atan(pos_sat_y/pos_sat_x);
   acc_tierra_y = gravedad*abs(sin(angulo));
   acc_tierra_x = gravedad*abs(cos(angulo));
(escala_tiempo is a variable used to make the simulation faster)

I'll obviously never be able to simulate it perfectly, but there's something I'm missing and I'm apologizing for any lack of respect to y'all by trying to simplify something so complex in some few potentially wrongly-designed equations.

Here you got a trajectory calculated with a mass of 14000 kg, an initial orbital speed of 300000 m/s (yes, 300 km/s!! :-\), and a starting distance of 35786 km. I tried to calculate this geosynchronous orbit with the regular orbital speed of 3 km/s but it just 'fell to earth'.

Hope you enjoy it, greetings from Argentina!!!


Online Jorge

  • Senior Member
  • *****
  • Posts: 6180
  • Liked: 1
Re: Simulate orbit using Matlab
« Reply #1 on: 11/22/2008 07:34 PM »
Hi!! I'm trying to make an orbit simulation with Matlab, and I'm having some trouble making it simulate real scenarios such as the ISS, the Moon or a sattelite in Geosynchronous orbit.

My model is extremely simple, I avoided so many facts I'm kind of embarassed presenting this here, but perhaps you could help me with some advices.

Basically what I did was just consider the gravitational attraction between two spots in a cartesian coordinate system. Every body is just a spot in the system. The Earth does not move, it's not affected by the satellite's gravitational pull.

The key is the orbital speed, I believe my mistake is here. The sattelite is given an initial distance from the center of the earth, an initial orbital speed and a mass.
Position is calculated as:
   pos_sat_x = pos_sat_x + vel_sat_x*escala_tiempo;
   pos_sat_y = pos_sat_y + vel_sat_y*escala_tiempo;
Orbital speed as:
   vel_sat_x  = vel_sat_x + acc_tierra_x*escala_tiempo;
   vel_sat_y  = vel_sat_y + acc_tierra_y*escala_tiempo;
And acceleration as:
   distancia = sqrt(pos_sat_y^2+pos_sat_x^2);
   gravedad  = -((G*(masa_tierra*masa_satelite))/(distancia^2));
   angulo = atan(pos_sat_y/pos_sat_x);
   acc_tierra_y = gravedad*abs(sin(angulo));
   acc_tierra_x = gravedad*abs(cos(angulo));
(escala_tiempo is a variable used to make the simulation faster)

Here may be your problem: gravedad is the magnitude of the gravitational acceleration but you're using the equation for gravitational force. The two are related by Newton's law (F=ma, or a=F/m, where m=masa_satelite). So you should divide the right hand side by masa_satelite, which makes it disappear from the equation:

   gravedad  = -G*masa_tierra/(distancia^2);

This is valid whenever masa_satelite << masa_tierra.

Quote
Here you got a trajectory calculated with a mass of 14000 kg, an initial orbital speed of 300000 m/s (yes, 300 km/s!! :-\), and a starting distance of 35786 km. I tried to calculate this geosynchronous orbit with the regular orbital speed of 3 km/s but it just 'fell to earth'.

With the extraneous masa_satelite factor in the equation, no wonder - you're making Earth's gravity 14,000 times stronger than it really is.

You can simplify the gravitational equation by substituting the gravitational parameter (mu_tierra=G*masa_tierra=398600 km^3/s^2).

You can also simplify the trig functions out of the equations by utilizing standard trig identities: sin(angulo)=pos_sat_y/distancia and cos(angulo)=pos_sat_x/distancia. That leaves you with:

   gravedad  = -mu_tierra/(distancia^3);
   acc_tierra_y = gravedad*pos_sat_y;
   acc_tierra_x = gravedad*pos_sat_x;

(note that distancia squared becomes distancia cubed.)
JRF

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Re: Simulate orbit using Matlab
« Reply #2 on: 11/22/2008 08:26 PM »
thank you so much!!! It was as you said, now all the parameters work!!

you said "This is valid whenever masa_satelite << masa_tierra.". What about the Moon?? it works in my model...

Online Jorge

  • Senior Member
  • *****
  • Posts: 6180
  • Liked: 1
Re: Simulate orbit using Matlab
« Reply #3 on: 11/22/2008 09:31 PM »
thank you so much!!! It was as you said, now all the parameters work!!

you said "This is valid whenever masa_satelite << masa_tierra.". What about the Moon?? it works in my model...

If masa_satelite !<< masa_tierra, then use:

G*(masa_tierra+masa_satelite)

or if you prefer, multiply through by G and substitute the gravitational parameters to get:

(mu_tierra+mu_satelite)

so the whole line would look like:

gravedad  = -(mu_tierra+mu_satelite)/(distancia^3);
Since mu_luna = 4903 km^3/s^2, the error of just using mu_tierra is only a little more than 1%. It will appear to work fine for short term propagation, though.

If you want to make the model even more accurate I would recommend using a higher-order numerical integration method. What you have now is essentially Euler's first-order method. But in order to make your code more amenable to swapping out integrators, you probably want to first reformulate pos_sat, vel_sat, and acc_tierra as vectors. That plays into the strengths of Matlab anyway, and also makes it trivial to extend your code to three dimensions.
JRF

Online Herb Schaltegger

  • I used to be a rocket scientist
  • Full Member
  • ****
  • Posts: 1054
  • Liked: 116
  • Murfreesboro TN
Re: Simulate orbit using Matlab
« Reply #4 on: 11/23/2008 01:37 AM »
Jorge -

I can tell why they let you teach this stuff to the astronauts. :)
Quote
Less talk, more rocket.
~ Robotbeat

Offline RedSky

  • Veteran
  • Full Member
  • ****
  • Posts: 449
  • Liked: 2
Re: Simulate orbit using Matlab
« Reply #5 on: 11/23/2008 02:27 AM »
Have you seen this site?  Once you get the hang of launching new planets (mouse down, pull back and release... like pinball), its really an addictive time sink.  Its actually quite hard to keep a stable solar system going.  You can change the sizes (masses) of the planets, and what's really neat is getting a small moon to orbit a planet.

http://dan-ball.jp/en/javagame/planet/

Offline m330533

  • STS Fanatic
  • Member
  • Posts: 67
  • Liked: 0
  • Duesseldorf, Germany
    • Air France Cargo - KLM Cargo
Re: Simulate orbit using Matlab
« Reply #6 on: 11/23/2008 09:49 AM »
Jorge -

I can tell why they let you teach this stuff to the astronauts. :)
Yeah, exactly THAT's what came up into my mind as well !

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Re: Simulate orbit using Matlab
« Reply #7 on: 11/23/2008 02:11 PM »
Thank you Jorge!!!

Now I'm moving to simulate an Earth-Moon trajectory. I know Wiki is not the best source but I couldn't find anywhere else. It states that post TLI burn the S-IVB+CSM+LM had a speed of 11.2 km/s is that correct?? When I load it to my simulation, the spacecraft describes an incredibly eccentric orbit, you can almost confuse it with a straight line.......I believe this should be OK since 11.2km/s is quite a speed.

Now I got a few questions:
Is it OK to consider the initial velocity vector orthogonal to the position vector?? or the TLI burn was carried out with the engines facing a different angle??
Was this the actual escape speed for Apollo missions??
Thid they made a burn to enter moon orbit??

I'm really embarassed about my ignorance :P thank you so much for your help!!!

Online Jorge

  • Senior Member
  • *****
  • Posts: 6180
  • Liked: 1
Re: Simulate orbit using Matlab
« Reply #8 on: 11/24/2008 01:20 AM »
Thank you Jorge!!!

Now I'm moving to simulate an Earth-Moon trajectory. I know Wiki is not the best source but I couldn't find anywhere else. It states that post TLI burn the S-IVB+CSM+LM had a speed of 11.2 km/s is that correct??

For Apollo trajectory data it's best to go straight to the source. The original Apollo mission reports can be downloaded from the NASA History Office's Apollo Lunar Surface Journal:

http://history.nasa.gov/alsj/alsj-mrs.html

Each report has a table of trajectory parameters following each burn. You'll need to convert them to metric. For example, Apollo 17 had a post-TLI velocity of 35589.6 ft/s, or 35589.6*0.3048 = 10847.7 m/s.

Quote
When I load it to my simulation, the spacecraft describes an incredibly eccentric orbit, you can almost confuse it with a straight line.......I believe this should be OK since 11.2km/s is quite a speed.

That's correct. The post-TLI orbits were highly eccentric ellipses, with apogee considerably beyond the moon's distance.

Quote
Now I got a few questions:
Is it OK to consider the initial velocity vector orthogonal to the position vector?? or the TLI burn was carried out with the engines facing a different angle??

At the end of the TLI burn, the velocity vector was not quite orthogonal to the position vector. In the mission reports, this is given by the "flight path angle". For example, for Apollo 17 the flight path angle was 6.947 degrees at TLI cutoff. That means the velocity vector was 6.947 degrees above the local horizontal plane, which in turn is orthogonal to the velocity vector.

If you're going to use the actual post-TLI velocities and flight path angles for your model, you should also use the altitudes listed in the mission reports. For example, Apollo 17 was in a 90 nmi circular orbit at TLI ignition, but by the time of TLI cutoff it had risen to 162.4 nmi. (Also note that these are nautical miles, not statute miles. 1 nmi = 1.852 km.)

There are other trajectory parameters in each table, such as the longitude, latitude, and heading angle. You don't need those now, but save them for when you extend your simulation to three dimensions.

Quote
Thid they made a burn to enter moon orbit??

Yes. The Lunar Orbit Insertion (LOI) burn was around 0.9 km/s, and placed the CSM/LM stack in a 60x170 nmi lunar orbit (except for Apollo 17, which lowered the pericynthion to 51 nmi).

Quote
I'm really embarassed about my ignorance :P thank you so much for your help!!!

No problem. Ask any time!
JRF

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Re: Simulate orbit using Matlab
« Reply #9 on: 11/24/2008 07:49 PM »
Thank you so much Jorge, I'm now checking for the Apollo 17 Mission report.

According to what you've said, for the post-TLI x and y velocity components, I should do vel_x = 10847.7*cos(6.947) and vel_y = 10847.7*sin(6.947) is that right?? I put it on the simulation and the result wasn't that convincing, I'm attaching it

I would really need you to help me understand the parameters of the LOI, why is it:
LOI ignition:  8110.2 ft/s  (-9.90º)
LOI cutoff:    5512.1 ft/s  (0.43º)

I'm still on cartesian system, and didn't convert the positions and velocity into vectors as you see =S I feel much more conmfortable this way because it's more intuitive

thanks again!!!!

Online Jorge

  • Senior Member
  • *****
  • Posts: 6180
  • Liked: 1
Re: Simulate orbit using Matlab
« Reply #10 on: 11/24/2008 08:22 PM »
Thank you so much Jorge, I'm now checking for the Apollo 17 Mission report.

According to what you've said, for the post-TLI x and y velocity components, I should do vel_x = 10847.7*cos(6.947) and vel_y = 10847.7*sin(6.947) is that right?? I put it on the simulation and the result wasn't that convincing, I'm attaching it

You're treating the flight path angle as if it were the angle between the position vector and the velocity vector, which isn't the case. Rather, it's the angle between the local horizontal plane (which is normal to the position vector) and the velocity vector.

So, assuming you're still formulating the initial position vector as:

pos_sat_x   = radio_tierra+distancia;      % m
pos_sat_y   = 0;               % m

then the normal to the position vector would be along the y-axis, so the velocity vector should be formulated as:

vel_sat_x    = velocidad*sin(gamma);            % m/s
vel_sat_y   = velocidad*cos(gamma);            % m/s

in other words, you should reverse the sin and cos terms.

Quote
I would really need you to help me understand the parameters of the LOI, why is it:
LOI ignition:  8110.2 ft/s  (-9.90º)
LOI cutoff:    5512.1 ft/s  (0.43º)

Keep in mind the LOI parameters are expressed in a moon-centered frame, as opposed to the TLI parameters which are in an Earth-centered frame.

Seen from the moon's perspective, the spacecraft is approaching along a hyperbolic trajectory. The spacecraft performs LOI such that at cutoff, it is in a 51x170 nmi orbit and is near the point of pericynthion. For an elliptical orbit, the periapsis and apoapsis are the only two points at which the flight path angle is zero, so that explains why the flight path angle is near zero at LOI cutoff. Prior to that point, the spacecraft is still dropping toward the moon so the flight path angle is negative.

If, on the other hand, the spacecraft did not perform a burn at all, it would whip around the moon and depart along the same hyperbola. So the spacecraft needs to slow down in order to enter lunar orbit. Hence the lower velocity at LOI cutoff.

Quote
I'm still on cartesian system, and didn't convert the positions and velocity into vectors as you see =S I feel much more conmfortable this way because it's more intuitive

Vectors can be expressed as cartesian (as opposed to magnitude/direction) as well, and in fact that's the way Matlab prefers to operate, so you don't have to give up the intuitiveness. The advantage of using vectors (besides more compact code) is that you can have Matlab do the heavy lifting of propagating the state (position and velocity) for you. But go ahead and stick with scalar variables if you learn better that way.
JRF

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Re: Simulate orbit using Matlab
« Reply #11 on: 11/25/2008 01:15 AM »
Thanks again Jorge....been working and got some news =)

After 1.2 million iterations I finally calculated the trajectory after TLI, i'm attaching it, I'd need to know if it's close to the actual trajectory, some posts earlier you said the apogee was considerably away from the Moon.....and my trajectory's apogee is almost in the lunar orbit....what I don't find weird because I see no point in consuming more fuel to get an orbit further away from my target....

To try to attempt a LOI burn simulation, I'll need to put very precisly the Moon's starting point so they can meet in some spot....it's doable you think???

To finish, I'll need some further explanation about the angles, do you know any source I could read??

Thank you!!!!!!

Sergio

Online Jorge

  • Senior Member
  • *****
  • Posts: 6180
  • Liked: 1
Re: Simulate orbit using Matlab
« Reply #12 on: 11/25/2008 03:41 AM »
Thanks again Jorge....been working and got some news =)

After 1.2 million iterations I finally calculated the trajectory after TLI, i'm attaching it, I'd need to know if it's close to the actual trajectory, some posts earlier you said the apogee was considerably away from the Moon.....and my trajectory's apogee is almost in the lunar orbit....what I don't find weird because I see no point in consuming more fuel to get an orbit further away from my target....

That doesn't look quite right... the ellipse should be nearly horizontal. Make sure you convert the flight path angle to radians (multiply by pi/180) before passing it to the sin() and cos() functions.

Quote
To try to attempt a LOI burn simulation, I'll need to put very precisly the Moon's starting point so they can meet in some spot....it's doable you think???

Very tricky. Remember that the moon's true orbit is elliptical, so it won't do to simply  start propagating the moon's orbit from perigee. You need to know where the moon was at the time of Apollo 17, and for that you need an ephemeris. But for now you can play with it using trial and error.

Quote
To finish, I'll need some further explanation about the angles, do you know any source I could read??

Most any good astrodynamics text will do, but I'll look for some references.
JRF

Offline toddbronco2

  • Member
  • Full Member
  • **
  • Posts: 268
  • Liked: 0
  • Pasadena, CA
Re: Simulate orbit using Matlab
« Reply #13 on: 11/29/2008 08:04 PM »
Thanks again Jorge....been working and got some news =)

After 1.2 million iterations I finally calculated the trajectory after TLI, i'm attaching it, I'd need to know if it's close to the actual trajectory, some posts earlier you said the apogee was considerably away from the Moon.....and my trajectory's apogee is almost in the lunar orbit....what I don't find weird because I see no point in consuming more fuel to get an orbit further away from my target....

To try to attempt a LOI burn simulation, I'll need to put very precisly the Moon's starting point so they can meet in some spot....it's doable you think???

To finish, I'll need some further explanation about the angles, do you know any source I could read??

Thank you!!!!!!

Sergio

Your simulations look great!  It might be time for you to graduate to a higher fidelity model.  Have you considered simulating the Circular Restricted Three-Body Problem (CR3BP)??

You may want to consider numerically integrating the relative equations of motion that describe the acceleration of two point sources moving in circular motion around one another.  In MATLAB this is easy to do with ode45.  You can use ode45 to call a function that including the equations of motion describing the acceleration felt by the spacecraft due to the influence of the Earth and the Moon.  For example, the way I simulate this is by numerically integrating this 1st order state space differential equations:


function x_dot = CR3BP2(t,x);

n = 1;
m1 = 398600.4418;
m2 = 4902.7949;

mStar = m1+m2;
mu = m2/mStar;

x1 = x(1);
y1 = x(2);
z1 = x(3);
x_dot(1) = x(4);
x_dot(2) = x(5);
x_dot(3) = x(6);

xdot = x_dot(1);
ydot = x_dot(2);
zdot = x_dot(3);

d = (((x1+mu)^2)+y1^2+z1^2)^.5;
r = ((x1-1+mu)^2+y1^2+z1^2)^.5;

x_dot(4) = (2*n*ydot)+((n^2)*x1)-((1-mu)*(x1+mu)/(d^3))-((mu*(x1-1+mu))/(r^3));
x_dot(5) = (-2*n*xdot)+((n^2)*y1)-(((1-mu)*y1)/(d^3))-((mu*y1)/(r^3));
x_dot(6) = -(((1-mu)*z1)/(d^3))-((mu*z1)/(r^3));

x_dot = x_dot';
return


I'll attach some plots to show the kinds of cool trajectories you can get when you include the Earth and Moon's gravity fields.  The trajectory I'm attaching is a repeating free-return trajectory that encounters the moon once every 5 spacecraft orbits around the Earth.  This trajectory is NOT like the one Apollo used (it doesn't intersect the Earth's atmosphere on the return and it passes 10,000 km from the moon, but still it's pretty cool in my opinion.  The Earth is near the origin in both plots and moon is fixed at x=384000 in the rotating view and is not shown in the inertial view (each of the little loops in the inertial view is a lunar flyby)

Offline C5C6

  • Regular
  • Full Member
  • **
  • Posts: 264
  • Liked: 1
  • Córdoba - Argentina
    • programaespacial.com
Re: Simulate orbit using Matlab
« Reply #14 on: 12/01/2008 01:39 PM »
Great pictures todd!! unfortunately I didn't understand the code much =S

I'm now stuck with 3d graphics...I made vectors, but im not sure this is what you meant by vectors jorge... check it out...

the problem is here i'm sure:
   distancia = sqrt((SAT(1,2)^2)+(SAT(1,1)^2)+(SAT(1,3)^2));

        gravedad  = -mu_tierra/distancia^2;

   anguloXZ = atan(SAT(1,3)/SAT(1,1));
   anguloXY = atan(SAT(1,2)/SAT(1,1));

   acc_tierra_x = gravedad*abs(cos(anguloXZ));
   acc_tierra_y = gravedad*abs(cos(anguloXY));
   acc_tierra_z = gravedad*abs(sin(anguloXZ));

both X and Z seems to work just fine but Y is higher then what it should.....see that i calculated 2 angles, I tried changing the angles and always fails the acceleration wich uses the different angle: here, as x and z uses the same angle, y fails
« Last Edit: 12/01/2008 04:31 PM by C5C6 »