Geology Reference
In-Depth Information
analysis (MEM), it is useful to have the pole positions uniformly spaced in time.
Since the co-ordinates of the VLBI pole path are given at unequally spaced times,
this will require interpolation.
One possibility is to invert the DFT found by singular value decomposition.
While this DFT quite accurately reproduces the original unequally spaced pole
positions, the large number of harmonics generate noise away from the measured
positions. Instead, we interpolate between the measured positions using natural
cubic splines, as described in Section 1.6. Again, we take the part of the record to
be interpolated to begin at 4.0 years into the VLBI series and to end at 30.25 years,
for a segment length of 26.25 years. There are then 3513 internal nodes at time
points 48 to 3560. Function values at the beginning and end of the record segment
are found by assuming constant slope in the first two and last two intervals. The
spline interpolation produces 3513 equally spaced pole positions at intervals of
26.25
3513 years, or approximately 2.73 days. The programme PPINTP performs
natural spline interpolation with calls to the subroutines SPMAT and INTPL. The
corresponding calendar dates are found by calls to the subroutine CALDAT. The
output file ppintp.dat contains the 3513 equally spaced pole positions as well as
their times and calendar dates.
/
C
PROGRAMME PPINTP.FOR
C
C PPINTP.FOR does a natural spline interpolation of the VLBI pole path
C onto 3513 equally spaced time points. The record segment is set
C to begin at 4.0 years into the VLBI series and to end at 30.25 years,
C for a segment length of 26.25 years. There are then 3513 internal nodes
C at time points 48 to 3560. Values at the beginning and end of
C the record segment are found by assuming constant slope in the first two
C and last two intervals. The corresponding calendar dates are found by
C calls to the subroutine CALDAT.
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION DAY(4000),XWOB(4000),YWOB(4000),C(4000,4000),
1 B(4000,8000),GTX(4000),GTY(4000)
OPEN(UNIT=1,FILE='ppy2009a.dat',STATUS='OLD')
OPEN(UNIT=2,FILE='ppintp.dat',STATUS='UNKNOWN')
C Set number of points, N, and number of interpolates, M.
N=3513
M=3513
AM=DFLOAT(M)
C Set segment length.
T=26.25D0
C Find time increment.
DT=T/AM
C Read in unused initial record segment, points 1 to 47.
DO 10 I=1,47
READ(1,11)K,IYYY,MM,ID,TD,X,Y,EX,EY
11 FORMAT(1X,4I5,F10.6,4F11.6)
10 CONTINUE
C Read in segment to be interpolated, points 48 to 3560.
DO 12 I=1,N
IP1=I+1
READ(1,11)K,IYYY,MM,ID,DAY(IP1),XWOB(IP1),YWOB(IP1),EXWOB,EYWOB
Search WWH ::
Custom Search