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))
=) 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.
pos_sat_x = radio_tierra+distancia; % mpos_sat_y = 0; % mvel_sat_x = velocidad*sin(gamma); % m/svel_sat_y = velocidad*cos(gamma); % m/s
Please check the other two 2D simulations of TLI, where I could make a nice and kinda-violent rendezvous =)
Here's my TLI attempt, enjoy!!
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.
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.
Quote from: Thiha Kyaw on 06/29/2014 11:44 amI 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_precessionSince 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.
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_precessionSince 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.