Index: /trunk/doc/pslib/ChangeLogADD.tex =================================================================== --- /trunk/doc/pslib/ChangeLogADD.tex (revision 3428) +++ /trunk/doc/pslib/ChangeLogADD.tex (revision 3429) @@ -28,3 +28,4 @@ \item Added section on FITS WCS. \item Renamed TDT to TT plus misc. cleanups in the Time section. +\item Updated definition of UT1. \end{itemize} Index: /trunk/doc/pslib/interp.f =================================================================== --- /trunk/doc/pslib/interp.f (revision 3429) +++ /trunk/doc/pslib/interp.f (revision 3429) @@ -0,0 +1,165 @@ + SUBROUTINE INTERP (RJD,X,Y,T,N,RJDINT,XINT,YINT,TINT) +C +C THIS SUBROUTINE TAKES A SERIES OF X, Y, AND UT1-UTC VALUES +C AND INTERPOLATES THEM TO AN EPOCH OF CHOICE. THIS ROUTINE +C ASSUMES THAT THE VALUES OF X AND Y ARE IN SECONDS OF +C ARC AND THAT UT1-UTC IS IN SECONDS OF TIME. AT LEAST +C ONE POINT BEFORE AND ONE POINT AFTER THE EPOCH OF THE +C INTERPOLATION POINT ARE NECESSARY IN ORDER FOR THE +C INTERPOLATION SCHEME TO WORK. +C +C PARAMETERS ARE : +C RJD - ARRAY OF THE EPOCHS OF DATA (GIVEN IN MJD) +C X - ARRAY OF X POLAR MOTION (ARCSEC) +C Y - ARRAY OF Y POLAR MOTION (ARCSEC) +C T - ARRAY OF UT1-UTC (SEC) +C N - NUMBER OF POINTS IN ARRAYS +C RJDINT- EPOCH FOR THE INTERPOLATED VALUE +C XINT - INTERPOLATED VALUE OF X +C YINT - INTERPOLATED VALUE OF Y +C TINT - INTERPOLATED VALUE OF UT1-UTC +C + DOUBLE PRECISION RJD(N), X(N), Y(N), T(N), + . RJDINT, XINT, YINT, TINT, CORX, CORY, CORT + CALL LAGINT (RJD,X,N,RJDINT,XINT) + CALL LAGINT (RJD,Y,N,RJDINT,YINT) + CALL LAGINT (RJD,T,N,RJDINT,TINT) + CALL RAY (RJDINT,CORX,CORY,CORT) + XINT = XINT + CORX + YINT = YINT + CORY + TINT = TINT + CORT + RETURN + END +C +C---------------------------------------------------------------- +C + SUBROUTINE LAGINT (X,Y,N,XINT,YOUT) +C +C THIS SUBROUTINE PERFORMS LAGRANGIAN INTERPOLATION +C WITHIN A SET OF (X,Y) PAIRS TO GIVE THE Y +C VALUE CORRESPONDING TO XINT. THIS PROGRAM USES A +C WINDOW OF 4 DATA POINTS TO PERFORM THE INTERPOLATION. +C IF THE WINDOW SIZE NEEDS TO BE CHANGED, THIS CAN BE +C DONE BY CHANGING THE INDICES IN THE DO LOOPS FOR +C VARIABLES M AND J. +C +C PARAMETERS ARE : +C X - ARRAY OF VALUES OF THE INDEPENDENT VARIABLE +C Y - ARRAY OF FUNCTION VALUES CORRESPONDING TO X +C N - NUMBER OF POINTS +C XINT - THE X-VALUE FOR WHICH ESTIMATE OF Y IS DESIRED +C YOUT - THE Y VALUE RETURNED TO CALLER +C + REAL*8 X(N),Y(N),XINT,YOUT,TERM + INTEGER N,I,J +C + YOUT = 0.0D0 + DO 5 I = 1,N-1 + IF ( XINT .GE. X(I) .AND. XINT .LT. X(I+1) ) K = I + 5 CONTINUE + IF ( K .LT. 2 ) K = 2 + IF ( K .GT. N-2 ) K = N-2 + DO 20 M = K-1,K+2 + TERM = Y(M) + DO 10 J = K-1,K+2 + IF ( M .NE. J ) THEN + TERM = TERM * (XINT - X(J))/(X(M) - X(J)) + END IF + 10 CONTINUE + YOUT = YOUT + TERM + 20 CONTINUE + RETURN + END + +C +C---------------------------------------------------------------- +C + SUBROUTINE RAY (RJD,CORX,CORY,CORT) +C +C THIS SUBROUTINE IMPLEMENTS THE RAY MODEL FOR +C DIURNAL/SUBDIURNAL TIDES. IT USES THE SIMON ET AL. +C FUNDAMENTAL ARGUMENTS. THE CORRECTIONS IN X AND Y ARE IN +C UNITS OF SEC. OF ARC AND UT1-UTC IN SEC. OF TIME. THESE +C CORRECTIONS SHOULD BE ADDED TO "AVERAGE" EOP VALUES TO GET +C ESTIMATES OF THE INSTANTANEOUS VALUES. +C +C PARAMETERS ARE : +C RJD - EPOCH OF INTEREST GIVEN IN MJD +C CORX - TIDAL CORRECTION IN X (SEC. OF ARC) +C CORY - TIDAL CORRECTION IN Y (SEC. OF ARC) +C CORT - TIDAL CORRECTION IN UT1-UTC (SEC. OF TIME) +C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DOUBLE PRECISION + . L, LPRIME + HALFPI = 1.5707963267948966d0 + T = (RJD - 51544.5D0)/36525.0D0 + L = -0.00024470d0*T**4 + 0.051635d0*T**3 + 31.8792d0*T**2 + . + 1717915923.2178d0*T + 485868.249036d0 + L = DMOD(L,1296000d0) + LPRIME = -0.00001149d0*T**4 - 0.000136d0*T**3 + . - 0.5532d0*T**2 + . + 129596581.0481d0*T + 1287104.79305d0 + LPRIME = DMOD(LPRIME,1296000d0) + CAPF = 0.00000417d0*T**4 - 0.001037d0*T**3 - 12.7512d0*T**2 + . + 1739527262.8478d0*T + 335779.526232d0 + CAPF = DMOD(CAPF,1296000d0) + CAPD = -0.00003169d0*T**4 + 0.006593d0*T**3 - 6.3706d0*T**2 + . + 1602961601.2090d0*T + 1072260.70369d0 + CAPD = DMOD(CAPD,1296000d0) + OMEGA = -0.00005939d0*T**4 + 0.007702d0*T**3 + . + 7.4722d0*T**2 + . - 6962890.2665d0*T + 450160.398036d0 + OMEGA = DMOD(OMEGA,1296000d0) + THETA = (67310.54841d0 + + . (876600d0*3600d0 + 8640184.812866d0)*T + + . 0.093104d0*T**2 - + . 6.2d-6*T**3)*15.0d0 + 648000.0d0 + ARG7 = DMOD((-L - 2.0D0*CAPF - 2.0D0*OMEGA + THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI + ARG1 = DMOD((-2.0d0*CAPF - 2.0d0*OMEGA + THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI + ARG2 = DMOD((-2.0d0*CAPF + 2.0d0*CAPD - 2.0d0*OMEGA + . + THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI + ARG3 = DMOD(THETA * + . 3.14159265D0/648000.0D0,6.28318530718D0) + . + HALFPI + ARG4 = DMOD((-L - 2.0d0*CAPF - 2.0D0*OMEGA + 2.0d0*THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) + ARG5 = DMOD((-2.0D0*CAPF - 2.0D0*OMEGA + 2.0d0*THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) + ARG6 = DMOD((-2.0d0*CAPF + 2.0d0*CAPD - 2.0d0*OMEGA + . + 2.0d0*THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) + ARG8 = DMOD((2.0d0*THETA) + . * 3.14159265D0/648000.0D0,6.28318530718D0) + CORX = - 0.026D0*DSIN(ARG7) + 0.006D0*DCOS(ARG7) + . - 0.133D0*DSIN(ARG1) + 0.049D0*DCOS(ARG1) + . - 0.050D0*DSIN(ARG2) + 0.025D0*DCOS(ARG2) + . - 0.152D0*DSIN(ARG3) + 0.078D0*DCOS(ARG3) + . - 0.057D0*DSIN(ARG4) - 0.013D0*DCOS(ARG4) + . - 0.330D0*DSIN(ARG5) - 0.028D0*DCOS(ARG5) + . - 0.145D0*DSIN(ARG6) + 0.064D0*DCOS(ARG6) + . - 0.036D0*DSIN(ARG8) + 0.017D0*DCOS(ARG8) + CORY = - 0.006D0*DSIN(ARG7) - 0.026D0*DCOS(ARG7) + . - 0.049D0*DSIN(ARG1) - 0.133D0*DCOS(ARG1) + . - 0.025D0*DSIN(ARG2) - 0.050D0*DCOS(ARG2) + . - 0.078D0*DSIN(ARG3) - 0.152D0*DCOS(ARG3) + . + 0.011D0*DSIN(ARG4) + 0.033D0*DCOS(ARG4) + . + 0.037D0*DSIN(ARG5) + 0.196D0*DCOS(ARG5) + . + 0.059D0*DSIN(ARG6) + 0.087D0*DCOS(ARG6) + . + 0.018D0*DSIN(ARG8) + 0.022D0*DCOS(ARG8) + CORT = + 0.0245D0*DSIN(ARG7) + 0.0503D0*DCOS(ARG7) + . + 0.1210D0*DSIN(ARG1) + 0.1605D0*DCOS(ARG1) + . + 0.0286D0*DSIN(ARG2) + 0.0516D0*DCOS(ARG2) + . + 0.0864D0*DSIN(ARG3) + 0.1771D0*DCOS(ARG3) + . - 0.0380D0*DSIN(ARG4) - 0.0154D0*DCOS(ARG4) + . - 0.1617D0*DSIN(ARG5) - 0.0720D0*DCOS(ARG5) + . - 0.0759D0*DSIN(ARG6) - 0.0004D0*DCOS(ARG6) + . - 0.0196D0*DSIN(ARG8) - 0.0038D0*DCOS(ARG8) + CORX = CORX * 1.0d-3 + CORY = CORY * 1.0d-3 + CORT = CORT * 0.1d-3 + RETURN + END Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 3428) +++ /trunk/doc/pslib/psLibADD.tex (revision 3429) @@ -1,6 +1,7 @@ -%%% $Id: psLibADD.tex,v 1.64 2005-03-16 00:29:31 jhoblitt Exp $ +%%% $Id: psLibADD.tex,v 1.65 2005-03-16 01:49:30 jhoblitt Exp $ \documentclass[panstarrs]{panstarrs} \usepackage{amsmath} +\usepackage{verbatim} % basic document variables @@ -961,9 +962,44 @@ \label{sec:ut1} -Univseral Time is a measure of the rotation angle of the Earth. When -corrected for polar motion it is referred to as UT1. This is distict -from UT0 which does not involve corrections for polar motion. UT1 may -be calculated from UTC through a table lookup of the appropriate value -of UTC - UT1. +UT1 is directly tied to the rotation of the Earth. Historically, time has been +measured with respect to the rising and setting of the Sun. However, in the +modern era of atomic clocks, the rotation of the Earth makes for a highly +unstable time standard. Tidal effects, changes in the angular momentum of the +atmosphere, seasonal changes in the polar ice caps, movement within the Earth's +core, and other effects all cause measurable changes in the Earth's rotation on +a daily basis. However, UT1 is still vitally important for determining the +orientation of the Earth with respect to the sky. + +Since 2003-01-01, UT1 has been defined as directly proportional to the Earth +Rotation Angle (see IERS Technical Note 32\footnote{IERS Technical Note 32 - +http://maia.usno.navy.mil/conv2003.html}). Previous to that date, a different +definition was in effect (see IERS Technical Note 21\footnote{IERS Technical +Note 21 - http://maia.usno.navy.mil/conventions.html}). We will always use the +post-2003 definition. + +UT1 is continuously measured by the International Earth Rotation +Service\footnote{IERS - http://maia.usno.navy.mil/}, and tabulated values of the +offset of UT1 from UTC are published at regular intervals, along with predicted +future values. IERS Bulletin A gives "rapid response" values necessary for +real-time and near real-time data analysis (such as Pan-STARRS Otis and IPP +subsystems). Bulletin B gives the results of a final, definitive data +reduction. An amalgam of Bulletin A and B values is published daily on the +IERS website\footnote{IERS Bulletin A \& B - +http://maia.usno.navy.mil/ser7/finals2000A.daily} along with a desciption of +the format\footnote{IERS finals2000A.daily table format - +http://maia.usno.navy.mil/ser7/readme.finals2000A}. + +The UT1 offsets should be interpolated using the prescription of IERS Gazette +\#13\footnote{IERS Gazette \#13 - http://maia.usno.navy.mil/iers-gaz13}. This +entails using a third order polynomial to interpolate the table values, and +then applying a model for diurnal and semi-diurnal fluctuations due to tidal +effects. An example implimentation\footnote{interp.f +ftp://maia.usno.navy.mil/dist/interp.f} written in Fortran is available. + +\verbatiminput{interp.f} + +The Polar motion X, and Y coordinates are also important for determining the +orientation of the Earth with respect to the sky. This is also given in the +IERS publications references above, and should be interpolated in the same way. \paragraph{Julian Day and Modified Julian Day}