IAU Rotation Models

earthThe 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 Toolkitiau_orientation_module [GitHub]
Tagged with: , , , ,

Leave a Reply

Your email address will not be published. Required fields are marked *

*