raytracing.f

最优化法射线追踪

      PARAMETER (M1=2000,ML=50,M6=100,M3=120,M7=2*ML+1,M36=M3*M6)
	DIMENSION A(ML),B(ML),C(ML),XR(M3),XS(M6),VS(ML),VP(ML)
	DIMENSION T(M3),XXD(M7),ZZD(M7),X0(M7),Y(M7),R(M7),
     2          RR(M7),AA(M7,M7),AF(M36,M1),WV(13)
	DATA WV/0.0,-0.3,-0.7,-0.1,0.3,0.8,1.5,0.8,0.3,-0.1,
     1       -0.7,-0.3,-0.1/
101	FORMAT(2X,F8.6,1X,F8.5,1X,F8.2,1X,F8.1,1X,F8.1,1X)
102   FORMAT(10(F7.1,1X))
103   FORMAT(10(F7.4,1X))
      OPEN(1,FILE='AAA',STATUS='UNKNOWN')
      OPEN(2,FILE='TIME',STATUS='UNKNOWN')
      OPEN(3,FILE='XZ.FIG',STATUS='UNKNOWN')
	OPEN(4,FILE='A',ACCESS='DIRECT',RECL=4*M1)
	WRITE(*,*)'ENTER NL:'
      READ(1,*) NL
	WRITE(2,*)'NL=',NL
	WRITE(*,*)'ENTER A,B,C,VP,VS:'
	WRITE(2,*)'     A      B      C      VP     VS'
      DO 5 K=1,NL
	READ(1,*)A(K),B(K),C(K),VP(K),VS(K)
5	WRITE(2,101)A(K),B(K),C(K),VP(K),VS(K)
      WRITE(*,*)'ENTER NSP,XS:'
	READ(1,*)NSP
	WRITE(2,*)'NSP=',NSP
	READ(1,*)(XS(I),I=1,NSP)
      WRITE(2,102)(XS(I),I=1,NSP)
	WRITE(*,*)'ENTER NR,XR'
	READ(1,*)NR
	WRITE(2,*)'NR=',NR
      READ(1,*)(XR(J),J=1,NR)
	WRITE(2,102)(XR(J),J=1,NR)
	WRITE(*,*)'ENTER X01,X02,Z01,Z02'
	READ(1,*)X01,X02,Z01,Z02
	WRITE(3,102)X01,X02,Z01,Z02
	WRITE(3,*)0,0
	N1=21
	DO 12 K=1,NL
	   DDX=(X02-X01)/FLOAT(N1-1)
	   DO 14 I=1,N1
            XI=X01+(I-1)*DDX
		  XXD(I)=XI
14	   ZZD(I)=(A(K)*XI+B(K))*XI+C(K) 
         WRITE(3,*)N1,0
	   WRITE(3,102)(XXD(J),ZZD(J),J=1,N1)
12    CONTINUE
      DO 22 KK=1,NSP*NR
	DO 22 II=1,M1
22	AF(KK,II)=0.0
      DO 10 I=1,NSP
	   XSK=XS(I)
	   DO 15 K=1,NL
	      N=K
		  NN=2*N
		  DO 20 L=1,NR
		     XE=XSK+XR(L)
			 IC=(I-1)*NR+L
			IF(XSK.LT.X01.OR.XE.LT.X01.OR.XSK.GT.X02.OR.XE.GT.X02)THEN 
			   T(L)=-1.0
			   GO TO 20
			 ENDIF
			 XC=(XSK+XE)/2.0
			 DO 25 J=1,N-1
			    ZI=(A(J)*XC+B(J))*XC+C(J)
				ZI1= (A(J+1)*XC+B(J+1))*XC+C(J+1)
				DELX=(ZI1-ZI)/ZI1*(XC-XSK)
				XXD(J+1)=XC-DELX
25			 XXD(NN-J+1)=XC+DELX
	         XXD(N+1)=XC
	         XXD(1)=XSK
	         XXD(NN+1)=XE
	         ISO=1
	         CALL QP(AA,A,B,C,XXD,ZZD,VP,VS,R,RR,Y,NN-1,X0,ISO)
               IF(ISO.EQ.0)THEN
	           T(L)=-1.0
                 GOTO 20
	         ENDIF
	         TTT=0.0
	         DO 30 J=1,N
	            XI=XXD(J+1)-XXD(J)
              	ZI=ZZD(J+1)-ZZD(J)
	            RI=SQRT(XI*XI+ZI*ZI)
	            ID=NN-J+1
	            XID=XXD(ID+1)-XXD(ID)
	            ZID=ZZD(ID+1)-ZZD(ID)
                  RID=SQRT(XID*XID+ZID*ZID)
	            TT=RI/VP(J)+RID/VS(J)
	            TTT=TTT+TT
30             CONTINUE
               T(L)=TTT
	         WRITE(3,*)NN+1,1
	         WRITE(3,102)(XXD(J),ZZD(J),J=1,NN+1)
	         IT=NINT(TTT/0.004)
	         IF(IT.GT.M1)THEN
	           WRITE(*,*)'IT>M1'
	           STOP
	         ENDIF
	         DO 35 JJ=1,13
35	         AF(IC,IT+JJ-1)=AF(IC,IT+JJ-1)+WV(JJ)
20          CONTINUE
            WRITE(*,*)'SHOT NO.=',I,'LAYER NO.=',K
	      WRITE(*,103)(T(J),J=1,NR)
15       CONTINUE
10    CONTINUE
      DO 50 L=1,NSP*NR
50    WRITE(4,REC=L)(AF(L,JJ),JJ=1,M1)
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      CLOSE(4)
	STOP
	END
c
      SUBROUTINE QP(AA,A,B,C,X,Z,VP,VS,R,RR,Y,NN,X0,ISO)
	PARAMETER (ML=50,M7=2*ML+1)
	DIMENSION AA(M7,M7),X(M7),Z(M7),Y(M7),X0(M7),A(ML),B(ML),C(ML),
     +          VP(ML),R(M7),RR(M7),VS(ML)
      HH=0.5
	W=0.1
	ESP=1.0E-15
	ESS=1.0E8
	NNS=150
	H=HH
	L=0
1	CALL GS(A,B,C,X,Z,NN)
      CALL FF(A,B,X,Z,VP,VS,R,RR,Y,NN)
	DO 10 I=1,NN
	   IF(ABS(Y(I)).GT.ESP)GOTO 15
10	CONTINUE 
      KEY=1
	RETURN
15    DO 30 J=1,NN
         X(J+1)=X(J+1)+H
	   CALL GS(A,B,C,X,Z,NN)
	   CALL FF(A,B,X,Z,VP,VS,R,RR,Y,NN)
	   DO 20 K=1,NN
20	   AA(K,J)=Y(K)
30    X(J+1)=X(J+1)-H
      CALL GS(A,B,C,X,Z,NN)
	CALL FF(A,B,X,Z,VP,VS,R,RR,Y,NN)
	CALL GOMP(AA,Y,NN,KEY)
	IF(KEY.EQ.0) GO TO 60
	P=1
	DO 40 I=1,NN
40	P=P-Y(I)
      IF(KEY.EQ.0) GOTO 70
	DO 50 I=1,NN
50    X0(I+1)=X(I+1)-(Y(I)*H)/P
      L=L+1
	IF(L.EQ.1)THEN
	  DO 51 I=1,NN
51	  X0(I+1)=X(I+1)
        H=H*W
	  GOTO 1
	  ELSE
	  F=0.0
	  DO 52 I=1,NN
52	  F=F+ABS(X0(I+1)-X(I+1))**2
	  IF(F.GT.ESS)THEN
	    ISO=0
	    RETURN
	  ENDIF
	  IF(F.GT.ESS)THEN
	    IF(L.GT.NNS)THEN
	       ISO=0
	       RETURN
	    ENDIF
	    DO 53 I=1,NN
53	    X0(I+1)=X(I+1)
          H=W*H
	    GOTO 1
	    ELSE
	    KEY=1
	    RETURN
	  ENDIF
	ENDIF
60	KEY=-1
	RETURN
70	KEY=-2
	RETURN
	END
C	
      SUBROUTINE GOMP(AA,Y,NN,KEY)
	PARAMETER(ML=50,M7=2*ML+1)
	DIMENSION AA(M7,M7),Y(M7)
	DO 18 M=1,NN
	   P=0.0
	   DO 11 I=M,NN
	      DO 11 J=1,NN
	         IF(ABS(AA(I,J)).LE.P)GOTO 11
	         P=ABS(AA(I,J))
	         K=I
	         L=J
11       CONTINUE
         KEY=0
	   IF(P.LE.0)RETURN
	   DO 12 J=1,NN
	      Q=AA(K,J)
	      AA(K,J)=AA(M,J)
12       AA(M,J)=Q
         Q=Y(K)
	   Y(K)=Y(M)
	   Y(M)=Q
	   DO 16 I=1,NN
	      Q=-AA(I,L)/AA(M,L)
	      IF(I.EQ.M) GOTO 16
	      DO 14 J=1,NN
14          AA(I,J)=AA(I,J)+AA(M,J)*Q
            Y(I)=Y(I)+Y(M)*Q
16	   CONTINUE
         DO 17 J=1,NN
17	   AA(M,J)=AA(M,J)/Q
18	Y(M)=Y(M)/Q
      DO 21 J=1,NN
	   DO 21 I=1,NN
	      IF(AA(I,J).LE.0.5) GOTO 21
	      DO 19 K=1,NN
	         Q=AA(I,K)
	         AA(I,K)=AA(J,K)
19          AA(J,K)=Q
            Q=Y(I)
	      Y(I)=Y(J)
	      Y(J)=Q
21    CONTINUE
      KEY=1
	RETURN
	END

C
      SUBROUTINE FF(A,B,X,Z,VP,VS,R,RR,Y,NN)
	PARAMETER (ML=50,M7=2*ML+1)
	DIMENSION  A(M7),B(M7),X(M7),Z(M7),Y(M7),VP(ML),R(M7),RR(M7),VS(ML)
      DO 20 J=1,NN
	   RR(J)=SQRT((X(J+1)-X(J+2))**2+(Z(J+1)-Z(J+2))**2)
20    RR(J)=SQRT((X(J+1)-X(J))**2+(Z(J+1)-Z(J))**2)
      DO 36 I=1,NN
	   IF(I-(NN+1)/2)30,32,35
30       Y(I)=(X(I+1)-X(I)+(2.0*A(I)*X(I+1)+B(I))*(Z(I+1)-Z(I)))
     1   /(R(I)*VP(I))+(X(I+1)-X(I+2)+(2.0*A(I)*X(I+1)+B(I))
     1   *(Z(I+1)-Z(I+2)))/(RR(I)*VP(I+1))
	   GOTO 36
32       Y(I)=(X(I+1)-X(I)+(2.0*A(I)*X(I+1)+B(I))*(Z(I+1)-Z(I)))
     1   /(R(I)*VP(I))+(X(I+1)-X(I+2)+(2.0*A(I)*X(I+1)+B(I))
     1   *(Z(I+1)-Z(I+2)))/(RR(I)*VS(I))
	   GOTO 36
35       Y(I)=(X(I+1)-X(I)+(2.0*A(I)*X(I+1)+B(I))*(Z(I+1)-Z(I)))
     1   /(R(I)*VP(I))+(X(I+1)-X(I+2)+(2.0*A(I)*X(I+1)+B(I))
     1   *(Z(I+1)-Z(I+2)))/(RR(I)*VS(NN-I+1))
36    CONTINUE
      RETURN
	END
C
      SUBROUTINE GS(A,B,C,X,Z,NN)
	PARAMETER (ML=50,M7=2*ML+1)
	DIMENSION X(M7),Z(M7),A(ML),B(ML),C(ML)
	Z(1)=0
      Z(NN+2)=0
	DO 11 I=1,NN
	   IF(I-(NN+1)/2)10,10,15
10         Z(I+1)=(A(I)*X(I+1)+B(I))*X(I+1)+C(I)
           GOTO 11
15         Z(I+1)=(A(NN+1-I)*X(I+1)+B(NN+1-I))*X(I+1)+C(NN+1-I)
11    CONTINUE
      RETURN
      END