API Question Vessel velocity in ECEF?

[GLS]

How can I get the vessel's velocity in the ECEF reference frame?

 

[MontBlanc2012]

How can I get the vessel's velocity in the ECEF reference frame?

I assume that what you want to do here is to transform a velocity vector in the ECI reference frame to the ECEF reference frame.

The general transformation of a position vector, [math]Q_{ECI}(t)[/math], in the ECI reference frame to a position vector in the ECEF reference frame, [math]Q_{ECEF}(t)[/math], is given by a rotation.

[math]Q_{ECEF}(t) = R(t).Q_{ECI}(t)[/math]
The time-dependent rotation matrix, [math]R(t)[/math], for an arbitrary body in Orbiter (and in particular in this case for the Earth), is as set out in From BCI to BCBF and back again.

This is the transformation for position vectors, but we can differentiate the expression to work out the transformation run for the velocity vector. So, differentiating both sides, we get:

[math]V_{ECEF}(t)= Q'_{ECEF}(t) = R'(t).Q_{ECI}(t)+R(t).Q'_{ECI}(t) = R'(t).Q_{ECI}(t)+R(t).V_{ECI}(t)[/math]
In other words, to calculate the velocity vector in the ECEF reference frame, we need to know how both its position and its velocity in the ECI reference frame. Now, we know [math]R(t)[/math] from the above link. The only thing that we don't know is [math]R'(t)[/math]. To work this out, let's look at the expression for [math]R(t)[/math]. From the referenced thread, we have:

[math] R(t) = R_5(\psi).R_4(\epsilon).R_3(\tau).R_2(\epsilon_p).R_1(\text{LAN}_p) [/math]
i.e., for an arbitrary gravitating body in Orbiter it is the product of five separate rotation matrices. However, two of these are time dependent - [math]R_3[/math] and [math]R_5[/math]. The [math]R_3[/math] matrix is just accounts for a very small precession effect of the Earth's axis; while the [math]R_5[/math] matrix accounts for the bulk of the Earth's rotation about its axis via the time-dependent parameter [math]\psi[/math].

If we assume that the earth precesses two slowly to worry about for the purpose of this then we can set its time derivative [math]R'_3(t)[/math] to zero. That means, when calculating the derivate of [math]R(t)[/math], we only need to worry about taking the derivative [math]R'_5(t)[/math]. An explicit expression for the [math]R_5(t)[/math] is given soaking the derivative should be pretty straightforward. in which case you just assemble all of the pieces in the above expressions to give you the result you want with.

[math] R'(t) = \psi'(t)\,\frac{\partial\,R_5(\psi)}{\partial\psi}.R_4(\epsilon).R_3(\tau).R_2(\epsilon_p).R_1(\text{LAN}_p) [/math]
If helpful, I'll write a piece of python code as a template. (Can't do this now - but I can do this a bit later.)

 

[GLS]

I'm not having much success getting this calculation to work... :uhh:
Vessel position in ECEF (I'm 99.9% sure) is like this:

vessel->GetGlobalPos( G_pos );
oapiGlobalToLocal( hEarth, &G_pos, &ECEF_pos );

For the other 3 (ECI_pos, ECI_vel and ECEF_vel) I'm kinda lost between GetGlobalVel(), GetRelativeVel(), GetRelativePos(), and now the math above, which I think I got right.... maybe the problem is that I'm not feeding the correct ECI_***...

 

[Urwumpe]

Remember, Orbiters coordinate systems are left-handed...

 

[GLS]

Remember, Orbiters coordinate systems are left-handed...

Yes, and in the ECEF_pos calculation above, this is missing:

ECEF_pos = _V( ECEF_pos.x, ECEF_pos.z, ECEF_pos.y );

But the issue seems to be me not understanding what frame is the "local frame".

 

[Urwumpe]

Yes, and in the ECEF_pos calculation above, this is missing:

ECEF_pos = _V( ECEF_pos.x, ECEF_pos.z, ECEF_pos.y );

But the issue seems to be me not understanding what frame is the "local frame".

Depends... local to what? Didn't I post you a link to a NASA document with all different coordinate definitions some months ago?

 

[GLS]

Depends... local to what? Didn't I post you a link to a NASA document with all different coordinate definitions some months ago?

The local frame mentioned in Orbiter's API.... in this case it would be local to Earth, but that apparently only refers to the center of the frame, leaving the orientation to make me confused... :facepalm:

 

[Urwumpe]

The local frame mentioned in Orbiter's API.... in this case it would be local to Earth, but that apparently only refers to the center of the frame, leaving the orientation to make me confused... :facepalm:

Noo... Local coordinates are VESSEL-local coordinates. Body coordinates would be a more suitable term.


Global coordinates and relative coordinates share the same axis vectors, but have a different origins.

 

[indy91]

Maybe VESSEL::GetGroundspeedVector using the FRAME_REFLOCAL frame of reference is doing what you want?

 

[GLS]

Maybe VESSEL::GetGroundspeedVector using the FRAME_REFLOCAL frame of reference is doing what you want?

That gives the velocity in relation to the ground below the vessel, but I want in relation to the center of the Earth, in the frame that rotates with Earth.
I could probably use that, plus the latitude and longitude to calculate what I need, but it would only work in a spherical Earth.

---------- Post added at 08:12 PM ---------- Previous post was at 08:11 PM ----------

Noo... Local coordinates are VESSEL-local coordinates. Body coordinates would be a more suitable term.

From OrbiterAPI.h

OAPIFUNC void oapiGlobalToLocal (OBJHANDLE hObj, const VECTOR3 *glob, VECTOR3 *loc);
 * \brief Maps a point from the global frame to a [B]local object frame[/B]......

If hObj is Earth, then it refers to Earth.

Global coordinates and relative coordinates share the same axis vectors, but have a different origins.

That seems to be the case.

---------- Post added 09-10-18 at 12:35 AM ---------- Previous post was 09-09-18 at 08:12 PM ----------

For future reference, compilation of what (I think) I know so far:

//// position in ECEF ////
VECTOR3 G_pos;// "global position", ecliptic, J2000
VECTOR3 ECEF_pos;// position in ECEF
vessel->GetGlobalPos( G_pos );
oapiGlobalToLocal( hEarth, &G_pos, &ECEF_pos );
ECEF_pos = _V( ECEF_pos.x, ECEF_pos.z, ECEF_pos.y );

//// position in ECI ////
VECTOR3 ECEJ_pos;// position in "ECEJ", a.k.a., Earth-centered, ecliptic, J2000
VECTOR3 ECI_pos;// position in ECI
MATRIX3 Ro;// rotation matrix from "Earth-centered"-to-global coordinates
vessel->GetRelativePos( hEarth, ECEJ_pos );
oapiGetPlanetObliquityMatrix( hEarth, &Ro );
ECI_pos = mul( Inverse( Ro ), ECEJ_pos );
ECI_pos = _V( ECI_pos.x, /*not sure if should change this signal*/ECI_pos.z, ECI_pos.y );

 

[martins]

Just a small side note: Since the rotation matrix returned by oapiGetPlanetObliquityMatrix is orthonormal, the inverse is given by the transpose, and Orbiter already provides the function "tmul" to multiply the transpose of a matrix with a vector with the same computational cost as "mul".
So instead of

ECI_pos = mul( Inverse( Ro ), ECEJ_pos );

write

ECI_pos = tmul( Ro, ECEJ_pos );

This should give the same result (check!) and avoids the explicit calculation of the inverse, which is the most expensive operation in this line.

Edit: Incidentally, this is true for all rotation matrices handled by Orbiter. You can always invert the rotation by use of tmul. In the entire Orbiter codebase, I don't think I have a single instance requiring the explicit computation of the inverse. So for squeezing out that last drop of performance from your code, you might want to look for any calls to "Inverse" and see if you can remove them.

BTW, where is "Inverse" implemented? I didn't think I had provided that in the API.

 

[martins]

Back to topic: expressing the velocity in the planet-relative frame (ECEF) requires two steps:

- rotating the velocity vector from the global to planet-local frame (this is the second term in MontBlanc's equation above)
- subtracting the planet's rotation component (this is the first term)

To get the planet's rotation matrix R(t), use oapiGetRotationMatrix. Note that Orbiter's definition of R is the inverse of MontBlanc's definition (rotation from local to global), so you need to use the tmul trick as discussed above for implementing the transformation.

To get the planet's rotation velocity [math]V_p[/math] at the vessel position, you need the magnitude [math]v_P = |V_p|[/math] and direction [math]\hat{V}_P = V_p/v_p[/math].

The magnitude is given by [math]v_P = \omega_P r[/math], where [math]\omega_P = 2\pi/T_P[/math] is the angular velocity (get T_P with oapiGetPlanetPeriod) and r is the distance of the vessel from the rotation axis, which, with your code above, is given by hypot(ECEF_pos.x, ECEF_pos.z).

The direction is given by rotating the the unit vector (0,0,1) from position (1,0,0) (i.e. lat=0, lon=0 on the unit sphere) to the current planet rotation:

[math]\hat{V}_P = \tilde{R}_\mathrm{rot}(t) [0,0,1]^T = [-\sin(s), 0, \cos(s)]^T[/math]
where [math]\tilde{R}_\mathrm{rot}(t)[/math] is given by Eq. 17 of Technotes\precession.pdf, and s is the current planet rotation angle, obtained by oapiGetPlanetCurrentRotation.

So putting everything together,

[math]V_\mathrm{ECEF} = R^T(t)V_\mathrm{glob} - v_P \hat{V}_P[/math]
Let me know if this works.

 

[GLS]

I'll check, thanks!
I was thinking about the "velocity problem" during the day and had an idea, but it probably won't work... :S

---------- Post added at 10:33 PM ---------- Previous post was at 07:41 PM ----------

Back to topic: expressing the velocity in the planet-relative frame (ECEF) requires two steps:

- rotating the velocity vector from the global to planet-local frame (this is the second term in MontBlanc's equation above)
- subtracting the planet's rotation component (this is the first term)

To get the planet's rotation matrix R(t), use oapiGetRotationMatrix. Note that Orbiter's definition of R is the inverse of MontBlanc's definition (rotation from local to global), so you need to use the tmul trick as discussed above for implementing the transformation.

To get the planet's rotation velocity [math]V_p[/math] at the vessel position, you need the magnitude [math]v_P = |V_p|[/math] and direction [math]\hat{V}_P = V_p/v_p[/math].

The magnitude is given by [math]v_P = \omega_P r[/math], where [math]\omega_P = 2\pi/T_P[/math] is the angular velocity (get T_P with oapiGetPlanetPeriod) and r is the distance of the vessel from the rotation axis, which, with your code above, is given by hypot(ECEF_pos.x, ECEF_pos.z).

The direction is given by rotating the the unit vector (0,0,1) from position (1,0,0) (i.e. lat=0, lon=0 on the unit sphere) to the current planet rotation:

[math]\hat{V}_P = \tilde{R}_\mathrm{rot}(t) [0,0,1]^T = [-\sin(s), 0, \cos(s)]^T[/math]
where [math]\tilde{R}_\mathrm{rot}(t)[/math] is given by Eq. 17 of Technotes\precession.pdf, and s is the current planet rotation angle, obtained by oapiGetPlanetCurrentRotation.

So putting everything together,

[math]V_\mathrm{ECEF} = R^T(t)V_\mathrm{glob} - v_P \hat{V}_P[/math]
Let me know if this works.

From the values I getting, both for your solution and mine*, the orbital speed of the Earth is missing.

*) So, I thought about converting the vessel's global velocity into a global position by adding the vessel's global position, and then converting that to ECEF (or ECI) using the logic from my post from yesterday, and finaly subtracting the vessel's ECEF (or ECI) position... but it's not working. :shrug:

 

[martins]

Ah yes, good point. V_glob has to be the relative velocity in the global frame, so

[math] V_\mathrm{glob} = V(\mathrm{vessel})_\mathrm{glob} - V(\mathrm{planet})_\mathrm{glob} [/math]
or simply use oapiGetRelativeVel to obtain the vessel's relative velocity in the planetocentric ecliptic frame.

 

[GLS]

It should be giving 0,0,0 (or something very very small) when landed and it isn't... :uhh:

(It should be weekend again, so I could look at this all day long. :lol

---------- Post added 09-11-18 at 12:32 AM ---------- Previous post was 09-10-18 at 11:11 PM ----------

Tomorrow at lunch I'll "upload" the precession.pdf file to my brain, and maybe tomorrow night I'll have a way to calculate both the ECEF and ECI velocities.

 

[martins]

My bad. In my musings above "s" should not have been the planet rotation angle, but the vessel's longitude in the local planet frame, so s = atan2(ECEF_pos.z, ECEF_pos.x)

Just to make sure I didn't have another bug in that algorithm, I implemented it and it seems to work. Put the following into DeltaGlider::clbkDrawHUD:

	// DEBUG
	VECTOR3 v, r;
	MATRIX3 R;
	OBJHANDLE hVessel = this->GetHandle();
	OBJHANDLE hPlanet = this->GetSurfaceRef();
	oapiGetRelativeVel(hVessel, hPlanet, &v);
	oapiGetRelativePos(hVessel, hPlanet, &r);
	oapiGetRotationMatrix(hPlanet, &R);
	double T = oapiGetPlanetPeriod(hPlanet);
	VECTOR3 vloc = tmul(R, v);
	VECTOR3 rloc = tmul(R, r);
	double s = atan2(rloc.z, rloc.x);
	VECTOR3 vsurf = _V(-sin(s), 0, cos(s)) * PI2 / T*hypot(rloc.x, rloc.z);
	VECTOR3 V = vloc - vsurf;
	char cbuf[256];
	sprintf(cbuf, "V.x=%+0.4lf, V.y=%+0.4lf, V.z=%+0.4lf", V.x, V.y, V.z);
	skp->Text(0, cy + 20, cbuf, strlen(cbuf));

 

[Urwumpe]

It should be giving 0,0,0 (or something very very small) when landed and it isn't... :uhh:

Rotating frame or non-rotating frame? And if rotating frame - proper axis and rotation rate?

 

[GLS]

Rotating frame or non-rotating frame? And if rotating frame - proper axis and rotation rate?

AFAIK, when landed (and stopped) on the surface:

• ECEF velocity = 0
• ECI velocity = rotation of the planet at that latitude (Z component = 0), so 0 at the poles and maximum at the equator
• ECEF position = constant
• ECI position = non-constant (matches ECEF position once a day)

 

[Urwumpe]

AFAIK, when landed (and stopped) on the surface:

• ECEF velocity = 0
• ECI velocity = rotation of the planet at that latitude (Z component = 0), so 0 at the poles and maximum at the equator
• ECEF position = constant
• ECI position = non-constant (matches ECEF position once a day)

Yes, that is correct.

But if you have for example different rotation frames between your model and Orbiter, you should get ECEF velocity >> 0

Especially, if you have for example right-handed vs left-handed coordinate systems - then you might even get ECEF velocity (wrong) = 2 x ECI velocity (correct).

 

[GLS]

Yes, that is correct.

But if you have for example different rotation frames between your model and Orbiter, you should get ECEF velocity >> 0

Especially, if you have for example right-handed vs left-handed coordinate systems - then you might even get ECEF velocity (wrong) = 2 x ECI velocity (correct).

I think "my model", a.k.a. ECEF and ECI, is well defined... :dry:

---------- Post added at 11:53 PM ---------- Previous post was at 11:04 AM ----------

My bad. In my musings above "s" should not have been the planet rotation angle, but the vessel's longitude in the local planet frame, so s = atan2(ECEF_pos.z, ECEF_pos.x)

Just to make sure I didn't have another bug in that algorithm, I implemented it and it seems to work. Put the following into DeltaGlider::clbkDrawHUD:

	// DEBUG
	VECTOR3 v, r;
	MATRIX3 R;
	OBJHANDLE hVessel = this->GetHandle();
	OBJHANDLE hPlanet = this->GetSurfaceRef();
	oapiGetRelativeVel(hVessel, hPlanet, &v);
	oapiGetRelativePos(hVessel, hPlanet, &r);
	oapiGetRotationMatrix(hPlanet, &R);
	double T = oapiGetPlanetPeriod(hPlanet);
	VECTOR3 vloc = tmul(R, v);
	VECTOR3 rloc = tmul(R, r);
	double s = atan2(rloc.z, rloc.x);
	VECTOR3 vsurf = _V(-sin(s), 0, cos(s)) * PI2 / T*hypot(rloc.x, rloc.z);
	VECTOR3 V = vloc - vsurf;
	char cbuf[256];
	sprintf(cbuf, "V.x=%+0.4lf, V.y=%+0.4lf, V.z=%+0.4lf", V.x, V.y, V.z);
	skp->Text(0, cy + 20, cbuf, strlen(cbuf));

Yes, it works!
Plus, ECI velocity:

VECTOR3 v;
vessel->GetRelativeVel( hEarth, v );
MATRIX3 Ro;
oapiGetPlanetObliquityMatrix( hEarth, &Ro );
ECI_vel = tmul( Ro, v );

It took a few days but it's done!
Tomorrow I'll create some functions for this code, put it in libUltra and post it here for future reference.