Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Feb 08, 2015

IAU Rotation Models

earth

The IAU Working Group on Cartographic Coordinates and Rotational Elements (WGCCRE) is the keeper of official models that describe the cartographic coordinates and rotational elements of planetary bodies (such as the Earth, satellites, minor planets, and comets). Periodically, they release a report containing the coefficients to compute body orientations, based on the latest data available. These coefficients allow one to compute the rotation matrix from the ICRF frame to a body-fixed frame (for example, IAU Earth) by giving the direction of the pole vector and the prime meridian location as functions of time. An example Fortran module illustrating this for the IAU Earth frame is given below. The coefficients are taken from the 2009 IAU report [Reference 1]. Note that the IAU models are also available in SPICE (as "IAU_EARTH", "IAU_MOON", etc.). For Earth, the IAU model is not suitable for use in applications that require the highest possible accuracy (for that a more complex model would be necessary), but is quite acceptable for many applications.

module iau_orientation_module

use, intrinsic :: iso_fortran_env, only: wp => real64

implicit none

private

!constants:
real(wp),parameter :: zero = 0.0_wp
real(wp),parameter :: one = 1.0_wp
real(wp),parameter :: pi = acos(-one)
real(wp),parameter :: pi2 = pi/2.0_wp
real(wp),parameter :: deg2rad = pi/180.0_wp
real(wp),parameter :: sec2day = one/86400.0_wp
real(wp),parameter :: sec2century = one/3155760000.0_wp

public :: icrf_to_iau_earth

contains

!Rotation matrix from ICRF to IAU_EARTH
pure function icrf_to_iau_earth(et) result(rotmat)

implicit none

real(wp),intent(in) :: et !ephemeris time [sec from J2000]
real(wp),dimension(3,3) :: rotmat !rotation matrix

real(wp) :: w,dec,ra,d,t
real(wp),dimension(3,3) :: tmp

d = et * sec2day !days from J2000
t = et * sec2century !Julian centuries from J2000

ra = ( - 0.641_wp * t ) * deg2rad
dec = ( 90.0_wp - 0.557_wp * t ) * deg2rad
w = ( 190.147_wp + 360.9856235_wp * d ) * deg2rad

!it is a 3-1-3 rotation:
tmp = matmul( rotmat_x(pi2-dec), rotmat_z(pi2+ra) )
rotmat = matmul( rotmat_z(w), tmp )

end function icrf_to_iau_earth

!The 3x3 rotation matrix for a rotation about the x-axis.
pure function rotmat_x(angle) result(rotmat)

implicit none

real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]

real(wp) :: c,s

c = cos(angle)
s = sin(angle)

rotmat = reshape([one, zero, zero, &
                  zero, c, -s, &
                  zero, s, c],[3,3])

end function rotmat_x

!The 3x3 rotation matrix for a rotation about the z-axis.
pure function rotmat_z(angle) result(rotmat)

implicit none

real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]

real(wp) :: c,s

c = cos(angle)
s = sin(angle)

rotmat = reshape([ c, -s, zero, &
                   s, c, zero, &
                   zero, zero, one ],[3,3])

end function rotmat_z

end module iau_orientation_module

References

  1. B. A. Archinal, et al., "Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009", Celest Mech Dyn Astr (2011) 109:101-135.
  2. J. Williams, Fortran Astrodynamics Toolkit - iau_orientation_module [GitHub]