Transform ECI to ECEF

Question

I need to transform ECI to ECEF coordinate. Here, I need the formula in math and algorithm. I have looking for the formula in the internet, but I can find one. I have x, y, z cartesian coord in ECI. And want to transform it to ECEF coord. Anybody knows the formula and algorithm?

Answer 1

The Explanatory Supplement to the Astronomical Almanac has all the equations you need. Take a look at chapters 3 and 4.

Keep in mind you need some clear definitions of what you mean by ECEF and ECI. Most people utilize WGS84 for ECEF, but that is not a requirement. Similarly ECI could be J2000 or ICRF

In general you will need 4 steps to convert from ECI to ECEF:

1. Calculate earth's precession
2. Calculate earth's nutation
3. Account for earth's rotation including UTC-UT1 offset
4. Account for polar motion
5. Convert to lan/lon if desired

The first two transformations are purely analytical, and can be obtained from books like the Explanatory Supplement. The third and fourth transformation rely on irregularly varying parameters that have to be updated for a specific date.

The help system for the Systems ToolKit (STK) has some excellent descriptions on various frames and transformation process, but not actual equations.

Answer 2

I know this is a old/dead post but I essentially took on the same challenge and when I googled the answer I kept being brought here. I wanted to add some information here so anyone in the future who Google's "How to compute ECI to ECEF or how to compute ECEF to ECI?" Can avoid some of the short falls I fell into or sifting through hundreds of papers on the subject. I winded up buying the very expensive 3rd edition of the "Explanatory Supplement to the Astronomical Almanac" edited by S.E.Urban and P.K.Seidelmann to help me (Great book by the way). Here's my solution. It just does position but I'm planning on understanding velocity and acceleration next.

What Carlos said is dead on. Except what was chapter 3 and 4 in the 2nd edition he linked is now chapter 6. If we define what we mean by ECI as the J2000.0 ... \begin{equation} x_{CRS} = BPNTW\cdot x_{TRS} \end{equation} from page 208 of the book. \begin{equation} B=Bias\space rotation \space matrix\\ P=Precession\space rotation \space matrix\\ N=Nutation\space rotation \space matrix\\ T=Earth\space rotation\space rotation \space matrix\\ W=Polar\space motion\space rotation\space matrix \end{equation} There are multiple ways to compute each and therefor they all change your definition of what ECI means.

Starting with the bias matrix. They give one in the book. There's also others which could be found on the internet like this paper if you're not completely bored this far in.

The next is the precession matrix. This one is a lot more interesting as it's a series of rotations that are dependent on time. It's also where it gets a lot more complicated and you can read about the more complicated form here or in the book. I'm going to describe the much simpler form as a function of zeta, theta, and z (page 217 if you buy the 3rd edition of the book. If you're serious about this I highly recommend it as it's great for all sorts of orbital fun). Precession could be writen as ... \begin{equation} P=R_3(-z_a)R_2(\theta_a)R_3(-\zeta_a)\\ R_3(-z_a) = \begin{pmatrix}cos(z_a) & -sin_(z_a) & 0 \\\ sin_(z_a) & cos(z_a) & 0 \\\ 0 & 0 & 1\end{pmatrix}\\ R_2(\theta_a) = \begin{pmatrix}cos(\theta_a) & 0 & sin_(\theta_a)\\\ 0 & 1 & 0 \\\ -sin_(\theta_a) & 0 & cos(\theta_a)\end{pmatrix}\\ R_3(-\zeta_a) = \begin{pmatrix}cos(\zeta_a) & -sin_(\zeta_a) & 0 \\\ sin_(\zeta_a) & cos(\zeta_a) & 0 \\\ 0 & 0 & 1\end{pmatrix}\\ \end{equation}

Expand it out and you should get the matrix for precession as ... \begin{equation} P=\begin{pmatrix} cos(z_a)cos(\theta_a)cos(\zeta_a)-sin(z_a)sin(\zeta_a) & -cos(z_a)cos(\theta_a)sin(\zeta_a)-sin(z_a)cos(\zeta_a) & -cos(z_a)sin(\theta_a) \\\ sin(z_a)cos(\theta_a)cos(\zeta_a)+cos(z_a)sin(\zeta_a) & -sin(z_a)cos(\theta_a)sin(\zeta_a)+cos(z_a)cos(\zeta_a) & -sin(z_a)sin(\theta_a) \\\ sin(\theta_a)cos(\zeta_a) & -sin(\theta_a)sin(\zeta_a) & cos(\theta_a) \end{pmatrix} \end{equation}

You can look up what each means on your own time as this is only the second matrix and will be a very long explanation if I go into too much detail. They're all functions of time in Julian centuries since J2000.0 which is exactly 36525 days. Here's one set of coefficients there are many others and they will give slightly different answers. They're all in arcseconds per century.

\begin{equation} \zeta_a(T) = 2.650545 + 2306.083227*T + 0.2988499 * T^2 + 0.1801828*T^3 + Higher\space Order\space Terms\\ z_a(T) = -2.650545 + 2306.077181*T + 1.0927348 * T^2 + 0.1826837*T^3 + Higher\space Order\space Terms\\ \theta_a(T) = 2004.191903 * T - 0.4294934 * T^2 - 0.04182264 * T^3 + Higher\space Order\space Terms\\ \end{equation}

Look in my github repo under toluene/config/config.yml for the complete definition of the coefficients.

Now we get into nutation. This one is very complicated and I suggest this reading. Page 45 and 46 in the reading have definitions for psi and epsilon you'll need. Has the added bonus of the nutation series you'll need at the bottom. Page 225 of the book has the nutation rotation matrix from Emerson 1973.

\begin{equation} N= \begin{pmatrix} cos(\Delta\psi) & -sin(\Delta\psi)cos(\epsilon_a) & -sin(\Delta\psi)sin(\epsilon_a) \\\ sin(\Delta\psi)cos(\epsilon) & cos(\Delta\psi)cos(\epsilon)cos(\epsilon_a)+ sin(\epsilon)sin(\epsilon_a) & cos(\Delta\psi)cos(\epsilon)sin(\epsilon_a)-sin(\epsilon)cos(\epsilon_a) \\\ sin(\Delta\psi)sin(\epsilon) & cos(\Delta\psi)sin(\epsilon)cos(\epsilon_a)- cos(\epsilon)sin(\epsilon_a) & cos(\Delta\psi)sin(\epsilon)sin(\epsilon_a) + cos(\epsilon)cos(\epsilon_a) \end{pmatrix} \end{equation}

where

\begin{equation} \epsilon=\epsilon_a+\Delta\epsilon \end{equation}

epsilon_a is a function of time in julian centuries and is given in arc seconds also. The config file I was talking about earlier has the coefficients for epsilon_a also.

\begin{equation} \epsilon_a(T) = 84381.406000 - 46.836769T - 0.0001831T^2 + 0.00200340T^3 - 0.000000576T^4 - 0.0000000434T^5 \end{equation}

Look in the paper for the definitions of Delta psi and Delta epsilon.

I won't really go into the matrix T as Jai Willems answer and the discussion beneath should be sufficient and this response is already dragging on.

Polar motion is pretty negligible in most cases. Example in matlab you need to tell ecef2eci to account for it manually as it doesn't do it by default. If you're really curious I won't drag on my answer any longer and tell you to look here. It completely relies on up to date observation data that could be found there.

As for ECEF to ECI you compute all the matricies the same way and take their transpose as they're orthogonal.

To anyone crazy enough to read the insanity I went through and still want to compute ECI to ECEF, I hope you find this helpful, and best of luck.

Answer 3

The ECI and ECEF frames have approximately the same origin and $z$ axis and only differ by an angular component on the $xy$ plane. This angular difference can be used to rotate a vector from one frame to the other.

The $x_{ECEF}$ axis rotates along the equatorial plane about $z_{ECI}=z_{ECEF}$. The angle that it makes with the $x_{ECI}$ axis is known as the Earth Rotation Angle (ERA) or Greenwich Sidereal Angle. To convert between ECI and ECEF frames, we must first gain the ERA denoted as $\gamma\in[0, 2\pi]$.

We begin at an epoch, say J2000, where the ERA is known to be $280.46$ degrees. Using the elapsed time since epoch, we can propagate the epoch ERA through time to a date of interest. Stated mathematically in the following block equation, the Earth rotation rate is given by $360.985...$ degrees per day, $\Delta T$ is the elapsed time in days since the epoch, and $280.46$ is the ERA at $\Delta T=0$.

\begin{equation}\label{eqERA} \gamma = 360.9856123035484\times\Delta T + 280.46 \end{equation}

To convert a vector between the ECI and ECEF frames requires rotating the vector about the $z$-axis by the Greenwich Sidereal Angle, $\gamma$, using the rotational matrix given in the following block equation. These conversions can be seen in equations the final two equations where $\mu_{C'}^C$ converts coordinates from frame $C$ to frame $C'$.

\begin{equation}\label{eqEz} E_3^\theta = \begin{pmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} \end{equation}

\begin{equation}\label{eqECEF2ECI} \vec{V}_{ECI} = \mu_{ECI}^{ECEF}\;\vec{V}_{ECEF} = E_3^\gamma\;\vec{V}_{ECEF} \end{equation} \begin{equation}\label{eqECI2ECEF} \vec{V}_{ECEF} = \mu_{ECEF}^{ECI}\;\vec{V}_{ECI} = E_3^{-\gamma}\;\vec{V}_{ECI} \end{equation}

Answer 4

Example code to calculate the ECI to ECEF rotation matrix is given below.

It's based on the algorithms given in the ESA "Earth Observation Mission CFI Software Conventions Document" but with a very simple approximation for the nutation series taken from the 3rd edition of the "Explanatory Supplement to the Astronomical Almanac".

Information about the code - (a) it's written in IDL which is easy to read and understand, but note that comment lines start with a semicolon and a 'd' is used with floating point numbers to indicate that they are double precision rather than the default single precision. (b) arrays with 9 elements are used for the 3 by 3 transformation matrices (c) the code has been tested against the ESA EOCFI software with very good results (d) the algorithms use UT1 time (in MJD 2000 format which is days after 2000-01-01 00:00:00) (e) it's common to use TAI or GPS for time values rather than UTC (to avoid leap second problems) and you can convert from TAI and UTC to UT1 using the information in IERS bulletin A (which also contains polar motion data if you want to use it in the calculation). Probably you could use UTC instead of UT1 with only a small loss of accuracy. (f) the algorithms use the "geocentric mean of 2000" (GM2000), "mean of date" (MOD), "true of date" (TOD), "pseudo Earth fixed" (PEF) and "Earth fixed" (EF) coordinate systems, but GM2000 is just the usual ECI (J2000) coordinate system and the others are obtained by successive corrections for precession, nutation, Earth rotation and polar motion. (g) do not rely on this code, it's just an example to show that calculation of a reasonably accurate transformation matrix can be done in approximately 50 lines of code.

Other information - (a) The full IAU 2000A nutation model contains nearly 1400 terms and is described in "United States Naval Observatory Circular No. 179". The simpler IAU 2000B model contains only 77 terms and code implementing it can be found in the IAU SOFA library. (b) the EOCFI software can be downloaded from the ESA website, but it's difficult to use and the full documentation is over 1000 pages.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; calculate eci to ecef rotation matrix
;   input - time (ut1 in mjd2000)
;   output - rotation matrix
; based on "EOCFI mission conventions"
; doc (EO-MA-DMS-GS-0001, v4.12) with
; simple nutation model from the 3rd
; edition of the "Explanatory Supplement
; to the Astronomical Almanac"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function calc_eci2ecef, ut1

  ; rotation matrix is W T N P
  ; where
  ;   P is precession matrix
  ;   N is nutation matrix
  ;   T is Earth rotation
  ;   W is polar motion
  ; the effect of polar motion is
  ; small and W is neglected

  ; constant for converting degrees to radians
  deg2rad = !dpi/180d

  ; precession (gm2000 to mod)
  t = (ut1 - 0.5d)/36525d
  zeta = t*( 0.6406161d + t*( 0.0000839d + 0.0000050d*t ) )
  z = t*( 0.6406161d + t*( 0.0003041d + 0.0000051d*t ) )
  theta = t*( 0.5567530d - t*( 0.0001185d + 0.0000116d*t ) )
  a = calc_rz( -!dpi*0.5d - z*deg2rad )
  b = calc_rx( theta*deg2rad )
  c = calc_rz( !dpi*0.5d - zeta*deg2rad )
  m1 = mult(a,mult(b,c))

  ; nutation (mod to tod)
  t = ut1 - 0.5d
  a = ( 125.0d - 0.05295d*t )*deg2rad
  b = ( 200.9d + 1.97129d*t )*deg2rad
  dpsi = ( -0.0048d*sin(a) - 0.0004d*sin(b) )*deg2rad
  deps = ( 0.0026d*cos(a) + 0.0002d*cos(b) )*deg2rad
  eps = 23.439291d*deg2rad
  dmu = dpsi*cos(eps)
  dnu = dpsi*sin(eps)
  a = calc_rz( -dmu )
  b = calc_rx( -deps )
  c = calc_ry( dnu )
  m2 = mult(a,mult(b,c))

  ; earth rotation (tod to pef)
  t = ut1
  g = 99.96779469d + t*( 360.9856473662860d + 0.29079d-12*t )
  h = g*deg2rad + dmu
  m3 = calc_rz(h)

  ; multiply the matrices
  ans = mult(mult(m3,m2),m1)
  return, ans
end

; rotation matrices around x, y and z axes

function calc_rx, theta
  cs = cos(theta)
  sn = sin(theta)
  return, [ 1d, 0, 0, 0, cs, sn, 0, -sn, cs ]
end

function calc_ry, theta
  cs = cos(theta)
  sn = sin(theta)
  return, [ cs, 0, -sn, 0, 1, 0, sn, 0, cs ]
end

function calc_rz, theta
  cs = cos(theta)
  sn = sin(theta)
  return, [ cs, sn, 0, -sn, cs, 0, 0, 0, 1 ]
end

; matrix multiplication

function mult, a, b
  ans = dblarr(9)
  x = a[0] & y = a[1] & z = a[2]
  ans[0] = x*b[0] + y*b[3] + z*b[6]
  ans[1] = x*b[1] + y*b[4] + z*b[7]
  ans[2] = x*b[2] + y*b[5] + z*b[8]
  x = a[3] & y = a[4] & z = a[5]
  ans[3] = x*b[0] + y*b[3] + z*b[6]
  ans[4] = x*b[1] + y*b[4] + z*b[7]
  ans[5] = x*b[2] + y*b[5] + z*b[8]
  x = a[6] & y = a[7] & z = a[8]
  ans[6] = x*b[0] + y*b[3] + z*b[6]
  ans[7] = x*b[1] + y*b[4] + z*b[7]
  ans[8] = x*b[2] + y*b[5] + z*b[8]
  return, ans
end

Answer 5

Two issues I found with the code provided here:

1. The following line really threw me as I was converting it to C. I now presume that 0.29079d-12 is equivalent to 0.29079e-12 in C syntax.

So my C code for this line now reads:

gg = 99.96779469 + t*(360.9856473662860 + 0.29079e-12 * t);

I checked this in the ESA document and it is correct.

2. I found that the Earth Rotation was in the wrong direction! So I modified the line

m3 = calc_rz(h)

to pass -h instead of h, as follows:

m3 = calc_rz(-h)