Geology Reference
In-Depth Information
C
PROGRAMME MEMSPEC.FOR
C
C MEMSPEC.FOR finds the maximum entropy method (MEM) spectral density for
C the VLBI pole path, interpolated on to 3513 equally spaced time points
C by natural spline interpolation, for a record segment of 26.25 years
C in length. The length of the prediction error filter is entered as input
C on request.
C
IMPLICIT DOUBLE COMPLEX(A-H,O-Z)
DIMENSION GT(4000),G(4000),PEF(4000),PER(4000),DFT(4000)
DOUBLE PRECISION DAY(4000),AM,PEP,PI,OMEGA,T,DT,AIM1,FREQ(4000),
1 AJ,ARG,SME(4000),DF
OPEN(UNIT=1,FILE='ppintp.dat',STATUS='OLD')
OPEN(UNIT=2,FILE='memspec.dat',STATUS='UNKNOWN')
C Set number of interpolated pole positions.
M=3513
AM=DFLOAT(M)
C Set values of pi and 2pi.
PI=3.141592653589793D0
OMEGA=2.D0*PI
C Set record length.
T=26.25D0
C Find sample interval.
DT=T/AM
C Set mean of pole path to zero.
GTM=(0.D0,0.D0)
C Set prediction error power to zero.
PEP=0.D0
C Read in interpolated pole path.
DO 10 I=1,M
READ(1,11)J,IYYY,IM,ID,DAY(I),GT(I)
11 FORMAT(1X,4I5,F12.7,2F12.6)
C Accumulate sum of co-ordinates.
GTM=GTM+GT(I)
10 CONTINUE
C Find mean of pole positions.
GTM=GTM/AM
C Reduce co-ordinates to zero mean.
DO 12 I=1,M
GT(I)=GT(I)-GTM
C Find sum of squares of pole positions.
PEP=PEP+GT(I)*DCONJG(GT(I))
12 CONTINUE
C Find initial value of prediction error power.
PEP=PEP/AM
C Enter length of prediction error filter.
WRITE(6,13)
13 FORMAT(1X,'Enter length of prediction error filter.')
READ(5,*)K
C Find prediction error filter.
DO 14 NP1=1,K
CALL BPEC(M,NP1,GT,G,PEF,PER)
C Update prediction error power.
PEP=PEP*(1.D0-G(NP1)*DCONJG(G(NP1)))
14 CONTINUE
C Find DFT of prediction error filter.
DF=0.01D0
DO 15 I=1,301
IM1=I-1
AIM1=DFLOAT(IM1)
FREQ(I)=-1.5D0+AIM1*DF
Search WWH ::
Custom Search