### Author Topic: Simulate orbit using Matlab  (Read 54673 times)

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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);
angulo = atan(pos_sat_y/pos_sat_x);
(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!!!

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### 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);
angulo = atan(pos_sat_y/pos_sat_x);
(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:

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:

(note that distancia squared becomes distancia cubed.)
JRF

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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...

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### 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:

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

#### Herb Schaltegger

##### 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.

#### RedSky

• Veteran
• Full Member
• Posts: 451
• Liked: 5
• Likes Given: 0
##### 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/

#### m330533

• STS Fanatic
• Member
• Posts: 66
• Duesseldorf, Germany
• Liked: 0
• Likes Given: 0
##### 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 !

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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??
Thid they made a burn to enter moon orbit??

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

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### 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 thank you so much for your help!!!

JRF

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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!!!!

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### 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_y   = 0;               % m

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

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

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### 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

#### toddbronco2

• Member
• Full Member
• Posts: 268
• Liked: 0
• Likes Given: 0
##### 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)

#### C5C6

• Regular
• Full Member
• Posts: 264
• Córdoba - Argentina
• Liked: 1
• Likes Given: 0
##### 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));

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

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 »

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### Re: Simulate orbit using Matlab
« Reply #15 on: 12/01/2008 07:54 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:

Whoa... let's try not to scare him off! I am belatedly realizing that he may not be familiar with the terminology we're throwing around. A step-by-step approach starting from his current code may be a better route. I'll try developing that over multiple posts.
« Last Edit: 12/02/2008 01:44 AM by Jorge »
JRF

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### Re: Simulate orbit using Matlab
« Reply #16 on: 12/01/2008 08:18 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));

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

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

Let's back up a bit. Your original code (with the graphics stripped out) looked like this:

Quote
pos_sat_y   = 0;               % m

vel_sat_x    = 0;               % m/s

...

pos_sat_x = pos_sat_x + vel_sat_x*escala_tiempo;
pos_sat_y = pos_sat_y + vel_sat_y*escala_tiempo;

distancia = sqrt(pos_sat_y^2+pos_sat_x^2);
angulo = atan(pos_sat_y/pos_sat_x);

if(pos_sat_x<0)
acc_tierra_x=-acc_tierra_x;
end
if(pos_sat_y<0)
acc_tierra_y=-acc_tierra_y;
end

vel_sat_x  = vel_sat_x + acc_tierra_x*escala_tiempo;
vel_sat_y  = vel_sat_y + acc_tierra_y*escala_tiempo;

and after the first round of fixes (gravitational acceleration vs force, use of trig identities to eliminate some unnecessary sin() and cos() calls), looked like this:

Quote
pos_sat_y   = 0;               % m

vel_sat_x    = 0;               % m/s

...

pos_sat_x = pos_sat_x + vel_sat_x*escala_tiempo;
pos_sat_y = pos_sat_y + vel_sat_y*escala_tiempo;

distancia = sqrt(pos_sat_x^2+pos_sat_y^2);

vel_sat_x  = vel_sat_x + acc_tierra_x*escala_tiempo;
vel_sat_y  = vel_sat_y + acc_tierra_y*escala_tiempo;

We can extend this to three dimensions by adding a z component, and for right now we can keep things simple by leaving z=0. Note that by sticking with cartesian components, you avoid the need for messy trig functions:

Quote
pos_sat_y   = 0;               % m
pos_sat_z   = 0;               % m

vel_sat_x    = 0;               % m/s
vel_sat_z    = 0;               % m/s

...

pos_sat_x = pos_sat_x + vel_sat_x*escala_tiempo;
pos_sat_y = pos_sat_y + vel_sat_y*escala_tiempo;
pos_sat_z = pos_sat_z + vel_sat_z*escala_tiempo;

distancia = sqrt(pos_sat_x^2+pos_sat_y^2+pos_sat_z^2);

vel_sat_x  = vel_sat_x + acc_tierra_x*escala_tiempo;
vel_sat_y  = vel_sat_y + acc_tierra_y*escala_tiempo;
vel_sat_z  = vel_sat_z + acc_tierra_z*escala_tiempo;

Now, one thing to point out: You don't need to be scared of switching to vectors, because you're already using them - you're just not taking advantage of them in Matlab. pos_sat, vel_sat, and acc_tierra are all vectors already so you just need to use subscripts rather than individual variables, like this:

Quote
pos_sat(2)   = 0;               % m
pos_sat(3)   = 0;               % m

vel_sat(1)    = 0;               % m/s
vel_sat(3)    = 0;               % m/s

...

pos_sat(1) = pos_sat(1) + vel_sat(1)*escala_tiempo;
pos_sat(2) = pos_sat(2) + vel_sat(2)*escala_tiempo;
pos_sat(3) = pos_sat(3) + vel_sat(3)*escala_tiempo;

distancia = sqrt(pos_sat(1)^2+pos_sat(2)^2+pos_sat(3)^2);

vel_sat(1)  = vel_sat(1) + acc_tierra(1)*escala_tiempo;
vel_sat(2)  = vel_sat(2) + acc_tierra(2)*escala_tiempo;
vel_sat(3)  = vel_sat(3) + acc_tierra(3)*escala_tiempo;

Now, you may look at the above code and think this really isn't saving you anything, and you'd be correct. The real savings comes from using Matlab's vector operations. It will probably speed up your code quite a bit since Matlab can use vector libraries that are optimized for your machine:

Quote
pos_sat   = [radio_tierra+distancia, 0, 0];      % m
vel_sat    = [0, velocidad, 0];               % m/s

...

pos_sat = pos_sat + vel_sat*escala_tiempo;
distancia = sqrt(pos_sat*pos_sat');
vel_sat  = vel_sat + acc_tierra*escala_tiempo;

Now the code is a lot shorter and probably runs a lot faster. But we can make it faster still, and more accurate... but I'll tackle that in a future post.
« Last Edit: 12/02/2008 01:48 AM by Jorge »
JRF

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### Re: Simulate orbit using Matlab
« Reply #17 on: 12/01/2008 08:50 PM »
Let's back up one step and look at the vectorized code, but using the individual components for clarity. I've rearranged it a bit to emphasize certain relationships:

Quote
pos_sat(2)   = 0;               % m
pos_sat(3)   = 0;               % m
vel_sat(1)    = 0;               % m/s
vel_sat(3)    = 0;               % m/s

...

distancia = sqrt(pos_sat(1)^2+pos_sat(2)^2+pos_sat(3)^2);

pos_sat(1) = pos_sat(1) + vel_sat(1)*escala_tiempo;
pos_sat(2) = pos_sat(2) + vel_sat(2)*escala_tiempo;
pos_sat(3) = pos_sat(3) + vel_sat(3)*escala_tiempo;
vel_sat(1)  = vel_sat(1) + acc_tierra(1)*escala_tiempo;
vel_sat(2)  = vel_sat(2) + acc_tierra(2)*escala_tiempo;
vel_sat(3)  = vel_sat(3) + acc_tierra(3)*escala_tiempo;

Check out the last six lines. The position and velocity vectors, pos_sat and vel_sat on the left hand side, can be treated as a single vector. We call this the "state vector". On the right hand side, vel_sat in the first three lines is the first derivative of pos_sat with respect to time, and acc_tierra in the last three lines is the first derivative of vel_sat with respect to time.

So, what we really have here is a system of six first-order ordinary differential equations (ODEs). By multiplying the derivative by the time step (escala_tiempo), and adding it back to the state, you're computing a numerical solution to this system of ODEs. Specifically, the method you used is Euler's first-order method.

Now, there's no reason to be intimidated by all the terminology. After all, you figured out the first-order solution yourself. That's good! But Euler's method has a weakness: it is only accurate if the derivative (vel_sat and acc_tierra) is reasonably constant over the time step (escala_tiempo). Obviously in this case, it isn't: velocity and acceleration are constantly changing throughout the orbit. We can improve accuracy by using shorter time steps, but ultimately if we want an accurate orbit we need to use a higher-order method.

There are many methods for the numerical solution to ODEs, and Matlab code for many of them are out there on the web. For example:

http://math.bu.edu/people/josic/classes/MA226/matlab/rk4.m

But it's not even necessary to use outside code, since Matlab has powerful ODE solvers built in. However, to use them, you need to write your code in a form that the ODE solvers can work with. Specifically, you need to express the state vector as a single vector, and you need to put the computation of the derivative into a separate function that the ODE solver can call. I'll discuss that in the next post.
JRF

#### Jorge

• Senior Member
• Posts: 6177
• Liked: 7
• Likes Given: 0
##### Re: Simulate orbit using Matlab
« Reply #18 on: 12/02/2008 01:40 AM »

But it's not even necessary to use outside code, since Matlab has powerful ODE solvers built in. However, to use them, you need to write your code in a form that the ODE solvers can work with. Specifically, you need to express the state vector as a single vector, and you need to put the computation of the derivative into a separate function that the ODE solver can call. I'll discuss that in the next post.

Matlab's ODE solvers are documented here:

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode23.html

Specifically, it requires that the derivative function take the form:

Quote
f = odefun(t,y)

where y is the state vector (a column vector), t is the time span for the integration, and f is the derivative of the state vector (also a column vector).

For our simple two-body problem, we can define the state vector like this:

Quote

Note that this is just pos_sat and vel_sat jammed together, with the elements separated by semicolons rather than commas to make Matlab treat it as a column vector.

We can define the derivative function like this:

Quote
function ydot = twobody (t, y)
global mu;
pos_sat   = y(1:3);
distancia   = sqrt(pos_sat'*pos_sat);
end

Since we need to put this in a separate file (twobody.m), we also need to declare mu as a global variable and define it like this in the main script:

Quote
global mu = G*(masa_tierra+masa_satelite);

When we're ready to actually propagate the state vector, we call ode45 like this:

Quote
[t,y]   = ode45(@twobody, 0:escala_tiempo:tf, y0);

In the above, tf is the stop time for the simulation. For Apollo 17, TLI cutoff was at 3:18:28 and LOI ignition was at 86:14:23, or a travel time of 82:55:55, or 298555 seconds. t is a column vector of output times and can be discarded. y contains the output state vectors; each row corresponds to a row in t and the columns correspond to the x,y,z position and x,y,z velocity components, respectively. Note that now you don't even need a main loop for the program; ode45 loops over the time span for you. You can plot the positions in 3D like this:

Quote
plot3(x(:,1), x(:,2), x(:,3))

Attached is a 3D plot of the 2-body Apollo 17 trajectory, using a 10-minute escala_tiempo. Once you've gotten more comfortable with this, you can start playing around with toddbronco2's restricted 3-body code - it looks ready to "plug and play" with ode45.
« Last Edit: 12/02/2008 03:40 AM by Jorge »
JRF

#### jurgen

• Member
• Posts: 13
• Liked: 0
• Likes Given: 0
##### Re: Simulate orbit using Matlab
« Reply #19 on: 12/02/2008 04:05 AM »
I did a project similiar to this in my freshman year.  The project mainly delt with planetary bodies in binary star systems.  I'm going to see if I can clean up the code and convert it to use matlabs built in vector functions to simplify the program.  I'll post it later this week.