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

C Calculate length of data set N.
N=IE-IB+1
C Find length of DFT, M, an odd integer less than or equal to N.
M=N
NB2=N/2
IF(N.EQ.2*NB2)M=N-1
C Find number of positive harmonics.
MP=(M-1)/2
C Find equivalent time sample interval.
AM=DFLOAT(M)
DT=T/AM
C Read in unused initial record segment.
IBM1=IB-1
DO 15 I=1,IBM1
READ(1,16)K,IYYY,MM,ID,TD,XWOB,YWOB,EXWOB,EYWOB
16 FORMAT(1X,4I5,F10.6,4F11.6)
15 CONTINUE
C Set sums of real and imaginary parts of co-ordinates to zero.
SUMX=0.D0
SUMY=0.D0
C Read in segment data set.
DO 17 I=1,N
READ(1,16)K,IYYY,MM,ID,DAY(I),XWOB,YWOB,EXWOB,EYWOB
C Accumulate sums of real and imaginary parts of co-ordinates.
SUMX=SUMX+XWOB
SUMY=SUMY+YWOB
C Convert co-ordinates to complex form.
G(I)=DCMPLX(XWOB,YWOB)
C Find variance.
SIGS(I)=EXWOB*EXWOB+EYWOB*EYWOB
17 CONTINUE
C Find means of real and imaginary parts of co-ordinates.
AN=DFLOAT(N)
SUMX=SUMX/AN
SUMY=SUMY/AN
C Reduce co-ordinates to zero mean and calculate their sum of squares.
RMS=0.D0
DO 18 I=1,N
G(I)=G(I)-DCMPLX(SUMX,SUMY)
RMS=RMS+G(I)*DCONJG(G(I))
C Window the data segment with the Parzen window.
CALL WINDOW(DAY(I),T0,T,W)
G(I)=G(I)*W
18 CONTINUE
C Find left side of Parseval's relation.
SUMS=RMS*DT
C Find root mean square of signal.
RMS=RMS/AN
RMS=DSQRT(RMS)
C Construct first column, A(J,1), of Toeplitz conditional equations matrix,
C and right-hand side vector B(J).
DO 19 J=1,M
AJ1=DFLOAT(J-1)
AJ2=DFLOAT(J-MP-1)
A(J,1)=(0.D0,0.D0)
B(J)=(0.D0,0.D0)
DO 20 I=1,N
C Form arguments of complex exponentials.
ARG=-2.D0*PI*DAY(I)/T
ARG1=ARG*AJ1
ARG2=ARG*AJ2