Orbiter-Forum  

Go Back   Orbiter-Forum > Far Side of the Moon > Math & Physics
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

Math & Physics Mathematical and physical problems of space flight and astronomy.

Reply
 
Thread Tools
Old 07-03-2018, 09:20 AM   #1
MontBlanc2012
Orbinaut
Default From BCI to BCBF and back again

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 x-y-z 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 J_2 (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:

Code:
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 $\tau$ and the other angle $\epsilon$, otherwise known as the obliquity. In Orbiter, the obliquity is assumed constant. However, the value of $\tau$ is a function of the Modified Julian Date (MJD) and is given by the explicit expression:

\tau = \text{LAN} + 2\,\pi\,\frac{\delta d}{T_p}

where

\delta d = \text{mjd} - \text{LAN\_MJD}

As you can see, these expressions are a function of 'LAN', 'LAN_MJD', T_p (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. \delta d is the number of days since the reference epoch given by 'LAN_MJD'. 'LAN', 'LAN_MJD' and T_p are all constants defined in the body's Orbiter configuration file. For example, the .cfg for Mars we have:

Code:
PrecessionPeriod = -63346652.48
LAN = 0.6210531483
LAN_MJD = 51544.5
The angle \tau describes how far the rotation axis has precessed about the precession axis since the start of the epoch. The angle \epsilon (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:

Code:
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 \psi. 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 \psi is:

\psi = 2\,\pi\,\frac{86400}{T_s}\,\delta d - 2\,\pi\,\frac{\cos\epsilon}{T_p}\,\delta d + \psi_0

This expression is a function of a two new parameters, T_s and \psi_0. T_s is the sidereal rotation period of the body; and \psi_0 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:

Code:
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:

Code:
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 (LAN_p) and the Precession Obliquity (\epsilon_p), 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 \delta d - i.e, the number of days since the start of epoch. Then we can calculate the angle \tau using the formula given above. The two angles \tau and the obliquity \epsilon are then used to make the transformation.

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

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

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)

and

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)

Then, we apply R_1 and R_2 to the vector X_{BCI} in succession - i.e.,

X_2 = R_2(\epsilon_p).R_1(\text{LAN}_p).X_{BCI}

To apply the second step, we do apply two further rotations, R_3(\tau) and R_4(\epsilon) defined as follows:

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

and

R_4(\epsilon) = \left(
\begin{array}{ccc}
 1 & 0 & 0 \\\\
 0 & \cos \epsilon & -\sin \epsilon \\\\
 0 & \sin \epsilon & \cos \epsilon \\
\end{array}
\right)

Then, we apply R_3 and R_4 to the vector X_2 in succession - i.e.,

X_4 = R_4(\epsilon).R_3(\tau).X_2

Finally, to apply the third step, we now rotate by a further angle \psi as defined above:

R_5(\psi) = \left(
\begin{array}{ccc}
 \cos \psi & \sin \psi & 0 \\\\
 -\sin \psi & \cos \psi & 0 \\\\
 0 & 0 & 1 \\
\end{array}
\right)

such that:

X_{BCBF} = R_5(\psi).X_4

The full BCI to BCBF transformation, then, is given by:

X_{BCBF} = R_5(\psi).R_4(\epsilon).R_3(\tau).R_2(\epsilon_p).R_1(\text{LAN}_p).X_{BCI}

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:

Code:
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:

Code:
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, X 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:

X_{BCI} = R_1(\text{-LAN}_p).R_2(-\epsilon_p).R_3(-\tau).R_4(-\epsilon).R_5(-\psi).X_{BCBF}

So, the equivalent Python code for the inverse transformation from the BCBF reference frame to the BCI reference frame is:

Code:
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:

X=\{4000000.0, 0.0, 0.0\}

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:

Code:
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:

X' = \{561155.82289003, 3535566.12080444, 1784622.18630623\}

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:

X' = \{561155.82289003, 3535566.12080444, 1784622.18630623\}

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:

X'_{Orb} = \{561155.82289003, 1784622.18630623, 3535566.12080444\}

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 is offline   Reply With Quote
Thanked by:
Reply

  Orbiter-Forum > Far Side of the Moon > Math & Physics


Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump


All times are GMT. The time now is 03:25 AM.

Quick Links Need Help?


About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2018, Jelsoft Enterprises Ltd.
Copyright 2007 - 2017, Orbiter-Forum.com. All rights reserved.