IAU Rotation Models
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
- 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.
- J. Williams, Fortran Astrodynamics Toolkit - iau_orientation_module [GitHub]