Three-body problem reference frames

[phate]

Hi everybody! I am new of the forum and I would like to ask you an opinion. I am setting a high precision simulator for the Earth - Moon system and I am integrating the equations of motion in the ECI reference frame (the earth centered inertial). I also have a spacecraft moving in the vicinity of one collinear lagrangian point and I know its state vector (position and velocity) in the synodic reference frame, that is a reference frame having origin in the earth-moon center of mass, x axis in the direction from the earth to the moon (hence, rotating in the eci frame) and z axis normal to the moon's orbital plane a y axis to form a right hand frame. Now I have to open points:

1. I know that the lagrangian points are solutions of the circular restricted three body problem so they actually do not exist in a more complex system, but is it possible to define "artificial" libration points for a elliptic restricted five body problem?

2. How can I transform the initial conditions for the integration of the spacecraft in the eci reference frame where the equations of motion are written?

I hope that someone of you can somehow help me or just have an opinion exchange.

Best regards and thanks for reading.

[Klaynos]

1. It depends on the 5 body problem, it might be a bit difficult to work out though, you'll need to sum teh vector fields, find the magnitude and see if it's possible to have a 0 magnitude at any point.

 

2. You'll need to know the transform from one frame to another, it shouldn't be too hard to work out in classical mechanics just work out how one frame is moving compared to the other.

[phate]

Klaynos thanks for the prompt reply. Concerning the definition of the lagrangian points I agree with you, that was the way I figured out but it seems to be an extremely hard work to cope with. For the second point, my doubt was about the transformation matrix. I mean, if the problem would have been planar the solution would have been a simple rotation about the Z axis of the inertial reference frame for the position and the composition of the velocities according to the relative motion theorem. But in my case, I use an ephemeris model so the problem is 3D and I find hard to define the angular rate vector of the Moon at the epoch in the ECI frame. It would be nice to be able to compute the rotation matrix starting from the Inertial coordinates of the Moon...any suggestion?

Thanks

 

Hi all!

After some time thinking about the problem, I came to this conclusion about the coordinate transformation: as far as the synodic reference frame is fixed with respect to the Moon (meaning that it is rotating with the Moon), position and velocity of the spacecraft in the rotating reference frame can be evaluated by subtracting the inertial state vector of the spacecraft and the inertial state vector of the Moon so that the velocity of the spacecraft in the inertial frame is the velocity of the moving frame (the Moon) plus the motion of the spacecraft with respect to the moving frame. Do you think I am right?

Regards

[Klaynos]

It does seem like alot of hard work, but I fear it might be the only method

 

And your solution to your frame situation seems to sound about right, although you must remember that orbiting frames arn't really inertial frames as they're accelerating But in most cases this can be pretty much ignored...

[D H]

phate,

 

I'll discuss the elliptical restricted three-body problem and then discuss generalizations to the real-world situation, where you obtain planetary positions from some ephemeris model.

 

The first generalization from the circular restricted three-body problem is the elliptical restricted problem. The location of the colinear libration points in the rotating frame obviously varies with time in the case of elliptical orbits. However, if you scale the reference frame axes with a time-varying scale factor so that the Earth-Moon distance is always one, the co-linear libration points once again take on constant locations.

 

The rotating frame in the elliptical is defined by

• [math]\mathbf r_{E\to M}[/math] Vector from the Earth to the Moon.
  
• [math]\mathbf v_{E\to M}[/math] Time derivative of [math]\mathbf r_{E\to M}[/math], non-rotating observer.
  
• [math]\mathbf l_{E\to M}[/math] Specific angular momentum of the Moon with respect to the Earth: [math]\mathbf l_{E\to M} = \mathbf r_{E\to M}\times\mathbf v_{E\to M}[/math]
  
• [math]\hat x[/math] Rotating frame x-axis unit vector: [math]\hat x = \mathbf r_{E\to M}/||\mathbf r_{E\to M}||[/math]
  
• [math]\hat z[/math] Rotating frame z-axis unit vector: [math]\hat z = \mathbf l_{E\to M}/||\mathbf l_{E\to M}||[/math]
  
• [math]\hat y[/math] Rotating frame y-axis unit vector: [math]\hat y = \hat z \times \hat x[/math]
  
• [math]\mathbf T_{I\to R}[/math] ECI-to-rotating frame transformation matrix: [math]\mathbf T_{\text{ECI}\to\text{rot}} = \bmatrix \hat x & \hat y & \hat z\endbmatrix[/math]

 

This will give a time-varying transformation matrix from ECI to the rotating frame. Not that since this transformation matrix is a function of the Moon's relative velocity, the time derivative of this transformation involves the Moon's acceleration. Combining with the scaling of the the axes by the inverse of the Earth-Moon separation yields

 

You can bring this concept of time-varying scaling of the reference frame axes forward to a real-world situation where the Earth and Moon no longer undergo even elliptical orbits about the Earth-Moon barycenter.

 

My wife is meeting me for lunch, so I have to run off now. I'll try to add more later.

[phate]

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

[D H]

phate,

 

First, you are dealing with a rotating frame. Velocities (or the time derivative of any vector quantity) do not simply transform from one frame to another. The time derivatives of some vector quantity [math]\mathbf q[/math] from the perspective of observers fixed with respect to the inertial and rotating frames are related via

 

[math]{\frac{d\mathbf q}{dt}}_{\text{inertial}} =

{\frac{d\mathbf q}{dt}}_{\text{rotating}} + \boldsymbol \omega \times \mathbf q[/math]

 

Explicitly incorporating the reference frames,

 

[math]\dot{\mathbf q}_I =

\mathbf T_{R \to I}\left(

\dot{\mathbf q}_R +

\boldsymbol \omega_R \times \mathbf q_R\right)

[/math]

 

where [math]\mathbf T_{R \to I}[/math] is the rotating-to-inertial transformation matrix and [math]\boldsymbol \omega_R[/math] is the angular velocity of the rotating frame with respect to the inertial frame expressed in rotating frame coordinates.

 

Second, if you do as I suggested and scale the spatial dimensions so that the distance between the Earth and Moon is always unity (this is a very common treatment of the elliptical and perturbed elliptical three-body body), you are dealing with a time-varying scaled reference frame. With this,

 

[math]\dot{\mathbf x}_I =

\mathbf T_{R \to I}\left(

\dot r_{E \to M} \mathbf x_S +

r_{E \to M}\left(\dot{\mathbf x}_S +

\boldsymbol \omega_R \times \mathbf x_S\right)\right)

[/math]

 

where [math]r_{E \to M}[/math] is the magnitude of the unscaled Earth-to-Moon distance, [math]\mathbf x_I[/math] is the position of some object in the unscaled inertial frame, and [math]\mathbf x_S[/math] is the position of the object in the scaled rotating frame.

 

Note

I moved these posts from stevo247's orbits thread to a separate thread.

 

To phlate: In the future, if you want to start a discussion on something unrelated to the topic of discussion in an existing thread, please start a new thread. It is not a good idea to have multiple threads of discussion in one thread.

Edited July 22, 2008 by D H

[phate]

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.

Edited July 22, 2008 by phate

[D H]

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?

If you want to deal with multiple observers, you have to be aware of the fact that the time derivatives of some vector quantity as seen from observers fixed with respect to two different reference systems are different vectors. The expression

 

[math]{\frac{d\mathbf q}{dt}}_{\text{inertial}} =

{\frac{d\mathbf q}{dt}}_{\text{rotating}} + \boldsymbol \omega \times \mathbf q[/math]

 

describes the as-observed in relationship between the two vectors [math]{\frac{d\mathbf q}{dt}}_{\text{inertial}}[/math] and [math]{\frac{d\mathbf q}{dt}}_{\text{rotating}}[/math]. It does not say anything about the as-represented in nature of the vectors. As-observed-in and as-represented-in are distinctly different concepts. If you want to deal with multiple observers, you must account for the as-observed-in relationship between the time derivatives of vector quantities. If you want to deal with the vectors numerically (which is exactly what you want to do) you must account for the as-represented-in relationship between the numerical representations of vectors.

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.

The specific angular momentum is easy: [math]\mathbf h_{E\to M} = \mathbf r_{E\to M} \times \mathbf v_{E\to M}[/math]. The angular velocity is a lot tougher. The (rotating frame) z-component of the angular velocity vector is just [math]h_{E\to M} / r_{E\to M}^{\;2}[/math]. The existence of other bodies (e.g. the Sun, and Jupiter to a much lesser extent) essentially act as a torque on the Earth-Moon rotating frame. One result of this torque is to add a radial component to the angular velocity vector.

 

Last, what vectors are Xs and Xr?

• [math]\mathbf x_I[/math] is the position of some object in non-rotating (e.g., ECI) coordinates.
  
• [math]\mathbf x_R[/math] is the position of the object in the Earth-Moon rotating frame: [math]\mathbf x_I = \mathbf T_{R\to I} \mathbf x_R[/math]
  
• [math]\mathbf x_S[/math] is the scaled position of the object in the Earth-Moon rotating frame: [math]\mathbf x_R = r_{E\to M} \mathbf x_S[/math]

[phate]

If you want to deal with multiple observers, you must account for the as-observed-in relationship between the time derivatives of vector quantities. If you want to deal with the vectors numerically (which is exactly what you want to do) you must account for the as-represented-in relationship between the numerical representations of vectors.

 

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?

 

The specific angular momentum is easy: [math]\mathbf h_{E\to M} = \mathbf r_{E\to M} \times \mathbf v_{E\to M}[/math]. The angular velocity is a lot tougher. The (rotating frame) z-component of the angular velocity vector is just [math]h_{E\to M} / r_{E\to M}^{\;2}[/math]. The existence of other bodies (e.g. the Sun, and Jupiter to a much lesser extent) essentially act as a torque on the Earth-Moon rotating frame. One result of this torque is to add a radial component to the angular velocity vector.

 

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

Edited July 22, 2008 by phate

[D H]

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?

This variation in the inclination (and right ascension; two angles are needed to describe the orientation of a plane in three-space) is exactly what I was talking about. The Earth and Moon are orbiting their common center of mass. If those third bodies were not present the Earth-Moon orbital plane would be fixed. The orbital plane is not fixed because of those third body effects. The in-plane component of the angular velocity encapsulates those effects.

 

Do you think that these effects linked to the other X and Y components of the angular velocity vector can be neglected?

To first order, yes. You can think of the circular restricted three-body problem as a zeroth-order description of behavior some object in the vicinity of the L1 or L2 points. The next step up in complexity is the elliptical restricted three-body problem. This captures a lot of the details missed by the circular problem. The last steps to make are to consider the Sun, and then the full solar system. I suggest you not to make these steps until you have worked on the elliptical problem.

 

Many thanks for your opinions and suggestions I really enjoy this debate, are you a professor?

No. I work at NASA. Modeling dynamics is my bread-and-butter.

[phate]

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

[D H]

The angular velocity depends on the gravitational acceleration of the Earth and Moon toward each other and toward all other bodies in the solar system, which shouldn't be too surprising since the angular velocity of the isolated Earth-Moon system is a function of the gravitational acceleration toward each other. I assume you are developing your equations of motion in the rotating frame, which means you have an omega-dot term in the equations of motion. This means your equations of motion will depend on gravitational jerk.

 

So, how to go about addressing this? You first need to develop an expression for [math]\boldsymbol \omega[/math]. You will need to manipulate the relation between the time derivative of a transformation matrix and the angular velocity vector:

 

[math]\dot{\mathbf T}_{I\to R}

= -\mathbf{Sk}(\boldsymbol \omega_R) \mathbf T_{I\to R}

= -\mathbf T_{I\to R}\mathbf{Sk}(\boldsymbol \omega_I)

[/math]

 

where [math]\mathbf{Sk}(\mathbf a)[/math] is the skew-symmetric cross product matrix generated from the vector [math]\mathbf a[/math].

 

Outline of approach:

• Generate an expression for [math]\dot{\mathbf T}_{I\to R}[/math] that incorporates third-body effects.
  
• Solve for [math]\mathbf{Sk}(\boldsymbol \omega_R)[/math] via [math]
  \mathbf{Sk}(\boldsymbol \omega_R) = - \dot{\mathbf T}_{I\to R}\mathbf T_{I\to R}^{\;T}[/math]
  
• This should be a skew-symmetric matrix (or you did something wrong).
  
• Extract [math]\boldsymbol \omega_R[/math] from this matrix.
  
• You should find that the y-component (rotating frame) of this vector is identically zero (or you did something wrong).
  
• Now that you have [math]\boldsymbol \omega_R[/math] you can differentiate to determine the angular acceleration.

[phate]

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

Edited July 23, 2008 by phate
multiple post merged

[D H]

First off, I never use Euler angles. Euler angles are evil. Quaternions, transformation matrices, Rodrigues parameters, and modified Rodrigues parameters are much better representation schemes for representing orientations. Modern navigation systems almost exclusively use quaternions or modified Rodrigues parameters. Modified Rodrigues parameters have three parameters (exactly the number needed to specify orientation in three-space) but still suffer from a singularity (at integral multiples of 360 degrees; plain-old-vanilla Rodrigues parameters have singularities at integral multiples of 180 degrees). Quaternions have no singularities but have one more parameter than is needed (and hence need to be normalized).

 

Can't say much more now, off to work.

[Arch2008]

Would this help?

http://www.space.edu/documents/McLaughlin2.pdf

[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

[phate]

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.

Edited July 30, 2008 by phate

[D H]

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.

This is going to take a while. What you really want are the angular velocity and its time derivative, as these are what you need to model the equations of motion in a rotating frame. You don't need quaternions to do derive these two vectors. In this post I will address your concern about the nature of the quaternion derivative. In the next post, I will help you get a start on deriving the angular velocity and angular acceleration of the rotating frame.

 

Regarding your concern on the time derivative of the quaternion, the relation between the time derivative of a left¹ transformation² unit quaternion and the angular rate is

 

[math]

\dot {\mathcal Q}_{I\to R} =

\bmatrix 0\\-\,\frac 1 2 \boldsymbol{\omega}_{:R}\endbmatrix

\, {\mathcal Q}_{I\to R}

[/math]

 

In short, there is no reason to expect the time derivative of the quaternion to be pure imaginary.

 

 

====================

Footnotes

 

¹ There are two ways to use unit quaternions to represent rotations/transformations in 3-space:

 

[math]

\bmatrix 0 \\ \boldsymbol{x}_{:B}\endbmatrix =

\mathcal Q_{A\to B}\,

\bmatrix 0 \\ \boldsymbol{x}_{:A}\endbmatrix

\mathcal Q_{A\to B}^{\ast}

[/math]

 

and

 

[math]

\bmatrix 0 \\ \boldsymbol{x}_{:B}\endbmatrix =

\mathcal Q_{A\to B}^{\ast}\,

\bmatrix 0 \\ \boldsymbol{x}_{:A}\endbmatrix

\mathcal Q_{A\to B}

[/math]

 

The only difference between the two is whether the quaternion is on the left (and its conjugate on the right) of the vector to be transformed, or vice-versa. I call the first form left quaternions and the latter, right quaternions. Very few (I have yet to see one!) books or papers on using quaternions to represent rotations/transformations in 3-space address this source of confusion. They simply pick one form as if that is the only way to do it. The advantage of left quaternions is (a) that the left quaternion form looks very much like a symmetry transformation, and (b) left transformation quaternions chain right-to-left, exactly how transformation matrices chain: [math]\mathcal Q_{A\to C} = \mathcal Q_{B\to C}\,\mathcal Q_{A\to B}[/math] for left quaternions. The chief advantage of right transformation quaternions is that they do not chain in the goofy way that transformation matrices do: [math]\mathcal Q_{A\to C} = \mathcal Q_{A\to B}\,\mathcal Q_{B\to C}[/math] for right quaternions. The chief advantage in knowing that there are two equally valid ways of doing things is that this enables one to act as an intermediary between two groups who don't know any better and think the other side is "wrong".

 

² I'm also very picky regarding the terms rotation and transformation. By rotation, I mean a physical rotation of some object from one orientation to another. A transformation on the other hand enables one to transform the representation of some thing as expressed in some system to the representation of the exact same thing on some other system. Rotation and transformation are duals of one another.

[phate]

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

[D H]

First off, a correction:

ECI-to-rotating frame transformation matrix: [math]\mathbf T_{\text{ECI}\to\text{rot}} = \bmatrix \hat x & \hat y & \hat z\endbmatrix[/math]

Oops! That is the rotating-to-inertial transformation matrix.

 

I'll start very generically here. Suppose you have at hand the representations in some n-dimensional Cartesian reference frame [math]A[/math] of the n unit vectors of some other Cartesian reference frame [math]B[/math]. Denoting these unit vectors as [math]\lbrace \hat{\mathbf u}_{B_i:A} \rbrace[/math], where the subscript [math]B_i[/math] indicates the i^th unit vector of reference frame [math]B[/math] and the subscript [math]:A[/math] indicates the unit vector is expressed in terms of reference frame [math]A[/math]. The matrix that transforms column vectors from frame A to frame B in terms of these unit vectors is

 

[math]\mathbf T_{A\to B} =

\bmatrix

\hat{\mathbf u}_{B_1:A}^{\;T} \\

\hat{\mathbf u}_{B_2:A}^{\;T} \\

\cdots \\

\hat{\mathbf u}_{B_n:A}^{\;T}

\endbmatrix

[/math]

 

As a sanity check, [math]\mathbf T_{A\to B}\,\hat{\mathbf u}_{B_i:A} = \hat {\mathbf e}_i\quad(\hat {\mathbf e}_{i_j} \equiv \delta_{ij})[/math], as desired.

 

Back to the problem at hand:

 

In [math]\mathcal R^3[/math], the time derivative of some vector [math]\mathbf q[/math] as observed by an observer fixed in reference frame A is related to the derivative as observed by an observer fixed in reference frame B via

 

[math]\frac{d\mathbf q}{dt}\Bigl|_A =

\frac{d\mathbf q}{dt}\Bigl|_B\, +\,\boldsymbol{\omega}_{A\to B}\times \mathbf q[/math]

 

where [math]\boldsymbol{\omega}_{A\to B}[/math] is the rotation rate of frame B with respect to frame A. The subscript [math]A\to B[/math] on [math]\boldsymbol{\omega}[/math] is going to make the math get ugly, so I will simply use [math]\boldsymbol{\omega}[/math] from this point on.

 

The above equation can be expressed in matrix form by means of the skew symmetric cross product matrix. Denoting [math]\boldsymbol{\omega}^{\mathbf X}[/math] as this matrix,

 

[math]\frac{d\mathbf q}{dt}\Bigl|_A =

\frac{d\mathbf q}{dt}\Bigl|_B\, +\,\boldsymbol{\omega}^{\mathbf X}\, \mathbf q[/math]

 

The natural way to express the derivative of some vector as observed by an observer fixed in some reference frame is to express the vector in terms of that reference frame and differentiate the vector element-by-element:

 

[math]\frac{d\mathbf q}{dt}\Bigl|_{A:A} =

\frac{d\mathbf q_{:A}}{dt}[/math]

 

With the vectors expressed in the reference frame of the observer, the above relation between the two time derivatives becomes

 

[math]\frac{d\mathbf q_{:A}}{dt} =

\mathbf T_{B\to A} \left(

\frac{d\mathbf q_{:B}}{dt} +\,\boldsymbol{\omega}_{:B}^{\;\mathbf X}\, \mathbf q_{:B}\right)[/math]

 

Another way to arrive at the relation is to differentiate the transformed vector [math]\mathbf q_{:A} =

\mathbf T_{B\to A}\,\mathbf q_{:B}[/math]

 

[math]\frac{d\mathbf q_{:A}}{dt} =

\mathbf T_{B\to A}\,\frac{d\mathbf q_{:B}}{dt} +

\dot{\mathbf T}_{B\to A}\,\mathbf q_{:B}[/math]

 

Equating the two results,

[math]

\mathbf T_{B\to A} \left(

\frac{d\mathbf q_{:B}}{dt} +\,\boldsymbol{\omega}_{:B}^{\;\mathbf X}\, \mathbf q_{:B}\right) =

\mathbf T_{B\to A}\,\frac{d\mathbf q_{:B}}{dt} +

\dot{\mathbf T}_{B\to A}\,\mathbf q_{:B}[/math]

 

Canceling the common term yields

 

[math]

\mathbf T_{B\to A}\,

\boldsymbol{\omega}_{:B}^{\;\mathbf X}\, \mathbf q_{:B} =

\dot{\mathbf T}_{A\to B}^{\;T}\,\mathbf q_{:B}[/math]

 

For the above to be true for any vector quantity [math]\mathbf q[/math],

 

[math]

\mathbf T_{B\to A}\,

\boldsymbol{\omega}_{:B}^{\;\mathbf X} =

\dot{\mathbf T}_{B\to A}[/math]

 

The above can be used to determine the rotation rate if the transformation matrix and its time derivatives are known:

 

[math]

\boldsymbol{\omega}_{:B}^{\;\mathbf X} =

\mathbf T_{B\to A}^{\;T}\,\dot{\mathbf T}_{B\to A}[/math]

 

Denoting [math]\hat u, \hat v, \hat w[/math] as the unit vectors for frame B expressed in terms of frame A,

 

[math]\boldsymbol{\omega}_{:B}^{\;\mathbf X}

= \bmatrix

\hat u\cdot\dot{\hat u}&\hat u\cdot\dot{\hat v}&\hat u\cdot\dot{\hat w}\\

\hat v\cdot\dot{\hat u}&\hat v\cdot\dot{\hat v}&\hat v\cdot\dot{\hat w}\\

\hat w\cdot\dot{\hat u}&\hat w\cdot\dot{\hat v}&\hat w\cdot\dot{\hat w}\endbmatrix

[/math]

 

This matrix is indeed skew symmetric since [math]\hat u_i\cdot \hat u_j = \delta_{ij}[/math] Differentiating wrt time, [math]\frac{d}{dt}(\hat u_i\cdot \hat u_j = hat u_i \cdot \dot{\hat u}_j + \dot{hat u}_i \cdot \hat u_j = 0[/math] This reduces to zero for [math]j=i[/math]. From the above,

 

[math]

\begin{aligned}

\omega_u &= \hat w\cdot\dot{\hat v} = -\hat v\cdot\dot{\hat w} \\

\omega_v &= \hat u\cdot\dot{\hat w} = -\hat w\cdot\dot{\hat u} \\

\omega_w &= \hat v\cdot\dot{\hat u} = -\hat u\cdot\dot{\hat v}

\end{aligned}

[/math]

 

The frame rotating with the Moon's orbit about the Earth is defined in terms of the displacement vector [math]\mathbf r_{E\to M:I}[/math] from the center of the Earth to the center of the Moon expressed in some inertial system; [math]\dot{\mathbf r}_{E\to M:I}[/math], the inertial time derivative of this vector; and [math]\mathbf h_{E\to M:I} = \mathbf r_{E\to M:I} \times \dot{\mathbf r}_{E\to M:I}[/math], the specific angular momentum of the Moon with respect to the Earth. With these, the rotating frame unit vectors expressed in inertial coordinates are

 

[math]

\begin{aligned}

\hat {\mathbf r} &= \frac{\mathbf r_{E\to M}}{||\mathbf r_{E\to M}||} \\

\hat {\boldsymbol \theta} &= \hat {\mathbf r} \times \hat {\mathbf h} \\

\hat {\mathbf r} &= \frac{\mathbf h_{E\to M}}{||\mathbf h_{E\to M}||}

\end{aligned}

[/math]

 

Applying the general expression for the components of the frame angular velocity vector to this situation,

 

[math]

\begin{aligned}

\omega_r &=

\hat h\cdot\dot{\hat \theta} = -\hat {\theta}\cdot\dot{\hat h} \\

\omega_{\theta} &=

\hat r\cdot\dot{\hat h} = -\hat h\cdot\dot{\hat r} \\

\omega_h &=

\hat {\theta}\cdot\dot{\hat r} = -\hat r\cdot\dot{\hat {\theta}}

\end{aligned}

[/math]

 

Since the Moon's angular momentum vector is normal to the velocity vector, the [math]\hat \theta[/math] component of the frame angular velocity vector is zero. With a little manipulating, the [math]\hat h[/math] component of the frame angular velocity vector is simple [math]h_{E\to M}/r^2_{E\to M}[/math]. The radial component involves the out-of-plane acceleration of the Earth and Moon toward the Sun, Jupiter, ...

[phate]

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