Earth’s rotation: observations and theory - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 279

With the Julian day input as JULIAN, the subroutine outputs the corresponding
year, month and day as IYYY, MM, ID, given that a specified Julian day starts at
noon on a calendar day. All of these conversions are applied by the programme
PPCY2009A listed below. The output file ppy2009a.dat contains 3560 pole posi-
tions and their standard errors running from 3 August 1979 to 22 October 2009.

C PROGRAMME PPCY2009A.FOR

C
C PPCY2009A.FOR converts pole co-ordinates to centiseconds of arc subtended
C at the geocentre, reverses the sign of the x2 co-ordinate to change
C to a right-handed co-ordinate system, changes the time base to years
C and converts Julian days to calendar days. Input is pp2009a.dat,
C extracted from the file 2009a.eops on the GSFC website with null values
C eliminated and time values differing by less than a day averaged.
C Output is ppy2009a.dat.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
OPEN(UNIT=1,FILE='pp2009a.dat',STATUS='OLD')
OPEN(UNIT=2,FILE='ppy2009a.dat',STATUS='UNKNOWN')
C Read in VLBI pole co-ordinates with null values eliminated and values
C with identical time tags averaged.
DO 10 I=1,3560
READ(1,11)J,DAY,XWOB,EXWOB,YWOB,EYWOB
11 FORMAT(1X,I4,F16.6,4F10.6)
C Change Julian days from beginning of the record to years.
TOLD=DAY/365.25636D0
C Change sign of x2 co-ordinate, scale up to centiseconds of arc.
XWOB=100.D0*XWOB
EXWOB=100.D0*EXWOB
YWOB=-100.D0*YWOB
EYWOB=100.D0*EYWOB
C Change Julian days to calendar days as year, month, day.
C Add initial Julian day to Julian days from beginning of the record.
DAY=DAY+44089.993750D0+2400000.D0
IDAY=INT(DAY)
CALL CALDAT(IDAY,MM,ID,IYYY)
C Write out converted times and co-ordinates.
WRITE(2,12)J,IYYY,MM,ID,TOLD,XWOB,YWOB,EXWOB,EYWOB

12 FORMAT(1X,4I5,F10.6,4F11.6)

10 CONTINUE

END

As an illustration of the converted pole path, we show a plot of the pole path for
the year 2004 in Figure 4.3. The standard error bars inflated by a factor of 100 are
overplotted to show estimated uncertainties.

4.2.2 Spectral analysis of the VLBI pole path
The co-ordinates of the VLBI pole path are given at unequally spaced times t j with
accompanying standard errors. To calculate the discrete Fourier transform (DFT) of
a segment of the VLBI pole path we then apply the methods of Section 2.3.5. If the
sequence of complex pole positions is g j = x 1 j + ix 2 j , j = 1,..., L , it is represented

by the sum (2.178) of 2 N + 1 complex sinusoids,