*  SOFA examples

      IMPLICIT NONE

*  Arcseconds to radians
      DOUBLE PRECISION AS2R
      PARAMETER ( AS2R = 4.848136811095359935899141D-6 )

*  2Pi
      DOUBLE PRECISION D2PI
      PARAMETER ( D2PI = 6.283185307179586476925287D0 )

      INTEGER IY, IM, ID, IH, MIN, J
      DOUBLE PRECISION SEC, XP, YP, DUT1,
     :                 DDP80, DDE80, DX00, DY00, DX06, DY06,
     :                 DJMJD0, DATE, TIME, UTC, DAT,
     :                 TAI, TT, TUT, UT1, RP(3,3), DP80, DE80,
     :                 DPSI, DEPS, EPSA, RN(3,3), RNPB(3,3),
     :                 EE, GST, RC2TI(3,3), RPOM(3,3),
     :                 RC2IT(3,3), X, Y, S,
     :                 RC2I(3,3), ERA, DP00, DE00, RB(3,3),
     :                 RPB(3,3), V1(3), V2(3), DDP00, DDE00
      DOUBLE PRECISION iau_OBL80, iau_EQEQ94, iau_ANP, iau_GMST82,
     :                 iau_ERA00, iau_SP00, iau_EE00, iau_GMST00,
     :                 iau_S06


*  UTC.
      IY = 2007
      IM = 4
      ID = 5
      IH = 12
      MIN = 0
      SEC = 0D0

*  Polar motion (arcsec->radians).
      XP = 0.0349282D0 * AS2R
      YP = 0.4833163D0 * AS2R

*  UT1-UTC (s).
      DUT1 = -0.072073685D0

*  Nutation corrections wrt IAU 1976/1980 (mas->radians).
      DDP80 = -55.0655D0 * AS2R/1000D0
      DDE80 =  -6.3580D0 * AS2R/1000D0

*  CIP offsets wrt IAU 2000A (mas->radians).
      DX00 =  0.1725D0 * AS2R/1000D0
      DY00 = -0.2650D0 * AS2R/1000D0

*  CIP offsets wrt IAU 2006/2000A (mas->radians).
      DX06 =  0.1750D0 * AS2R/1000D0
      DY06 = -0.2259D0 * AS2R/1000D0

*  TT (MJD).
      CALL iau_CAL2JD ( IY, IM, ID, DJMJD0, DATE, J )
      TIME = ( 60D0*(60D0*DBLE(IH) + DBLE(MIN)) + SEC ) / 86400D0
      UTC = DATE + TIME
      CALL iau_DAT ( IY, IM, ID, TIME, DAT, J )
      TAI = UTC + DAT/86400D0
      TT = TAI + 32.184D0/86400D0

*  UT1.
      TUT = TIME + DUT1/86400D0
      UT1 = DATE + TUT

*  =============
*  IAU 1976/1980
*  =============

*  IAU 1976 precession matrix, J2000.0 to date.
      CALL iau_PMAT76 ( DJMJD0, TT, RP )

*  IAU 1980 nutation.
      CALL iau_NUT80 ( DJMJD0, TT, DP80, DE80 )

*  Add adjustments: frame bias, precession-rates, geophysical.
      DPSI = DP80 + DDP80
      DEPS = DE80 + DDE80

*  Mean obliquity.
      EPSA = iau_OBL80 ( DJMJD0, TT )

*  Build the rotation matrix.
      CALL iau_NUMAT ( EPSA, DPSI, DEPS, RN )

*  Combine the matrices:  PN = N x P.
      CALL iau_RXR ( RN, RP, RNPB )

*  Equation of the equinoxes, including nutation correction.
      EE = iau_EQEQ94 ( DJMJD0, TT ) + DDP80 * COS ( EPSA )

*  Greenwich apparent sidereal time (IAU 1982/1994).
      GST = iau_ANP ( iau_GMST82 ( DJMJD0+DATE, TUT ) + EE )

*  Form celestial-terrestrial matrix (no polar motion yet).
      CALL iau_CR ( RNPB, RC2TI )
      CALL iau_RZ ( GST, RC2TI )

*  Polar motion matrix (TIRS->ITRS, IERS 1996).
      CALL iau_IR ( RPOM )
      CALL iau_RX ( -YP, RPOM )
      CALL iau_RY ( -XP, RPOM )

*  Form celestial-terrestrial matrix (including polar motion).
      CALL iau_RXR ( RPOM, RC2TI, RC2IT )

*  ====================
*  IAU 2000A, CIO based
*  ====================

*  CIP and CIO, IAU 2000A.
      CALL iau_XYS00A ( DJMJD0, TT, X, Y, S )

*  Add CIP corrections.
      X = X + DX00
      Y = Y + DY00

*  GCRS to CIRS matrix.
      CALL iau_C2IXYS ( X, Y, S, RC2I )

*  Earth rotation angle.
      ERA = iau_ERA00 ( DJMJD0+DATE, TUT )

*  Form celestial-terrestrial matrix (no polar motion yet).
      CALL iau_CR ( RC2I, RC2TI )
      CALL iau_RZ ( ERA, RC2TI )

*  Polar motion matrix (TIRS->ITRS, IERS 2003).
      CALL iau_POM00 ( XP, YP, iau_SP00(DJMJD0,TT), RPOM )

*  Form celestial-terrestrial matrix (including polar motion).
      CALL iau_RXR ( RPOM, RC2TI, RC2IT )

*  ========================
*  IAU 2000A, equinox based
*  ========================

*  Nutation, IAU 2000A.
      CALL iau_NUT00A ( DJMJD0, TT, DP00, DE00 )

*  Precession-nutation quantities, IAU 2000.
      CALL iau_PN00 ( DJMJD0, TT, DP00, DE00,
     :                EPSA, RB, RP, RPB, RN, RNPB )

*  Transform dX,dY corrections from GCRS to mean of date.
      V1(1) = DX00
      V1(2) = DY00
      V1(3) = 0D0
      CALL iau_RXP ( RNPB, V1, V2 )
      DDP00 = V2(1) / SIN ( EPSA )
      DDE00 = V2(2)

*  Corrected nutation.
      DPSI = DP00 + DDP00
      DEPS = DE00 + DDE00

*  Build the rotation matrix.
      CALL iau_NUMAT ( EPSA, DPSI, DEPS, RN )

*  Combine the matrices:  PN = N x P.
      CALL iau_RXR ( RN, RPB, RNPB )

*  Greenwich apparent sidereal time (IAU 2000).
      GST = iau_ANP ( iau_GMST00 ( DJMJD0+DATE, TUT, DJMJD0, TT )
     :                + iau_EE00 ( DJMJD0, TT, EPSA, DPSI ) )

*  Form celestial-terrestrial matrix (no polar motion yet).
      CALL iau_CR ( RNPB, RC2TI )
      CALL iau_RZ ( GST, RC2TI )

*  Polar motion matrix (TIRS->ITRS, IERS 2003).
      CALL iau_POM00 ( XP, YP, iau_SP00(DJMJD0,TT), RPOM )

*  Form celestial-terrestrial matrix (including polar motion).
      CALL iau_RXR ( RPOM, RC2TI, RC2IT )

*  =========================
*  IAU 2006/2000A, CIO based
*  =========================

*  CIP and CIO, IAU 2006/2000A.
      CALL iau_XYS06A ( DJMJD0, TT, X, Y, S )

*  Add CIP corrections.
      X = X + DX06
      Y = Y + DY06

*  GCRS to CIRS matrix.
      CALL iau_C2IXYS ( X, Y, S, RC2I )

*  Earth rotation angle.
      ERA = iau_ERA00 ( DJMJD0+DATE, TUT )

*  Form celestial-terrestrial matrix (no polar motion yet).
      CALL iau_CR ( RC2I, RC2TI )
      CALL iau_RZ ( ERA, RC2TI )

*  Polar motion matrix (TIRS->ITRS, IERS 2003).
      CALL iau_POM00 ( XP, YP, iau_SP00(DJMJD0,TT), RPOM )

*  Form celestial-terrestrial matrix (including polar motion).
      CALL iau_RXR ( RPOM, RC2TI, RC2IT )

*  ===========================================
*  IAU 2006/2000A, CIO based, using X,Y series
*  ===========================================

*  CIP and CIO, IAU 2006/2000A.
      CALL iau_XY06 ( DJMJD0, TT, X, Y )
      S = iau_S06 ( DJMJD0, TT, X, Y )

*  Add CIP corrections.
      X = X + DX06
      Y = Y + DY06

*  GCRS to CIRS matrix.
      CALL iau_C2IXYS ( X, Y, S, RC2I )

*  Earth rotation angle.
      ERA = iau_ERA00 ( DJMJD0+DATE, TUT )

*  Form celestial-terrestrial matrix (no polar motion yet).
      CALL iau_CR ( RC2I, RC2TI )
      CALL iau_RZ ( ERA, RC2TI )

*  Polar motion matrix (TIRS->ITRS, IERS 2003).
      CALL iau_POM00 ( XP, YP, iau_SP00(DJMJD0,TT), RPOM )

*  Form celestial-terrestrial matrix (including polar motion).
      CALL iau_RXR ( RPOM, RC2TI, RC2IT )

      END
