From BCI to BCBF and back again

[MontBlanc2012]

This note sets out the transformations that Orbiter uses to convert points in its 'body-centred inertial' (BCI)reference frame to its 'body-centred, body-fixed' reference frame. (The BCI reference frame is a generalisation of the Earth-centred inertial (ECI) reference frame applying to any gravitating body not just the Earth. Likewise, the BCBF reference frame is a generalisation of the Earth-centred Earth-fixed (ECEF) reference frame again applying to any gravitating body.)

These transformations have a number of uses. In subsequent posts I will focus on using them to develop high fidelity numerical integration tools. In turn, I will use these integration tools to solve the targeting (rendezvous) problem. In this note, however, I will just define the transformations from the BCI to BCBF reference frame and its inverse.

You should note that the standard Orbiter 2016 documentation includes a description of these forward and backward transformations. However, I find this description to be a little opaque. Hopefully, this note may serve to augment the Orbiter documentation.

In this note, I will also include some explicit Python code that carries out the forward and backward transformations.


A quick recap on BCI and BCBF reference frames
To briefly summarise, the BCI reference frame is a reference frame whose origin lies at the centre of mass of the gravitating body but does not co-rotate with that gravitating body. The BCI reference frame is the natural reference frame in which to express the equations of motion of a satellite. Because the BCI reference frame is not an inertial reference frame, we need to add non-inertial tidal forces to arising from other gravitating bodies to the equations of motion. For the Earth, it is common to add tidal contributions from the Moon and from the Sun.

The BCBF reference frame is also centred on the gravitating body's centre of mass but its [imath]x[/imath]-[imath]y[/imath]-[imath]z[/imath] axes rotate with the instantaneous rotation axes of the gravitating body. This means that for any point on the surface of the gravitating body has fixed coordinates in the BCBF reference frame. The BCBF reference frame is the natural reference frame in which to express non-spherical gravitational contributions arising from the [imath]J_2[/imath] (and possibly higher order) gravitational field term.


Rotations included in the transformations
In Orbiter, to transform a vector from the BCI reference frame to the BCBF reference frame, we need to model:

1. the rotation of the gravitating body about its spin axis; and
2. the precession of the spin axis.

To do this, we need five angles:

1. two angles to define the body's precession axis - i.e., the axis about which the gravitating body's rotation axis is assumed to precess;
2. two angles to define the orientation of the rotation axis relative to the precession axis on any given date; and
3. one addition angle to define how far the gravitating body has rotated about the rotation axis since some defined starting epoch.

In Orbiter, the two fixed angles used to define the orientation of the precession axis are called "PrecessionLAN" and "PrecessionObliquity" in the body's Orbiter configuration files. For example, from the .cfg file for Mars we have:

PrecessionLAN = 4.005081124
PrecessionObliquity = 0.03224369545

Next we need two Euler angles to define the orientation of the rotation axis relative to the precession axis. In this note, I shall call one angle [imath]\tau[/imath] and the other angle [imath]\epsilon[/imath], otherwise known as the obliquity. In Orbiter, the obliquity is assumed constant. However, the value of [imath]\tau[/imath] is a function of the Modified Julian Date (MJD) and is given by the explicit expression:

[math]\tau = \text{LAN} + 2\,\pi\,\frac{\delta d}{T_p}[/math]
where

[math]\delta d = \text{mjd} - \text{LAN\_MJD}[/math]
As you can see, these expressions are a function of 'LAN', 'LAN_MJD', [imath]T_p[/imath] (the precession period of the rotation axis about the precession axis) and 'mjd'. Now, 'mjd' is simply the Modified Julian Date for which you wish to calculate the BCI to BCBF transformation. [imath]\delta d[/imath] is the number of days since the reference epoch given by 'LAN_MJD'. 'LAN', 'LAN_MJD' and [imath]T_p[/imath] are all constants defined in the body's Orbiter configuration file. For example, the .cfg for Mars we have:

PrecessionPeriod = -63346652.48
LAN = 0.6210531483
LAN_MJD = 51544.5

The angle [imath]\tau[/imath] describes how far the rotation axis has precessed about the precession axis since the start of the epoch. The angle [imath]\epsilon[/imath] (or the obliquity) defines the 'tilt' of the body's rotation axis relative to the body's precession axis. In Orbiter, this tilt is treated as a constant and is defined in the body's configuration file. For example, the .cfg file for Mars we have:

Obliquity = 0.4397415938

This accounts for four of the five angles. There is one final angle that we need to calculate. In this note, I shall call this angle [imath]\psi[/imath]. It is the angle of rotation that the body has rotated about its instantaneous rotation angle since the start of the epoch given by 'LAN_MJD'. The explicit formula for [imath]\psi[/imath] is:

[math]\psi = 2\,\pi\,\frac{86400}{T_s}\,\delta d - 2\,\pi\,\frac{\cos\epsilon}{T_p}\,\delta d + \psi_0[/math]
This expression is a function of a two new parameters, [imath]T_s[/imath] and [imath]\psi_0[/imath]. [imath]T_s[/imath] is the sidereal rotation period of the body; and [imath]\psi_0[/imath] is the angle of rotation at the start of the epoch given by 'LAN_MJD' and in Orbiter known as the Sidereal rotation offset. Again, these two values are treated as constants in Orbiter and are defined in the body's configuration file. For example, the .cfg file for Mars has:

SidRotOffset = 5.469523488
SidRotPeriod = 88642.66435

In total, eight parameters are needed to define a body's rotation in Orbiter. Collecting them all together, we have for Mars, for example:

PrecessionLAN = 4.005081124
PrecessionObliquity = 0.03224369545
PrecessionPeriod = -63346652.48
LAN = 0.6210531483
LAN_MJD = 51544.5
Obliquity = 0.4397415938
SidRotOffset = 5.469523488
SidRotPeriod = 88642.66435

These eight values are used to define five rotation angles. From these five rotation angles, we define the BCI to BCBF transformation and its inverse.


The BCI to BCBF transformation
The general strategy for transforming a vector from the BCI reference frame to the BCBF reference frame is a three step process:

* step 1: rotate the vector into a coordinate system that is aligned with the precession axis of the body. Two angles, the Precession LAN ([imath]LAN_p[/imath]) and the Precession Obliquity ([imath]\epsilon_p[/imath]), are needed to do this.

* step 2: rotate from a coordinate system aligned with the body's precession axis to a coordinate system aligned with the instantaneous rotation axis of the body. To do this we need to calculate the value [imath]\delta d[/imath] - i.e, the number of days since the start of epoch. Then we can calculate the angle [imath]\tau[/imath] using the formula given above. The two angles [imath]\tau[/imath] and the obliquity [imath]\epsilon[/imath] are then used to make the transformation.

* step 3: rotate the new vector about the rotation axis by an amount [imath]\psi[/imath]. The result is the transformed vector in the BCBF reference frame.

Let's suppose that a vector [imath]X_{BCI}[/imath] is given to us in BCI coordinates. To carry out step 1, we need to define two rotation matrices:

[math]R_1(\text{LAN}_p) = \left( \begin{array}{ccc} \cos \text{LAN}_p & \sin \text{LAN}_p & 0 \\\\ -\sin \text{LAN}_p & \cos \text{LAN}_p & 0 \\\\ 0 & 0 & 1 \\ \end{array} \right)[/math]
and

[math]R_2(\epsilon_p) = \left( \begin{array}{ccc} 1 & 0 & 0 \\\\ 0 & \cos \epsilon_p & -\sin \epsilon_p \\\\ 0 & \sin \epsilon_p & \cos \epsilon_p \\ \end{array} \right)[/math]
Then, we apply [imath]R_1[/imath] and [imath]R_2[/imath] to the vector [imath]X_{BCI}[/imath] in succession - i.e.,

[math]X_2 = R_2(\epsilon_p).R_1(\text{LAN}[I]p).X[/I]{BCI}[/math]
To apply the second step, we do apply two further rotations, [imath]R_3(\tau)[/imath] and [imath]R_4(\epsilon)[/imath] defined as follows:

[math]R_3(\tau) = \left( \begin{array}{ccc} \cos \tau & \sin \tau & 0 \\\\ -\sin \tau & \cos \tau & 0 \\\\ 0 & 0 & 1 \\ \end{array} \right)[/math]
and

[math]R_4(\epsilon) = \left( \begin{array}{ccc} 1 & 0 & 0 \\\\ 0 & \cos \epsilon & -\sin \epsilon \\\\ 0 & \sin \epsilon & \cos \epsilon \\ \end{array} \right)[/math]
Then, we apply [imath]R_3[/imath] and [imath]R_4[/imath] to the vector [imath]X_2[/imath] in succession - i.e.,

[math]X_4 = R_4(\epsilon).R_3(\tau).X_2[/math]
Finally, to apply the third step, we now rotate by a further angle [imath]\psi[/imath] as defined above:

[math]R_5(\psi) = \left( \begin{array}{ccc} \cos \psi & \sin \psi & 0 \\\\ -\sin \psi & \cos \psi & 0 \\\\ 0 & 0 & 1 \\ \end{array} \right)[/math]
such that:

[math]X_{BCBF} = R_5(\psi).X_4[/math]
The full BCI to BCBF transformation, then, is given by:

[math]X_{BCBF} = R_5(\psi).R_4(\epsilon).R_3(\tau).R_2(\epsilon_p).R_1(\text{LAN}[I]p).X[/I]{BCI}[/math]
Now, this transformation looks fairly formidable. But actually it can be encoded quite simply. In Python, for example, let's fix the eight required parameters of the transformation using the Mars .cfg file values:

pLan = 4.005081124    # Precession LAN
pObl = 0.03224369545  # Precession obliquity
Tp   = -63346652.48   # Precession period
LAN  = 0.6210531483   # LAN
mjdL = 51544.5        # LAN MJD
Ts   = 88642.66435    # Sidereal rotation period
eps  = 0.4397415938   # Obliquity
psi0 = 5.469523488    # Sidereal rotation offset

Having done this, we can define a BCI to BCBF transformation function as follows:

import numpy as np

def BCI2BCBF( X, mjd ):
 
    deld = mjd - mjdL
    tau  = LAN + 2 * np.pi * deld / Tp
    psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

    x    = np.empty(6)
    y    = np.empty(6)
    z    = np.empty(6)
 
    x[0] = X[0]
    y[0] = X[1]
    z[0] = X[2]
 
    c    = np.cos(pLan)
    s    = np.sin(pLan)
    x[1] = +x[0] * c + y[0] * s
    y[1] = -x[0] * s + y[0] * c
    z[1] = +z[0]

    c    = np.cos(pObl)
    s    = np.sin(pObl)
    x[2] = +x[1]
    y[2] = +y[1] * c - z[1] * s
    z[2] = +y[1] * s + z[1] * c
 
    c    = np.cos(tau)
    s    = np.sin(tau)
    x[3] = +x[2] * c + y[2] * s
    y[3] = -x[2] * s + y[2] * c
    z[3] = +z[2]
 
    c    = np.cos(eps)
    s    = np.sin(eps)
    x[4] = +x[3]
    y[4] = +y[3] * c - z[3] * s
    z[4] = +y[3] * s + z[3] * c
 
    c    = np.cos(psi)
    s    = np.sin(psi)
    x[5] = +x[4] * c + y[4] * s
    y[5] = -x[4] * s + y[4] * c
    z[5] = +z[4]

    return np.array([x[5],y[5],z[5]])

Given a vector, [imath]X[/imath] in BCI coordinates and a given Modified Julian Date ('mjd'), this routine applies the five rotations in sequence and returns the vector in the BCBF reference frame.


The BCBF to BCI transformation
Since the BCBF to BCI transformation is the inverse of the BCI to BCBF transformation, and since the transformation is simply defined as a sequence of elemental rotations, we can generate the inverse of the transformation by simply reversing the order in which we apply the elemental rotation operations and by changing the signs of all of the rotation angles.

In other words, we have:

[math]X_{BCI} = R_1(\text{-LAN}[I]p).R_2(-\epsilon_p).R_3(-\tau).R_4(-\epsilon).R_5(-\psi).X[/I]{BCBF}[/math]
So, the equivalent Python code for the inverse transformation from the BCBF reference frame to the BCI reference frame is:

import numpy as np

def BCBF2BCI( X, mjd ):
 
    deld = mjd - mjdL
    tau  = LAN + 2 * np.pi * deld / Tp

    psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

    x    = np.empty(6)
    y    = np.empty(6)
    z    = np.empty(6)
 
    x[5] = X[0]
    y[5] = X[1]
    z[5] = X[2]

    c    = np.cos(psi)
    s    = np.sin(psi)
    x[4] = +x[5] * c - y[5] * s
    y[4] = +x[5] * s + y[5] * c
    z[4] = +z[5]

    c    = np.cos(eps)
    s    = np.sin(eps)
    x[3] = +x[4]
    y[3] = +y[4] * c + z[5] * s
    z[3] = -y[4] * s + z[5] * c

    c    = np.cos(tau)
    s    = np.sin(tau)
    x[2] = +x[3] * c - y[3] * s
    y[2] = +x[3] * s + y[3] * c
    z[2] = +z[3]

    c    = np.cos(pObl)
    s    = np.sin(pObl)
    x[1] = +x[2]
    y[1] = +y[2] * c + z[2] * s
    z[1] = -y[2] * s + z[2] * c

    c    = np.cos(pLan)
    s    = np.sin(pLan)
    x[0] = +x[1] * c - y[1] * s
    y[0] = +x[1] * s + y[1] * c
    z[0] = +z[1]

    return np.array([x[0],y[0],z[0]])

An example
As a quick example, let's apply these two transformation routines. Let's suppose that we are considering Mars-centred reference frames. And let's suppose that we start with a vector in the Mars-centred inertial reference frame:

[MATH]X=\{4000000.0, 0.0, 0.0\}[/MATH]
Although the transformation is not restricted to position vectors, we can think of this as a point in space 4000 km from the centre of Mars. We can map this point to the Mars-centred Mars-fixed reference frame by applying the Python transformation for MJD = 52644.5:

X = np.array([4000000.0, 0.0, 0.0])
Y = BCBF2BCI(X, 52644.5)

print(Y)

If we apply the transformation, we get the vector:

[MATH]X' = \{561155.82289003, 3535566.12080444, 1784622.18630623\}[/MATH]
These transformation rules are quite general in that they should apply to any body in Orbiter's Solar System since the planetary rotation and precession is consistently defined by the same body-dependent parameters. Simply by extracting the relevant eight constant parameters from the body configuration file, one can readily obtain the transformation rule appropriate for that body.

Although the transformations have been written in Python (for my convenience!), it should be straightforward to translate these algorithms into C or Lua without much difficulty.


A note on Orbiter's left-handed co-ordinate system
In the above definitions, I have assumed a conventional right-handed coordinate system. Orbiter, however, uses a left-handed coordinate system. In my view, the best way of dealing with Orbiter's curious coordinate convention is to swap the 'y' and 'z' values of any state vector extracted from Orbiter. This converts from a left-handed to a right-handed coordinate system. Then, carry out whatever conventional operations are needed in the right-handed coordinate system. When completed, swap the 'y' and 'z' values again thereby transforming the solution back to Orbiter's left-handed coordinate system. These values can then be exported back into Orbiter.

In the above example, the transformed vector was given by:

[MATH]X' = \{561155.82289003, 3535566.12080444, 1784622.18630623\}[/MATH]
These are the values in the conventional right-handed coordinate system. So, to convert back to Orbiter's left-handed convention, we need to swap the 'y' and 'z; values. Thus, the Orbiter compatible result is the vector:

[MATH]X'_{Orb} = \{561155.82289003, 1784622.18630623, 3535566.12080444\}[/MATH]
These results should agree with Orbiter's estimation of the same transformation nigh-on exactly.


Next steps
In a subsequent post, I want to use these transformations to build a straightforward (but high fidelity) numerical integration scheme.

 

[MontBlanc2012]

This is a quick update on the my preceding posts in this thread.

Of late, I've been going though a process and updating a number of coordinate transformations so that they are based on quaternion multiplications rather than rotation matrices - see, for example, Quaternions, rotations and orbital elements. It occurred to me that the BCI to BCBF transformations (and their inverses) can also be simplified using quaternions. And indeed they can. So, without further ado, and using the notation of posts presented earlier, here is a python code snippet to carry out the transformations using quaternions

import numpy as np
import quaternion

pLan = 4.005081124    # Precession LAN
pObl = 0.03224369545  # Precession obliquity
Tp   = -63346652.48   # Precession period
LAN  = 0.6210531483   # LAN
mjdL = 51544.5        # LAN MJD
Ts   = 88642.66435    # Sidereal rotation period
eps  = 0.4397415938   # Obliquity
psi0 = 5.469523488    # Sidereal rotation offset

def ECI2ECEFNew( X, mjd ):
    
    deld = mjd - mjdL
    tau  = LAN + 2 * np.pi * deld / Tp
    psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

    R1   = np.quaternion(
                +np.cos(pLan/2) * np.cos(pObl/2), 
                -np.cos(pLan/2) * np.sin(pObl/2), 
                -np.sin(pLan/2) * np.sin(pObl/2),
                +np.sin(pLan/2) * np.cos(pObl/2))
    
    R2   = np.quaternion(
                +np.cos((tau + psi)/2) * np.cos(eps/2), 
                -np.cos((tau - psi)/2) * np.sin(eps/2), 
                -np.sin((tau - psi)/2) * np.sin(eps/2), 
                +np.sin((tau + psi)/2) * np.cos(eps/2))

    Q    = R1 * R2
    Y    = np.conj(Q) * np.quaternion(0, X[0], X[1], X[2]) * Q
        
    return np.array([Y.x, Y.y, Y.z])


def ECEF2ECINew( X, mjd ):
    
    deld = mjd - mjdL
    tau  = LAN + 2 * np.pi * deld / Tp
    psi  = 2 * np.pi * 86400.0 / Ts * deld - 2 * np.pi * np.cos(eps) / Tp * deld + psi0

    R1   = np.quaternion(
                +np.cos(pLan/2) * np.cos(pObl/2), 
                -np.cos(pLan/2) * np.sin(pObl/2), 
                -np.sin(pLan/2) * np.sin(pObl/2),
                +np.sin(pLan/2) * np.cos(pObl/2))
    
    R2   = np.quaternion(
                +np.cos((tau + psi)/2) * np.cos(eps/2), 
                -np.cos((tau - psi)/2) * np.sin(eps/2), 
                -np.sin((tau - psi)/2) * np.sin(eps/2), 
                +np.sin((tau + psi)/2) * np.cos(eps/2))

    Q    = R1 * R2
    Y    = Q * np.quaternion(0, X[0], X[1], X[2]) * np.conj(Q)
        
    return np.array([Y.x, Y.y, Y.z])

Enjoy.