NASASpaceFlight.com Forum

General Discussion => General Discussion => Topic started by: C5C6 on 11/22/2008 06:37 pm

Title: Simulate orbit using Matlab
Post by: C5C6 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!!!

Title: Re: Simulate orbit using Matlab
Post by: Jorge 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.)
Title: Re: Simulate orbit using Matlab
Post by: C5C6 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...
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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.
Title: Re: Simulate orbit using Matlab
Post by: Herb Schaltegger on 11/23/2008 01:37 am
Jorge -

I can tell why they let you teach this stuff to the astronauts. :)
Title: Re: Simulate orbit using Matlab
Post by: RedSky 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/
Title: Re: Simulate orbit using Matlab
Post by: m330533 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 !
Title: Re: Simulate orbit using Matlab
Post by: C5C6 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!!!
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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!
Title: Re: Simulate orbit using Matlab
Post by: C5C6 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!!!!
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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.
Title: Re: Simulate orbit using Matlab
Post by: C5C6 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
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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.
Title: Re: Simulate orbit using Matlab
Post by: toddbronco2 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)
Title: Re: Simulate orbit using Matlab
Post by: C5C6 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
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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.
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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));

        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

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

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

vel_sat_x    = 0;               % m/s
vel_sat_y   = velocidad;            % 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);
   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));

   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_x   = radio_tierra+distancia;      % m
pos_sat_y   = 0;               % m

vel_sat_x    = 0;               % m/s
vel_sat_y   = velocidad;            % 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);

   gravedad  = -((G*(masa_tierra+masa_satelite))/(distancia^3));

   acc_tierra_x = gravedad*pos_sat_x;
   acc_tierra_y = gravedad*pos_sat_y;

   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_x   = radio_tierra+distancia;      % m
pos_sat_y   = 0;               % m
pos_sat_z   = 0;               % m

vel_sat_x    = 0;               % m/s
vel_sat_y   = velocidad;            % 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);

   gravedad  = -((G*(masa_tierra+masa_satelite))/(distancia^3));

   acc_tierra_x = gravedad*pos_sat_x;
   acc_tierra_y = gravedad*pos_sat_y;
   acc_tierra_z = gravedad*pos_sat_z;

   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(1)   = radio_tierra+distancia;      % m
pos_sat(2)   = 0;               % m
pos_sat(3)   = 0;               % m

vel_sat(1)    = 0;               % m/s
vel_sat(2)   = velocidad;            % 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);

   gravedad  = -((G*(masa_tierra+masa_satelite))/(distancia^3));

   acc_tierra(1) = gravedad*pos_sat(1);
   acc_tierra(2) = gravedad*pos_sat(2);
   acc_tierra(3) = gravedad*pos_sat(3);

   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');
   gravedad  = -((G*(masa_tierra+masa_satelite))/(distancia^3));
   acc_tierra = gravedad*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.
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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(1)   = radio_tierra+distancia;      % m
pos_sat(2)   = 0;               % m
pos_sat(3)   = 0;               % m
vel_sat(1)    = 0;               % m/s
vel_sat(2)   = velocidad;            % m/s
vel_sat(3)    = 0;               % m/s

...

   distancia = sqrt(pos_sat(1)^2+pos_sat(2)^2+pos_sat(3)^2);
   gravedad  = -((G*(masa_tierra+masa_satelite))/(distancia^3));
   acc_tierra(1) = gravedad*pos_sat(1);
   acc_tierra(2) = gravedad*pos_sat(2);
   acc_tierra(3) = gravedad*pos_sat(3);

   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.
Title: Re: Simulate orbit using Matlab
Post by: Jorge 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
y0   = [radio_tierra+distancia; 0; 0; 0; velocidad; 0];

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);
   gravedad   = -mu/distancia^3;
   ydot   = [y(4:6); gravedad*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.
Title: Re: Simulate orbit using Matlab
Post by: jurgen 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.
Title: Re: Simulate orbit using Matlab
Post by: C5C6 on 12/02/2008 06:34 pm
Jorge thanks for everything!!! you say that Euler's first isn't accurate if acceleration varies over time, but my simulations are producing beautiful and accurate-like graphs...take a look:
Title: Re: Simulate orbit using Matlab
Post by: Jorge on 12/04/2008 06:12 pm
Good work! Euler's method will work fine for generating pretty pictures. You can pick up the higher-order methods if you decide later that you need more accuracy.
Title: Re: Simulate orbit using Matlab
Post by: robertross on 12/04/2008 08:29 pm
I want to say to both of you, that following along this thread has really shown the true nature of cooperation on this site. I commend Jorge for his excellent help and step-by-step approach. Maybe we can all take a page from this and be better for it to all the other posters. :)

Thanks for the priveledge for following along!
Title: Re: Simulate orbit using Matlab
Post by: wjbarnett on 12/04/2008 08:34 pm
What Robert said! Here's to Jorge and his kind, step by step explanations! I don't have Matlab, but I learned a lot just by reading. Thanks!
Title: Re: Simulate orbit using Matlab
Post by: C5C6 on 12/04/2008 09:15 pm
=) yeap nice cooperation!!!

I got a question thinkng about higher-order euler method: you said that first order method wasn't that accurrate due to huge variations in the velocity and accelerations. I believe you refer to the case of using high escala_tiempo: in the first simulations I used high escala_tiempo to accellerate the simulation process, but I found out that it produced huge differences so I reduced it to 1, and that actually 'solved the problem' because the orbits looked much like the expected ones.

The graphics I posted were made with an escala_tiempo 1, so the simulations used approximately as many cicles as seconds one orbit period lasted (95*60 for ISS, 60*60*24*29 for moon).

I want to make sure I understand the necessity of higher euler orders: using ODEs, would I be simulating as if I used an infinitesimal escala_tiempo??

----------------------------------

Now some doubts about TLI: we have discussed about what Flight-path angle was, but I'm afraid I got it wrong so I'm attaching an image, and I marked as beta what I believe that angle is. I'm also not sure how to introduce it in the simulation: what I did to the post TLI burn angle was to multiply the post-TLI velocity factors x and y with sin or cos of that angle. Please check the other two 2D simulations of TLI, where I could make a nice and kinda-violent rendezvous =)

THANKS!!
Title: Re: Simulate orbit using Matlab
Post by: C5C6 on 12/05/2008 12:56 am
Here's my TLI attempt, enjoy!!
Title: Re: Simulate orbit using Matlab
Post by: Jorge on 12/05/2008 03:20 am
=) yeap nice cooperation!!!

I got a question thinkng about higher-order euler method: you said that first order method wasn't that accurrate due to huge variations in the velocity and accelerations. I believe you refer to the case of using high escala_tiempo: in the first simulations I used high escala_tiempo to accellerate the simulation process, but I found out that it produced huge differences so I reduced it to 1, and that actually 'solved the problem' because the orbits looked much like the expected ones.

The graphics I posted were made with an escala_tiempo 1, so the simulations used approximately as many cicles as seconds one orbit period lasted (95*60 for ISS, 60*60*24*29 for moon).

I want to make sure I understand the necessity of higher euler orders: using ODEs, would I be simulating as if I used an infinitesimal escala_tiempo??

Kinda. With a higher order method, you can achieve much higher accuracy at the same step size, or you can achieve the same accuracy at much larger step sizes. And the built-in ODE solvers are highly optimized; even though they're more complex than a first-order Euler method, they will actually run faster at the same step size (Don't forget - Matlab is *interpreting* your hand-coded Euler integrator, but the built-in ODE solvers are *compiled* code).

And a note on terminology - not all methods come from Euler. Matlab's ode45 and ode23 solvers are based on Runge-Kutta methods, for example.

Quote
Now some doubts about TLI: we have discussed about what Flight-path angle was, but I'm afraid I got it wrong so I'm attaching an image, and I marked as beta what I believe that angle is.

It's not obvious from that picture. Flight path angle (gamma) is the angle between the velocity vector and the local horizontal. The local horizontal is, in turn, always perpendicular to the position vector. So gamma=0 means the velocity vector is perpendicular to position, and gamma=90 means the velocity vector is along the position vector (i.e. straight up).

Quote
I'm also not sure how to introduce it in the simulation: what I did to the post TLI burn angle was to multiply the post-TLI velocity factors x and y with sin or cos of that angle.

That's right. If you're putting the initial position all in x, the velocity should look like this:

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

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

(Make sure that the flight path angle, gamma, is in radians.)

Your TLI2.bmp picture looks correct for the Apollo 17 trajectory numbers I gave you.

If you want to go 3-D with the Apollo 17 post-TLI trajectory, you'll need to convert all six parameters from the mission report (latitude, longitude, altitude, speed, flight path angle, and azimuth) into cartesian position and velocity vectors. It's a bit complicated to explain so I'll just attach my code.

Quote
Please check the other two 2D simulations of TLI, where I could make a nice and kinda-violent rendezvous =)
Title: Re: Simulate orbit using Matlab
Post by: Jorge on 12/05/2008 03:27 am
Here's my TLI attempt, enjoy!!

It looks like it's going around the correct side of the moon - now you just need to find a way to slow it down once you get there. :)
Title: Re: Simulate orbit using Matlab
Post by: jurgen on 12/05/2008 04:41 am
One thing you can notice for a Euler method is that the error increases as the two planetary bodies get closer to each other.  With a fixed time step this can cause one of the bodies to move farther than it should.  So one method that you can use to increase accuracy is to implement a variable timestep that is a function of position between the two bodies.
Title: Re: Simulate orbit using Matlab
Post by: Apollo-phill on 12/26/2008 10:03 am
Hi

Get a copy of Donald C. Lundy book called "Lunar Landing And Return" - a simplified physics and mathematical investigation of the "Apollo-11 Saga".

It goes through the whole process of a Saturn V - Apollo lunar landing mission from launch to splashdown and everything in between and uses the same maths (simplified albeit) that was used in a typical Apollo mission.

ISBN 0-75961-858-5


It is a large format (A4 size) paperback.

It steps you through the calculations required for each major stage of the mission and will answer most of the questions the original enquirer wants.

Additionally, you could read my article on the Apollo onboard computers that appear in the Apollo Lunar Surface Journal :-

http://history.nasa.gov/afj/compessay.htm

and my Apollo microcircuit article at ALSJ at:-

http://www.hq.nasa.gov/office/pao/History/alsj/apollo-ic.html


Phill Parker
UK
Title: Re: Simulate orbit using Matlab
Post by: Apollo-phill on 12/26/2008 11:23 am
Hi

Ignore my previous entry - it was menat for another topic in the forum where I person wanted simplified details of Apollo type missions

Apologies - too much christmas pud !

Phill
Title: Re: Simulate orbit using Matlab
Post by: Marya on 04/10/2009 08:48 am
Great!

Very nice conversation about Orbit Simulation using MATLAB.

Thanks so much Jorge and C5C6.

marya
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/23/2009 03:20 am
First off, let me say thanks for reading and looking into my post.

Anyone who is interested can download and play with these files, but I am having a few problems getting this code from my teacher to run.  It is an older matlab file dealing with a simplified 2-D static orbit simulation problem.

I would like to see the results and charts, but I cannot get it to run completely. 

I keep getting an error saying no default parameters set.  I will upload all the files from my teacher.

The M file I am trying to run is 'orbit.m' with different cases to select.

Can anyone please help me figure what to add or change?

I put all the files in the same folder to run them so I did not have to change the directory.

*UPDATE*  to further clarify, I believe you also need orbitode, because in the orbit m file it is called upon....I think.  I see some people have downloaded it, anyone have any ideas yet?  I may need to change the parameters in orbitode?  Not sure though.
Thanks in advance everyone.
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/23/2009 06:26 pm
*UPDATE*

I believe this is the line that needs to be changed.
[t,y,te,ye,ie] = ode45( 'orbitode', [0; tend], [], options, y0, mu, fric );

More specifically I think my problem is  the first part of ODE45, the 'orbitode'

My question.  How do I modify the 'orbitode.m' that now works, to take the values set from 'orbit.m' and plug them in, differeniate, and then plot. 

I think I have two separate programs.  I need to modify the function in 'orbitode.m' to reflect orbit.m    So how do I modify orbitode so that when I use the above code in 'orbit.m' it works?


*********************************************************

I figured out the problem with 'orbitode.m'  Change odeibm or whatever to ode45, this then plots the figure eight.  However, I am still unable to get the file I want to work 'orbitm'  What do I need to change and/or add?

Anyone?


First off, let me say thanks for reading and looking into my post.

Anyone who is interested can download and play with these files, but I am having a few problems getting this code from my teacher to run.  It is an older matlab file dealing with a simplified 2-D static orbit simulation problem.

I would like to see the results and charts, but I cannot get it to run completely. 

I keep getting an error saying no default parameters set.  I will upload all the files from my teacher.

The M file I am trying to run is 'orbit.m' with different cases to select.

Can anyone please help me figure what to add or change?

I put all the files in the same folder to run them so I did not have to change the directory.

*UPDATE*  to further clarify, I believe you also need orbitode, because in the orbit m file it is called upon....I think.  I see some people have downloaded it, anyone have any ideas yet?  I may need to change the parameters in orbitode?  Not sure though.
Thanks in advance everyone.
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/24/2009 08:23 pm
Ok.  So I decided to try and do my own version.  So I started a new file that is not complete but should be written well enough for everyone to follow.  I will also upload a project description of what I am trying to do. 

If anyone is willing to help me get this running and working I would be so appreciative.  Not only is this an extra credit project, it is also something that I generally take interest in since I am an aero major. 

Thanks everyone for looking and helping.  Please post some comments.

I know this is a over simplified model, but it is still fun to work on.
Title: Re: Simulate orbit using Matlab
Post by: toddbronco2 on 04/25/2009 03:31 pm
I tried running the codes and I get the same error as you.  I've never seen ode45 fed arguments like that so I don't really know what the script is trying to do.  I'm sorry I'm not more help on that, but there are few things less fun than trying to debug code written by somebody else.

On the dynamics side though, your assignment is to simulate the "stationary" three body problem, right?  The differential equations that you've got in your "orbitode" script are those of the Circular Restricted Three Body Problem.  You'll need to modify those differential equations be removing the terms relating the centripetal and Coriolis accelerations relating to the Earth and Moon's circular motion around their common barycenter.  That's actually not too hard to do.  I think those are just the first two terms in both your dydt(3) and dydt(4) terms.
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/25/2009 05:11 pm
Thanks.  I decided to start a new file, which are those new codes attached to my above post. 

I know you and everyone else is busy with there own work, but can someone help me finish writing this code seeing as I am not very experienced with matlab code.

I asked my teacher and he said that I needed to:

"put the If-statement for the thrust in the function where it's used as part of the RHS. Also re-think the direction of this force [the thrust force in the ode file]. You need a unit vector in the tangential direction, how do you best construct this?"

After this, I think I pretty much have it.  Thanks for your help.

Here are the two files once again with my comments.
Title: Re: Simulate orbit using Matlab
Post by: prince1 on 04/25/2009 07:44 pm
hello guys,

can anyone help me to get the matlab files that simulate Molniya Orbit?
(3D simulation and the ground track)

i would really appreciat that.
Title: Re: Simulate orbit using Matlab
Post by: C5C6 on 04/25/2009 07:56 pm
Here you my code that simulates a molniya....it doesn't have a ground track, you can see the earth and the orbit...
Title: Re: Simulate orbit using Matlab
Post by: prince1 on 04/25/2009 10:30 pm
thank you soooooooooooooooo much C5C6
 i can see the earth and the orbit, that's awesome

i am still trying to get the ground track
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/25/2009 11:51 pm
So,  after looking at the files where should I place the if statement exactly?  I have my function but not sure where to place it or how it should be written. If everyone would prefer I can just paste the code instead of uploading the files.
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/26/2009 10:10 pm
driverss3matt.m

------
% Script Solution for 3-body problem (earth, moon, spaceship)
% 2-D
% Spaceship thruster is fired for a certain period of time from a...
% geostationary orbit
% Figure eight solution

clear;
close all;
clc;
global me mm ms G xm thrust thrustduration

% -----Constants------%
            me = 59742e24;  % mass of the earth [kg]
            mm = 7.35e22;   % mass of the moon [kg]
            ms = 1.5e5;     % mass of the spaceship [kg]
            G  = 6.673e-11; % Gravitational Constant [m^3 kg^-1 s^-2]
            xm = 384e6;     % distance from the earth to the moon [m]

% -----Parameters------%
thrustduration = 452;   
        tfinal = 5e4;   
thrust = 888e5; % amount of thrust [lb]

            x0 =  -6.578e6 ; % initial position in x
            y0 =  0        ; % initial position in y
            v0 = sqrt((G*me)/(6.4e6+2e5));    % initial velocity
           vx0 = 0; % initial velocity in x
           vy0 = sqrt((G*me)/(6.4e6+2e5))+(thrustduration*thrust)/ms; % initial velocity in y

%              t =  ; % time
%       t_begin =  0; % time for initiate thrust?
%       t_end   =  10;
%       t_travel =  ; % time for spaceship to travel

         tspan =  [0,tfinal] ; % time span, entire?
%          tsize =  size(tspan); % time size?
         
%----If Statement for Thrust----%
       

            y0 = [x0,y0,vx0,vy0]; % not sure what all needs to be here

         [t,y] = ode45('mattscode', tspan, y0);

             x = y(:,1) ; % position in x after time t
             y = y(:,2) ; % position in y after time t
%           vx1 = y(:,3) ; % velocity in x after time t
 %          vy1 = y(:,4) ; % velocity in y after time t

%----Plotting Velocity Analysis----%
% figure(1);
% subplot(2,1,1), plot(t,vx1);
% subplot(2,1,2), plot(t,vy1);
% xlabel('Time (sec.)');
% ylabel('Velocity (m/sec.)');


%----Plotting Earth and Moon----%

%Earth%
figure(2);
h=0; k=0; r=6.4e6; N=256;
t=(0:N)*2*pi/N;
fill( r*cos(t)+h, r*sin(t)+k,'r');
text(0,0,'E')
axis('equal')
hold on

%Moon%
h=xm; k=0; r=1.75e6; N=256;
t=(0:N)*2*pi/N;
fill( r*cos(t)+h, r*sin(t)+k,'b');
text(xm,0,'M')
axis('equal')

%----Plotting Orbit----%
plot(x,y,'k');
title('Spaceship Orbit');
ylabel('y(t)');
xlabel('x(t)');
axis('equal')
hold off

-------

mattscode.m

--------function ydot = mattscode(t,var)

global me mm ms G xm thrust thrustduration

x  = var(1);    % x position [m]
y  = var(2);    % y position [m]
vx = var(3);    % x velocity [m/s]
vy = var(4);    % y velocity [m/s]

% -----2-D problem => 4 ODE's------%
if t <= thrustduration
   
yd1 = vx;
yd2 = vy;
yd3 = 1/ms*(-G*me*ms/(x^2+y^2)*x/sqrt(x^2+y^2) + G*mm*ms/((xm-x)^2+y^2)*x/sqrt((xm-x)^2+y^2) + thrust*x/sqrt(x^2+y^2)); % unit vectors instead of angles
yd4 = 1/ms*(-G*me*ms/(x^2+y^2)*y/sqrt(x^2+y^2) - G*mm*ms/((xm-x)^2+y^2)*y/sqrt((xm-x)^2+y^2) + thrust*y/sqrt(x^2+y^2)); % unit vectors instead of angles
ydot = [yd1;yd2;yd3;yd4];
else
yd1 = vx;
yd2 = vy;
yd3 = 1/ms*(-G*me*ms/(x^2+y^2)*x/sqrt(x^2+y^2) + G*mm*ms/((xm-x)^2+y^2)*x/sqrt((xm-x)^2+y^2)); % unit vectors instead of angles
yd4 = 1/ms*(-G*me*ms/(x^2+y^2)*y/sqrt(x^2+y^2) - G*mm*ms/((xm-x)^2+y^2)*y/sqrt((xm-x)^2+y^2)); % unit vectors instead of angles

ydot = [yd1;yd2;yd3;yd4];
end



end
----------

if you run this as is, you will see that a spaceship goes around the earth and hits the moon.

there is something not taking into account in the equations so that gravity of moon affects spaceship.

Also, I think something about the parameters may be wrong.  Can anyone help or spot the problem.

**From last post, my teacher said to put the thrust as a tangential vector, so how do I update that from what I have now**

Sorry for the strikethough, some of it is commented out, but some of it is a comment made or not finished.
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 04/27/2009 02:42 am
I am so close now...if and when i get it I will upload.

The last thing is that dang thrust.

My teacher told me: You need a vector which always points in the tangential direction. Then you need to build a unit vector from this ...

I am assuming its the n-t coordinate system, but I am unsure how to formulate this into what I need for matlab. 

Any help?
Title: Re: Simulate orbit using Matlab
Post by: thundergod on 07/08/2009 06:24 pm
So I'm back.

I managed to finish and get most of the credit for my project.  Thanks everyone for your help.

My new personal project is to make this a 3d simulation such as an earlier member was doing.

Anyone want to help?

Right now I only have the figure and earth/moon plotted, but I am working on the code for the orbit.

I have included my final 2d version--it still crashes into the moon before making the figure eight, maybe someone could figure out why for me

Also the beginning of my 3d version.
Title: Re: Simulate orbit using Matlab
Post by: ersin2323 on 11/02/2009 03:31 pm
any one has standart LEO Matlab Codes? I need a simulation...
Title: Re: Simulate orbit using Matlab
Post by: Art_UA on 07/03/2010 03:37 pm
Hello!
Sorry for bumping such old topic. I REALLY need some help.
I'm trying to make a 2D model of restricted three-body problem, so my problem is relative to problems being discussed in this topic and I decided to reply here because of this.
My goal is to make a model of figure eight transfer between the Earth and the Moon. First at all, I decided to make a model of the static Earth and Moon moving around and then add a spaceship to this model. I have changed the first C5C6's program with Jorge's correctives, but what I got doen't actually looks like Moon's orbit. If someone will help me with this I could be very much obliged.
I've added a model of figure eight transfer between the Earth and the Moon with static Earth and Moon (from www.mathworks.com). Maybe modifying this model to make Moon move will be easier than step-by-step C5C6's program modifying?
Title: Re: Simulate orbit using Matlab
Post by: Art_UA on 07/04/2010 12:54 pm
I've added a satellite to my model, but it does't work correctly. There are some troubles witn calculating of the gravitational acceleration I think. I calculate it like this:

km = (pos_sat_y - pos_Moon_y)/(pos_sat_x - pos_Moon_x); % angle coefficient of vector between moon & satellite
    ke = (pos_sat_y)/(pos_sat_x); % angle coefficient of vector between earth & satellite
    cosalpha = abs(km*ke+1)/((1+km^2)^1/2)*(1+ke)^(1/2); %cos of angle between vectors form satellite to moon & from satellite to earth
    se = (G*Earth_mass)/(pos_sat_x^2+pos_sat_y^2)^(1/2); %vector between satellit & earth
    sm = ((G*mass_Moon)/((pos_sat_x - pos_Moon_x)^2 + (pos_sat_y - pos_Moon_y)^2)^(1/2)); %vect between satellite & moon
    grav_acc_sat = se^2 + sm^2 + 2*se*sm*cosalpha;
    anglesat = atan(pos_sat_y/pos_sat_x);
   sat_acc_y = grav_acc_sat*abs(sin(anglesat));
   sat_acc_x = grav_acc_sat*abs(cos(anglesat));

How do you think is it correct?
Thanks for any help
Title: Re: Simulate orbit using Matlab
Post by: Thiha Kyaw on 06/29/2014 11:44 am
Hi all friend.
 I am a student studying space system engineering at Myanmar Aerospace Engineering University. I am interested in orbit. ACE orbit is special. The line of apsides ( line between apogee and perigee)rotates 360 degree per year. I want know how this happen. If you know, share me some information or link of free book.
And I am also trying to simulate an orbit and ground track for an nanosatellite of eccentricity(0.65) , inclination(19 deg), semimajor axis ( 40600km). I want know what  you consider to simulate an orbit( 2D or 3D). How to consider about rotation of earth. If you have some m-files for molniya. Thank you all.
Title: Re: Simulate orbit using Matlab
Post by: Thiha Kyaw on 07/01/2014 11:58 am
The following attachment M-file is to simulate the ground track of a satellite operating any orbit using Matlab. You can change classical orbital elements. These are semi-major axis (a), eccentricity (e), inclination (iota), argument of perigee (ome), right ascension of as ascending node (omega). You can also change number of period (lab) to get successive ground track per orbit.   You can run your orbit to simulate ground track of satellite. The M-file is written for an ellipse orbit to pass over Myanmar country as much duration as possible.  If you change the inputs as following, it will provide ground track of polar orbit.   a = 6878e3; e=0; iota=97.404*rad; omega=75.308*rad; ome=0*rad; lab=15.29;
Title: Re: Simulate orbit using Matlab
Post by: Thiha Kyaw on 07/01/2014 12:08 pm
polar ground track simulated by my upper post
Title: Re: Simulate orbit using Matlab
Post by: Andrewwski on 07/02/2014 03:33 am
Hi all friend.
 I am a student studying space system engineering at Myanmar Aerospace Engineering University. I am interested in orbit. ACE orbit is special. The line of apsides ( line between apogee and perigee)rotates 360 degree per year. I want know how this happen. If you know, share me some information or link of free book.

I assume from your description you're talking about an Apogee at Constant time of day Equatorial orbit (that's the only orbit that fits your description that ACE makes sense in).  This is just due to nodal precession.  Wikipedia explains it decently: http://en.wikipedia.org/wiki/Nodal_precession

Since you want one rotation of the line of apsides per year, you can calculate the rate of procession and plug that in, along with your other known variables, to solve for whatever you're looking for.

Quote
And I am also trying to simulate an orbit and ground track for an nanosatellite of eccentricity(0.65) , inclination(19 deg), semimajor axis ( 40600km). I want know what  you consider to simulate an orbit( 2D or 3D). How to consider about rotation of earth. If you have some m-files for molniya. Thank you all.

The code you've posted above does that.  If you want to create a ground track of a Molniya orbit, put in those parameters.

The general process would be to convert your position vector in Earth-Centered Inertial (ECI) coordinates to Earth-Centered Earth-Fixed (ECEF) coordinates.  This transformation is a function of the current time (thus accounting for the rotation of the earth).  Then, convert from ECEF to geodetic (latitude, longitude, elevation) coordinates.  MATLAB has a function for this (ecef2geodetic), although you could write your own if you want a better understanding how it works.
Title: Re: Simulate orbit using Matlab
Post by: Thiha Kyaw on 07/04/2014 10:29 am


I assume from your description you're talking about an Apogee at Constant time of day Equatorial orbit (that's the only orbit that fits your description that ACE makes sense in).  This is just due to nodal precession.  Wikipedia explains it decently: http://en.wikipedia.org/wiki/Nodal_precession

Since you want one rotation of the line of apsides per year, you can calculate the rate of procession and plug that in, along with your other known variables, to solve for whatever you're looking for.



The code you've posted above does that.  If you want to create a ground track of a Molniya orbit, put in those parameters.

The general process would be to convert your position vector in Earth-Centered Inertial (ECI) coordinates to Earth-Centered Earth-Fixed (ECEF) coordinates.  This transformation is a function of the current time (thus accounting for the rotation of the earth).  Then, convert from ECEF to geodetic (latitude, longitude, elevation) coordinates.  MATLAB has a function for this (ecef2geodetic), although you could write your own if you want a better understanding how it works.
Thank you Mr.Andrewwski. Now, I am trying to understand nodal precession  with the help of your link. And then, I am learning the build-in function (ecef2geodetic). But this is difficult to me.  If you don’t mind, and have some m-files that use the function (ecef2geodetic), please share me. Thanks for your time.
Title: Re: Simulate orbit using Matlab
Post by: Andrewwski on 07/07/2014 06:12 pm
I don't have any simple files offhand, but it's really quite simple to use.  Take a look at the MATLAB documentation for the function: http://www.mathworks.com/help/map/ref/ecef2geodetic.html

The syntax for usage is:

[phi,lambda,h] = ecef2geodetic(x,y,z,ellipsoid)

The inputs x, y, and z are the coordinates in the ECEF frame.  All axes are centered at the center of the earth.  The "x" axis intersects the point where the equator and prime meridian cross, the "z" axis goes through the north pole, and the "y" axis is orthogonal to x and z to form a right-handed coordinate system (so it goes through the point where the equator intersects the 90 degree east line of longitude).

Ellipsoid is a referenceEllipsoid that defines the shape of the earth.  You can just do something like ellipsoid=referenceEllipsoid('earth'); using MATLAB's catalog - although for the purposes of this ground track you likely could just create a sphere with the radius of the earth.

The outputs phi and lambda are latitude and longitude in radians, respectively, and h is height above the surface (which isn't used in creating the ground track).

So once you have the position in ECEF coordinates, it's rather easy to plot the ground track.  The harder part will be getting your time-series of positions in ECEF coordinates.  You'll likely have them in ECI coordinates which you'll have found from some propagation of an initial state.  There's a number of sources that a quick Google search turns up on this - this one is pretty straightforward: http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Appendix.pdf

There's also a function in the MATLAB file exchange that I've used before and confirmed works correctly.

The m-file you posted seems to be a more generalized implementation of this by creating the latitude/longitude as theoretical rotations of the earth and projecting the current position onto this, but I'd have to dive deeper into it to figure out exactly what it's doing.  It's rather cumbersome and difficult to decipher code that isn't commented.