!用三次樣條函數(shù)S（X）對給定的結點（Xi,Yi）(i=1,2,…,N)進行分段插值。
 !下面是求出各段表達式中的系數(shù)Bi,Ci和Di的子程序。  
    SUBROUTINE SPL(N,X,Y,B,C,D);
    DIMENSION X(N),Y(N),B(N),C(N),D(N);
    NM1=N-1;
    IF(N<2)RETURN;
    IF(N<3)GOTO 60;
    D(1)=X(2)-X(1);
    C(2)=(Y(2)-Y(1))/D(1);
    DO 10 I=2,NM1;
    D(I)=X(I+1)-X(I);
    B(I)=2*(D(I-1)+D(I));
    C(I+1)=(Y(I+1)-Y(I))/D(I);
10  C(I)=C(I+1)-C(I);
    B(1)=-D(1);
    B(N)=-D(N-1);
    C(1)=0;
    C(N)=0;
    IF (N.EQ.3)GOTO 20;
    C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1));
    C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3));
    C(1)=C(1)*D(1)**2/(X(N)-X(N-3));
20  DO 30 I=2,N;
    T=D(I-1)/B(I-1);
    B(I)=B(I)-T*D(I-1);
30  C(I)=C(I)-T*C(I-1);
    C(N)=C(N)/B(N);
    DO 40 I1=1,NM1;
    I=N-I1;
40  C(I)=(C(I)-D(I)*c(I+1))/B(I);
    B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2*C(N));
    DO 50 I=1,NM1;
    B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2*C(I));
    D(I)=(C(I+1)-C(I))/D(I);
50  C(I)=3*C(I);
    C(N)=3*C(N);
    D(N)=D(N-1);
    RETURN;
60  B(1)=(Y(2)-Y(1))/(X(2)-X(1));
    C(1)=0;
    D(1)=0;
    B(2)=B(1);
    C(2)=0;
    D(2)=0;
    RETURN;
    END;
!本函數(shù)段對插值點x，調用三次樣條函數(shù)插值系數(shù)子程序后，可計算出其
!對應的函數(shù)值s(x)。
    FUNCTION SEVA(N,U,X,Y,B,C,D);
    DIMENSION X(N),Y(N),B(N),C(N),D(N);
    I=1;
    IF(I.GE.N)I=1;
    IF(U<X(I))GOTO 11;
    IF(U<=X(I+1))GOTO 31;
11  I=1;
    J=N+1;
21  K=INT((I+J)/2);
    IF(U<X(K))J=K;
    IF(U.GE.X(K))I=K;
    IF(J>(I+1))GOTO 21;
31  DX=U-X(I);
    SEVA=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)));
    RETURN;
    END;
!下面是本程序的主程序部分。  
    PROGRAM MAIN
    PARAMETER(L=11,M=400);
!上行中L是插值點數(shù)，M是地面上任意給定的自激自收點的橫坐標。
    DIMENSION X(L),Y(L),B(L),C(L),D(L),S(5001),E(5001),F(5001),Z(5001);
!上行中F數(shù)組是地下界面上的插值點到自激自收的距離。
    OPEN(110,FILE='2.TXT')
!輸入節(jié)點號及其相應的深度值。
    READ (110,*)X,Y;
    CLOSE(110);
!調用計算各段表達式系數(shù)的子程序。
    CALL SPL(L,X,Y,B,C,D);
    OPEN(120,FILE='22.TXT');
    OPEN(130,FILE='222.TXT');
!輸出各段表達式的系數(shù)，插值點的深度及其到自激自收點的距離。
    WRITE (120,150);
150 FORMAT(10X,1HI,15X,4HB(I),15X,4HC(I),15X,4HD(I));
    DO 180 I=1,L;
    WRITE(120,170)I,B(I),C(I),D(I);
170 FORMAT(8X,I3,10X,F10.5,2(9X,F10.5));
180 CONTINUE;  
    
    DO 190 P=0,5000;
	Q=P*0.1;
    S(P+1)=SEVA(L,Q,X,Y,B,C,D);
    E(P+1)=SQRT(S(P+1)**2+(M-Q)**2);
	
    WRITE(120,200)Q,S(P+1),Q,E(P+1);
!重新將插值點到自激自收點處的距離值輸出到文本文檔222中。
    WRITE(130,*)E(P+1);
200 FORMAT(10X,2HS(,F6.1,2H)=,F10.5,5X,2HE(,F6.1,2H)=,F10.5);
190 CONTINUE;
    Z(1)=0.;
	WRITE(130,240)Z(1);
240	FORMAT(3X,F8.6);
    DO 1000 I=1,5000;
	Z(I+1)=SQRT(0.01+(S(I+1)-S(I))**2);
	WRITE(130,*)Z(I+1);
1000 CONTINUE;	 
    CLOSE(120);
    CLOSE(130);
    OPEN(150,FILE='222.TXT');
    READ(150,*)F;
!讀取距離值后，選取出其中最小值R，這就是自激自收情況下地震波的最小旅行時路徑。
    R=F(1);
    DO 250 I=1,5001;
    IF(F(I).LE.R) R=F(I);
250 CONTINUE;
    T=R+7.875;
!上行中7.875是菲涅耳帶的距離波動范圍，即3150/（4*100）=7.875。
!地震波的頻率設為100Hz，波速度為3150 m/s。
    OPEN(350,FILE='2222.txt');
	W=0;
    DO 300 J=1,5001;
!輸出菲涅耳帶橫坐標的范圍，由于插值點的間距是0.1M，最終菲涅耳帶精度為0.1M。
    IF(F(J).LE.T)WRITE(350,*)(J-1)*0.1;
	IF(F(J).LE.T)W=W+Z(J);	
300 CONTINUE;
!計算并輸出菲涅爾帶的長度。
    WRITE(350,1100)M,W;
1100   FORMAT(7HLENGTH(,I5,2H)=,F10.2);
    CLOSE(350);
    STOP;   
    END;