* 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