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

0

-10

-20

-30

-40

-50

-60

-70

-30 -20 -10 0 10 20 30 40

m 1 in 0.01 Greenwich

Figure 4.5 Interpolated pole path for 2004. Pole positions measured by the VLBI
technique are shown as crosses. Continuous curve is the interpolated path.

density estimate for the single segment is given by the normalised squared mag-
nitude of the discrete Fourier transform (2.350). The programme INTPSPEC reads
in the 3513 interpolated pole path positions and reduces them to zero mean. They
are then Parzen windowed to suppress finite record e ects by a call to the sub-

routine WINDOW. The generalised DFT (2.150) is found by a call to the subroutine
FFTN, which in turn calls the subroutine FFT2N. We use the DFT expressed by
(2.148) which is obtained from the generalised DFT through expressions (2.160).
The single segment spectral density estimate is given as a function of frequency in
the output file intpspec.dat. As before, the confidence interval is of fixed length on
logarithmic plots of the spectral density. For a 95% confidence interval the fixed
length is 2.1634. A histogram plot of the spectral density is shown in Figure 4.6.
C

ff

PROGRAMME INTPSPEC.FOR

C
C INTPSPEC.FOR finds the spectral density estimate for the VLBI pole path,
C interpolated onto 3513 equally spaced time points by natural spline
C interpolation, for a record segment of 26.25 years in length.
C The fast Fourier transform is found by calls to the subroutines FFTN
C and FFT2N and then converted to the DFT, after the sequence is reduced
C to zero mean, and Parzen windowed to prevent finite record effects by
C a call to the subroutine WINDOW. The spectral density is given by
C the normalized squared magnitude of the DFT. The output file intpspec.dat
C gives the spectral density estimate as a function of frequency.
C

IMPLICIT DOUBLE COMPLEX(A-H,O-Z)
DOUBLE PRECISION FREQ(4000),SD(4000),DAY(4000),OMEGA,ARG,AI,
1 PI,T,T0,AN,DT,ANP,W
DIMENSION GT(4000),WK1(20000),WK2(20000),GF(4000)
OPEN(UNIT=1,FILE='ppintp.dat',STATUS='OLD')
OPEN(UNIT=2,FILE='intpspec.dat',STATUS='UNKNOWN')
C Set values of pi and 2pi.