Introduction and theoretical background - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 86

DO 11 K=1,N1
SUM1=SUM1+C(J,K)*S(K)
SUM2=SUM2+C(J+1,K)*S(K)
11 CONTINUE
C Find interpolate YV at XV.
YV=S(J)+(S(J+1)-S(J))*(XV-R(J))/(R(J+1)-R(J))

1 +(XV-R(J))*(XV-R(J+1))*(SUM1*(2.D0-(XV-R(J))/(R(J+1)-R(J)))

2 +SUM2*(1.D0+(XV-R(J))/(R(J+1)-R(J))))/6.D0

RETURN
END

erentiation of a tabulated function can easily be
accomplished with spline approximation. With some rearrangement, expression
(1.548) gives the derivative of the approximating function as

As well as interpolation, di ff

f

− f i

f i + 1 − f i ( x − x i + 1 )( x − x i )

s ( x )

=

x i +

i + 1

x i + 1 − 6 ( x i + 1 − x i )

1 f 2 f

(2 x − x i + 1 − x i )
6

x − x i
x i + 1 −

x − x i
x i + 1 −

+ + + − , (1.568)

i + 1 i

x i x i

for the subinterval containing x . The subroutine DERIV calculates the derivative at
x after SPMAT has been called:

SUBROUTINE DERIV(XV,YV,N1,C,R,S,M1)

C
C DERIV calculates the derivative YV at XV, using a cubic spline,
C given the function S(I) at N1 nodal points R(I). C is
C the coefficient matrix, from subroutine SPMAT, relating
C the second derivatives to the function values at the nodal points.
C M1 is the maximum dimension for N1 specified in the main programme.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION C(M1,M1),R(M1),S(M1)
NPM1=N1-1
J=NPM1
C Locate subinterval.
DO 10 L=1,NPM1
IF(XV.GE.R(L).AND.XV.LT.R(L+1)) J=L
10 CONTINUE
C Find second derivative SUM1 at beginning of subinterval
C and second derivative SUM2 at the end of subinterval.
SUM1=0.D0
SUM2=0.D0
DO 11 K=1,N1
SUM1=SUM1+C(J,K)*S(K)
SUM2=SUM2+C(J+1,K)*S(K)
11 CONTINUE
C Find derivative YV at XV.
YV=(S(J+1)-S(J))/(R(J+1)-R(J))

1 +(XV-R(J))*(XV-R(J+1))*(SUM2-SUM1)/(6.D0*(R(J+1)-R(J)))

2 +(2.D0*XV-R(J+1)-R(J))*(SUM1*(2.D0-(XV-R(J))/(R(J+1)-R(J)))

3 +SUM2*(1.D0+(XV-R(J))/(R(J+1)-R(J))))/6.D0

RETURN
END