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

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 - I can tell why they let you teach this stuff to the astronauts.

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!!!

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

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

Quote from: C5C6 on 11/25/2008 01:15 AMThanks 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!!!!!!SergioYour 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:

Great pictures todd!! unfortunately I didn't understand the code much =SI'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

pos_sat_x = radio_tierra+distancia; % mpos_sat_y = 0; % mvel_sat_x = 0; % m/svel_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;

pos_sat_x = radio_tierra+distancia; % mpos_sat_y = 0; % mvel_sat_x = 0; % m/svel_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;

pos_sat_x = radio_tierra+distancia; % mpos_sat_y = 0; % mpos_sat_z = 0; % mvel_sat_x = 0; % m/svel_sat_y = velocidad; % m/svel_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;

pos_sat(1) = radio_tierra+distancia; % mpos_sat(2) = 0; % mpos_sat(3) = 0; % mvel_sat(1) = 0; % m/svel_sat(2) = velocidad; % m/svel_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;

pos_sat = [radio_tierra+distancia, 0, 0]; % mvel_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;

pos_sat(1) = radio_tierra+distancia; % mpos_sat(2) = 0; % mpos_sat(3) = 0; % mvel_sat(1) = 0; % m/svel_sat(2) = velocidad; % m/svel_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;

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.

f = odefun(t,y)

y0 = [radio_tierra+distancia; 0; 0; 0; velocidad; 0];

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

global mu = G*(masa_tierra+masa_satelite);

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

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