phate

This is a great surprise to me. I thought JPL was able to know both position and velocity with high accuracy. Even because with such a problem about the initial velocities vector to propagate the state vector of a spacecraft AND the solar system bodies with high precision is not possible at all. I mada a rude check on the velocities and of course they match with those from the keplerian motion, but I can't check small errors... You are right. I use the standard MATLAB/SIMULINK rk4(5) integrator. I studied several other methods, but I have to think about how it could be possibile to "integrate" them in the SIMULINK environment. I don't know much about C++ or JAVA, do you think I should move to another programming language or try to implement another integrator in SIMULINK? I had the same idea before I started (even if I don't now anything about the Chebyshev polynomial coefficients, I thought to get the position of each planet from the ephemerides file in cartesian coordinates) but I moved to a more "dynamical" approach. Following the way you suggest, I should also take into account the possibility to insert a spacecraft in my scenario. In this case, I should get the position of the planets at each instant of time from the ephemerides, then calculate the perturbing force acting on the spacecraft to be added to the keplerian motion of the spacecraft around the primary of motion. Would it be correct? Do you suggest me to change the integrator in this case as well? Again, thank you very much for your patience and your willingness.

Hi everybody, in the last days I tried to manage with the problem. This is the situation to date. I used the equation for N-Body problem ([math]\frac{d^{2}\textbf{r}_{i}}{dt^{2}}=G\sum\frac{m_{j}}{r^{3}_{ij}}(\textbf{r}_{j}-\textbf{r}_{i}) [/math]) to integrate the state of the 9 planets of the solar system with initial condition from the JPL ephemerides tool HORIZONS. Then I compared the numerical data against the ephemerides data (assuming that JPL ephemerides are correct... :D). After a great effort I still have big differences between the predicted data and the real ones even after short integration times: I made some tests, and propagated the Earth "switching off" all the perturbing bodies, that is reducing the model to a simple two body problem, and the results are consistent with this simple model. Now my question is: which kind of perturbing actions JPL put into the ephemerides? I wonder if I am looking for a precision that is out of my possibility... The work is going ahead, I will keep you informed. In the meantime if you all have any kind of suggestion I will greatly appreciate. Regards

DH, so you are somehow saying that the drift I observed in my simulator is correct? I agree with the fact that the Sun acts as a perturbing force on the 2-body problem: my doubt is about the module of this drift, because my feeling is that approx 5000 km (compared to the keplerian orbit) in only 60 days is quite a big number. :confused: Do you have any number to compare? I thought to take the ephemerides of the Moon for one month and compare the propagated solution with the ephemerides, do you think it is a good check? Regards

Sisyphus, I think that using the N-body equations of motion with proper initial conditions should take into consideration the real eccentricities of the bodies. I agree, something is surely wrong... DH, do you suggest me to use the Lagrange's Planetary equations? As far as I know, the should give the rates for the orbital elements variations, but to be honest I am not so confident with them. I will have a look to the theory behind, in the meantime have you got any kind of suggestions? Thank you everybody for your help. Regards

Hi everybody, I would like to ask your opinion about some results from my solar sistem simulator. First, let me explain the method I used. I started with the classical N-Body problem and wrote the equations of motions in the usual form: [math]\frac{d^{2}\textbf{r}_{i}}{dt^{2}}=G\sum\frac{m_{j}}{r^{3}_{ij}}(\textbf{r}_{j}-\textbf{r}_{i})[/math] where, of course, the sum is valid for j = 1,..,N and j not equal to i. To semplify the implementation I integrate the equations of motion in the inertial reference frame centered on the Sun, being the initial conditions for velocity and position given by the JPL ephemerides software HORIZONS (even if my final target is to perform the integration in the Earth - Moon barycenter..). To date, my software only takes into account the presence of the Sun, the Earth, The Moon and Jupiter and in order to validate my results I focused my attention on the motion of the Moon around the Earth center. After a short time integration (60 days) I noticed a drift of the Moon orbit of about 5000 km approximately in the direction of the Sun. Now, I am quite surprised about that behavior and I would like to ask you if the gravitational effect of the Sun is strong enough to cause such a perturbation wrt the keplerian unperturbed orbit. I performed the integration in a SIMULINK environment, using a variable step RK4(5) method with a maximum step size of 10 seconds and a relative tolerance of 1e-13, so my hope is that if there is any problem it will not be related to the solver (which means that some how I can lukily manage it ). What do you think about it? Looking forward to hear from you, I thank you for your attention. Regards

DH, many thanks for the great explanation. You made it very clear and linear. I went through the mathematics and then I implemented this procedure in my SIMULINK model and it works fine. I got values between 2.96e-6 and 2.54-6 for the Z component and between 10e-9 and 10e-11 for the X component, do you think they are consistent? I made a check with the orbital period and I got 27.37 days calculated with the average value of the angular velocity over one orbit around the Earth. Looking forward to hear from you, I thank you for your kind help. Regards

Hi D H, many thanks for the hints and the clarification concerning the quaternion time derivative. At the moment, I evaluate the angular velocity of the rotating frame as omega = cross(r, v) / (norm®^2), where r and v are the position and the velocity of the Moon in ECI and the angular acceleration using the previous step value and the time step. Doing like this I should find the angular velocity and acceleration of the rotating frame in the inertial frame, am I right? I hope you have some suggestion about this. Cheers

DH, I spent the last days trying to manage with the problem of the determination of the Moon's angular velocity vector. The path I followed uses quaternions as you suggested. What I did is: 1. propagate the Moon's position in the ECI frame 2. compute the "attitude" of the Moon from the (X,Y,Z) coordinates in ECI using Euler angles with sequence 3-1-3 (I know, Euler angles are evil but I couldn't imagine how to compute quaternions from rectangular coordinates, so the intermediate step through the Euler angles and the rotation matrix between ECI and the rotating frame was necessary) 3. convert the rotation matrix to quaternion 4. evaluate the time derivative of the quaternion using the difference quotients definition f(a)_dot = lim(h->0) (f(a+h)-f(a))/h, because the time step h for the integration is very small 5. use the relation between the quaternion time derivative and the angular rate vector. Unfortunately, I got some strange results that make me think that something is wrong in the procedure or in the implementation, because the angular rate in quaternion notation has non zero scalar component, while I expected a vector in quaternion notation with only the "imaginary" parts non zero. Now I have some doubt concerning the evaluation of the quaternions time derivative: I considered the time step as a scalar, but I wonder if I have to transform the time step in quaternion notation and use the quaternion division. Moreover, do you think this procedure is correct? Looking forward to hear from you, I thank you for your help.

Hi everybody! I would like to ask you an opinion. I have a file with time variable quaternions and I would like to numerically evaluate the time derivative, that is given q(i), q(i-1) t(i) and t(i-1) I express the quaternion time derivative as q_dot = (q(i) - q(i-1)) / (t(i) - t(i-1)) where q denotes a quaternion and t is the time. Now my question is: do I have to express the time in a quaternion form [t(i) 0 0 0] and use the quaternion division or can I divide a quaternion by a scalar? I know there is a formula linking the angular velocity to the time derivative of the quaternion, but I don't know the angular rate which has to be determined using this formula. I hope some of you can help me or suggest me another way to cope with this problem. Many thanks for your help. Phate

Arch2008, I saw the link and found a good paper about the CRTBP, with just some hints about the perturbations acting on the system. I think this will be another good starting point for a reflection. Thanks

Ok, I will try to develop an expression of the time derivative of the transformation matrix which includes the third-body effects. This seems to be quite challenging at first sight, can you suggest me a paper so that I can learn more about this topic and start thinking of the problem? Regards Hi, I spent the morning thinking about the problem and finally I pointed out this solution. At each time-step the position of the rotating frame in the inertial one is known (bacause the motion of the Moon is propagated and the rotating frame is always towards the Moon), therefore the attitude of the rotating frame with respect to the inertial frame is known as well, for example via the Euler angles. Using the frame's attitude at the previous time-step it is possible to compute the time derivatives of the Euler angles and from these values the components of the angular velocity vector of the rotating frame (taking into consideration that the perturbation due to the Sun on the Moon is very slow). Then, the vector obtained is the angular velocity vector expressed in the rotating frame at that specific time. What do you think about this approach? I know it is quite a rude and not elegant solution, but do you think it would be an acceptable approximation, taking into account that the simulation has navigation purpose? Best regards

Hi DH, thanks for your reply. I managed the elliptical three body problem in several cases, but I have always studied the problem in the synodic rotating reference frame. Now I am trying to build a more detailed simulator using the ephemerides of the Sun and of the Moon and I perform the integration in the ECI frame. This is why I faced with the transformation between rotating and inertial reference frame for the first time. I am also experienced with the evaluation of the WSB, but all these studies have been made in the rotating frame and with a planar geometry, that is the plane of motion of the Moon around the Earth. Now my wish is to increase the accuracy of my simulations in order to have a better model of the "real world" and evaluate the efficiency of my trajectories taking into consideration also the smallest effects and the perturbations due to the third bodies and the 3D geometry. Besides, what I am doing now lets me really understand the difference between the academic world and the more detailed applications used in the research field. I know that there's a hard work to be done, but I am curious... That said, do you think I may study the perturbations due to the Sun and Jupiter on the Moon's orbit with the planetary equations (and add a time dependent variation of inclination and right ascension) ? At the moment this is the only way I can see, even if I have to think deeply about the question, but I think this depends basically on my knoledge. Do you have any suggestion to compute the X and Y components of the angular velocity vector? Looking forward to hear from you, I thank you for your patience. Regards

So the point is that on one hand there is how the same vector is seen by a fixed and a rotating observer and on the other hand there is how the vector is described in one frame or another. Did I get the difference? I agree with the definition of the z component of the angular velocity (which is responsable of the motion of the Moon on its orbital plane) but my doubt is about the other components of the angular velocity vector. The Sun ( or Jupiter to a much smaller extent) acts has a perturbing "third" body on the Moon's orbit around the Earth, so this should imply a variation of the orbital plane inclination that is related to the other components of the angular velocity vector, am I right? Do you think that these effects linked to the other X and Y components of the angular velocity vector can be neglected? Otherwise I have to admit that the motion of the Moon is as beautifull as complicated... To face the modelling of the real earth - moon system is a very challenging task, I am really surprised of the great amount of details you have to take into consideration. That's why I love mission analysis! Many thanks for your opinions and suggestions I really enjoy this debate, are you a professor? Cheers

D H, many thanks for the prompt reply. I just have a question concerning the basic equation of kinematics. If this equation states a relation between the rotating frame and the inertial one, why it is necessary to trasform the vector via the transformation matrix? Secondly, what about the angular velocity of the rotating frame? You said that it has to be expressed in the rotating frame, that is in the direction of the actual orbital moment, but with which module? The motion is highly non keplerian, so the common rules of the two bodies to obtain the orbital elements (such has the true anomaly) are no more valid and I don't know how to evaluate this quantity. Last, what vectors are Xs and Xr? Again, many thanks for your help. P.S. I apologize for the choice of the topic, but I read "orbits" and I thought it was a generic one about orbit mechanics.

DH, thanks for the replay, I greatly appreciated your post. So you suggest me to compute at each step of the integration: 1. position and velocity of the Moon with respect to the Earth 2. Evaluate the Moon's orbital momentum 3. Define the three unity vectors of the rotating frame 4. Build the rotation matrix having as columns the unity vectors of the rotating frame I just have one question: in this way I should be able to transform positions from one reference frame to another, but what about velocities? If I simply rotate the frame I'd neglect the velocity of the Moon about the Earth, am I wrong? About the position of the libration points, I just saw on Farquhar that for L1 and L2 two parameters exist so that the ratio between the actual Earth - Moon distance and the actual Moon - Libr. Point is constant. Do you referer to this parameters? Many thanks for your help, I hope we will have the possibility to keep on talking of the problem. Regards