C
C Program GBOPT to calculate quantities for optimization of the shape of C 2-D gaussian beams for prestack depth migrations C C Version: 5.50 C Date: 2001, June 10 C C Coded by: Karel Zacek C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: zacek@karel.troja.mff.cuni.cz C C....................................................................... C C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings, read by list directed (free C format) input. The strings have thus to be enclosed in C apostrophes. The interactive * external unit may be redirected to C the file containing the data. C (1) 'SEP'/ C 'SEP'...Name of the file with input parameters. C Description of file SEP C If blank, default values of the corresponding data are C considered. C Default: 'SEP'=' '. C C C Input file 'SEP' in the SEP format: C The file has the form of a SEP parameter file. C For the description of the SEP format refer to file C 'sep.for'. C Names of input files with the results of ray tracing from a receiver: C RAYALL='string'... Name of the input file CRT-R with the C quantities stored along rays (see C.R.T.5.5.1). C Default: RAYALL='r01.out' (usually sufficient) C INIALL='string'... Name of the input file CRT-I with the C quantities at the initial points of rays, corresponding C to file RAYALL (see C.R.T.6.1). C Default: INIALL='r01i.out' (usually sufficient) C Names of input formatted files with quantities corresponding to the C common source, interpolated at the nodes of input grid: C NUM='string'...Name of the input file with N1*N2*N3 array C of integer values, containing the number of arrivals at C each gridpoint of the target grid. C Default: NUM='mtt-num.out' C MTT='string'...Name of the input file with the array of real C values, containing for each gridpoint the travel times C of all determined arrivals. C Default: MTT='mtt-tt.out' C AMP='string'...Name of the input file with the array of real C values, containing for each gridpoint the amplitudes. C Default: AMP='mtt-ap.out' C Names of output files: C GBOUT1='string'... Name of the output file containing optimum C value Y0=SQRT(C11/C22) C of the imaginary part of parameter M0. C Default: GBOUT1='y0.out' C GBOUT2='string'... Name of the output file containing optimum C value R0=-C21/C22 of the real part of parameter M0. C Default: GBOUT2='r0.out' C GBOUT3='string'... Name of the output file containing initial C value SIGMA=SQRT(Y0/C22) of a standard deviation for C smoothing. C Default: GBOUT3='sigma.out' C Data specifying the dimensions of the input grid: C O1=real, O2=real, O3=real ... Coordinates of the origin of the C grid. C Default: O1=0., O2=0., O3=0. C D1=real... Grid interval along the X1 axis. C Default: D1=1. C D2=real... Grid interval along the X2 axis. C Default: D2=1. C D3=real... Grid interval along the X3 axis. C Default: D3=1. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C C Input unformatted file RAYALL: C See the description within source code file 'writ.for'. C Description of files CRT-R C C Input unformatted file INIALL: C See the description within source code file 'writ.for'. C Description of files CRT-I C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL AP00 C AP00... File 'ap.for'. EXTERNAL AP28 C AP28... File 'ap.for'. EXTERNAL RARRAY C RARRAY...File 'forms.for'. EXTERNAL RARRAI C RARRAI...File 'forms.for'. EXTERNAL ERROR C ERROR...File 'error.for'. EXTERNAL RSEP1 C RSEP1...File 'sep.for'. EXTERNAL RSEP3T C RSEP3T...File 'sep.for'. EXTERNAL RSEP3I C RSEP3I...File 'sep.for'. EXTERNAL RSEP3R C RSEP3R...File 'sep.for'. C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C pointc.inc C None of the storage locations of the common block are altered. C INCLUDE 'sep.inc' C sep.inc C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL O1,O2,O3,D1,D2,D3,NORM0,NORM1,NORM2,FAK1,FAK2,FAK3,DET REAL TAU,TT,FX,FFX REAL C(10),DC(7),SDC(7),SSDC(3),FN1(7),FFN1(3),T(8) REAL H11,H12,H13,H21,H22,H23,H31,H32,H33,V1,V2,E11,E12,E22 REAL CA,CB,CC11,CC12,CC22 INTEGER LU1,LU2,LU3,LU4,LU5,LU6,LU7,IPUF,INR,IP,IRR,IRRN INTEGER IFN1(7),IFN2(7),IFFN1(3),IFFN2(3),IRAM(MRAM) INTEGER K1,K2,K3,INDX,ISH,IUF,ISM,IGRD,ITAU,N1,N2,N3,N1N2N3 INTEGER ISRF1,NFN1,NFFN1,NFN2,NFFN2 PARAMETER (LU1=21,LU2=22,LU3=23,LU4=24,LU5=25,LU6=26,LU7=27) PARAMETER (LU8=28,LU9=29,LU10=30,LU11=31) CHARACTER*80 FILSEP,RAYALL,INIALL,NUM,MTT,AMP CHARACTER*80 GBOUT1,GBOUT2,GBOUT3,GBOUT4,GBOUT5 EQUIVALENCE (IRAM,RAM) REAL UNDEF,YMIN PARAMETER (UNDEF=-9.999E9,YMIN=1E-9) PARAMETER (INR=26) PARAMETER (IRR=181) C C....................................................................... C C Reading main input data: FILSEP=' ' WRITE(*,'(A)')'+GBOPT: Enter filename : ' READ(*,*) FILSEP IF(FILSEP.EQ.' ') THEN C GBOPT-01 CALL ERROR('GBOPT-01: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF WRITE(*,'(A)')'+GBOPT: Reading, calculating... ' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading filenames of the files with computed rays and C multivalued travel times CALL RSEP3T('RAYALL',RAYALL,'r01.out') CALL RSEP3T('INIALL',INIALL,'r01i.out') CALL RSEP3T('NUM',NUM,'mtt-num.out') CALL RSEP3T('MTT',MTT,'mtt-tt.out') CALL RSEP3T('AMP',AMP,'mtt-ap.out') C C Reading filenames of the output files CALL RSEP3T('GBOUT1',GBOUT1,'y0.out') CALL RSEP3T('GBOUT2',GBOUT2,'r0.out') CALL RSEP3T('GBOUT3',GBOUT3,'sigma.out') C CALL RSEP3T('GBOUT4',GBOUT4,'time.out') C CALL RSEP3T('GBOUT5',GBOUT5,'h21.out') C CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) OPEN(LU2,FILE=RAYALL,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU3,FILE=INIALL,FORM='UNFORMATTED',STATUS='OLD') OPEN(LU4,FILE=NUM,FORM='FORMATTED',STATUS='OLD') OPEN(LU5,FILE=MTT,FORM='FORMATTED',STATUS='OLD') OPEN(LU6,FILE=AMP,FORM='FORMATTED',STATUS='OLD') OPEN(LU7,FILE=GBOUT1,FORM='FORMATTED') OPEN(LU8,FILE=GBOUT2,FORM='FORMATTED') OPEN(LU9,FILE=GBOUT3,FORM='FORMATTED') C OPEN(LU10,FILE=GBOUT4,FORM='FORMATTED') C OPEN(LU11,FILE=GBOUT5,FORM='FORMATTED') N1N2N3=N1*N2*N3 ISH=2*N1N2N3+1 FX=0. FFX=0. IPUF=0 IRRN=1 C CALL RARRAI(LU4,' ','FORMATTED',.FALSE.,0,N1N2N3,IRAM) IRAM(N1N2N3+1)=0 DO 10 IGRD=1,N1N2N3 IRAM(N1N2N3+IGRD+1)=IRAM(N1N2N3+IGRD)+IRAM(IGRD) 10 CONTINUE IF(3*N1N2N3+2*IRAM(ISH)+1.GT.MRAM) THEN C GBOPT-02 CALL ERROR('GBOPT-02: Too small array RAM(MRAM)') C Too small array RAM(MRAM). If possible, increase dimension MRAM C in include file ram.inc. END IF CALL RARRAY(LU5,' ','FORMATTED',.FALSE.,0.,IRAM(ISH), * RAM(ISH+1)) CALL RARRAY(LU6,' ','FORMATTED',.FALSE.,1.,IRAM(ISH), * RAM(ISH+IRAM(ISH)+1)) DO 30 IGRD=1,N1N2N3 IF (IRAM(IGRD).GE.1) THEN NORM1=0. NORM2=0. DO 20 ITAU=1,IRAM(IGRD) NORM0=RAM(ISH+IRAM(ISH)+IRAM(N1N2N3+IGRD)+ITAU)**2 NORM1=NORM1+RAM(ISH+IRAM(ISH)+IRAM(N1N2N3+IGRD)+ITAU)**2 NORM2=NORM2+NORM0*RAM(ISH+IRAM(N1N2N3+IGRD)+ITAU) 20 CONTINUE IF (NORM1.EQ.0.) THEN RAM(ISH+2*IRAM(ISH)+IGRD)=-1. ELSE RAM(ISH+2*IRAM(ISH)+IGRD)=NORM2/NORM1 END IF ELSE RAM(ISH+2*IRAM(ISH)+IGRD)=-1. END IF 30 CONTINUE C Loop for the points of rays 40 CONTINUE C Reading the results of the complete ray tracing CALL AP00(0,LU2,LU3) IF(IWAVE.LT.1)THEN C End of rays IRRN=IRRN-1 IF (IPUF.LT.INR) THEN IPUF=INR-IPUF DO 41 IP=1,IPUF C WRITE(LU7,'(1(G16.6,1X))') UNDEF C WRITE(LU8,'(1(G16.6,1X))') UNDEF C WRITE(LU9,'(1(G16.6,1X))') UNDEF C WRITE(LU10,'(1(G16.6,1X))') UNDEF C WRITE(LU11,'(1(G16.6,1X))') UNDEF WRITE(LU7,'(I7,A)') 1,'*' WRITE(LU8,'(I7,A)') 1,'*' WRITE(LU9,'(I7,A)') 1,'*' 41 CONTINUE ELSE C GBOPT-03 CALL ERROR('GBOPT-03: Too low parameter INR') END IF IF (IRRN.LT.IRR) THEN IRRN=INR*(IRR-IRRN) DO 42 IP=1,IRRN C WRITE(LU7,'(1(G16.6,1X))') UNDEF C WRITE(LU8,'(1(G16.6,1X))') UNDEF C WRITE(LU9,'(1(G16.6,1X))') UNDEF C WRITE(LU10,'(1(G16.6,1X))') UNDEF C WRITE(LU11,'(1(G16.6,1X))') UNDEF WRITE(LU7,'(I7,A)') 1,'*' WRITE(LU8,'(I7,A)') 1,'*' WRITE(LU9,'(I7,A)') 1,'*' 42 CONTINUE END IF GO TO 100 END IF IF (IPT.LE.1) THEN C New ray IF (IRAY.GT.IRRN) THEN IRRN=INR*(IRAY-IRRN) DO 45 IP=1,IRRN C WRITE(LU7,'(1(G16.6,1X))') UNDEF C WRITE(LU8,'(1(G16.6,1X))') UNDEF C WRITE(LU9,'(1(G16.6,1X))') UNDEF C WRITE(LU10,'(1(G16.6,1X))') UNDEF C WRITE(LU11,'(1(G16.6,1X))') UNDEF WRITE(LU7,'(I7,A)') 1,'*' WRITE(LU8,'(I7,A)') 1,'*' WRITE(LU9,'(I7,A)') 1,'*' 45 CONTINUE END IF IRRN=IRAY+1 IF (IRAY.GT.1) THEN IF (IPUF.LT.INR) THEN IPUF=INR-IPUF DO 48 IP=1,IPUF C WRITE(LU7,'(1(G16.6,1X))') UNDEF C WRITE(LU8,'(1(G16.6,1X))') UNDEF C WRITE(LU9,'(1(G16.6,1X))') UNDEF C WRITE(LU10,'(1(G16.6,1X))') UNDEF C WRITE(LU11,'(1(G16.6,1X))') UNDEF WRITE(LU7,'(I7,A)') 1,'*' WRITE(LU8,'(I7,A)') 1,'*' WRITE(LU9,'(I7,A)') 1,'*' 48 CONTINUE ELSE C GBOPT-04 CALL ERROR('GBOPT-04: Too low parameter INR') END IF END IF IPUF=0 NFN1=7 NFFN1=3 NFN2=7 NFFN2=3 DO 50 ISM=1,7 IFN1(ISM)=ISM FN1(ISM)=0. IFN2(ISM)=ISM 50 CONTINUE DO 55 ISM=1,3 IFFN1(ISM)=ISM FFN1(ISM)=0. IFFN2(ISM)=ISM 55 CONTINUE END IF C Numerical quadrature: DC(1)=Y(12)*Y(20)+Y(13)*Y(21) DC(2)=Y(16)*Y(20)+Y(17)*Y(21) DC(3)=Y(20)*Y(20)+Y(21)*Y(21) DC(4)=Y(12)*Y(24)+Y(13)*Y(25) DC(5)=Y(16)*Y(24)+Y(17)*Y(25) DC(6)=Y(20)*Y(24)+Y(21)*Y(25) DC(7)=Y(24)*Y(24)+Y(25)*Y(25) CALL AP28(7,SDC,1,1,0.01,FX,ISRF1,NFN1,IFN1,FN1,NFN2,IFN2,DC) C(4)=SDC(1) C(5)=SDC(2) C(6)=SDC(3) C(7)=SDC(4) C(8)=SDC(5) C(9)=SDC(6) C(10)=SDC(7) DC(4)=Y(12) DC(5)=Y(13) DC(6)=Y(16) DC(7)=Y(17) DET=C(6)*C(10)-C(9)*C(9) IF (DET.NE.0.) THEN C DC(4)=Y(12)-((Y(20)*C(10)-Y(24)*C(9))*C(4)+(Y(24)*C(6)- C * Y(20)*C(9))*C(7))/DET C DC(5)=Y(13)-((Y(21)*C(10)-Y(25)*C(9))*C(4)+(Y(25)*C(6)- C * Y(21)*C(9))*C(7))/DET C DC(6)=Y(16)-((Y(20)*C(10)-Y(24)*C(9))*C(5)+(Y(24)*C(6)- C * Y(20)*C(9))*C(8))/DET C DC(7)=Y(17)-((Y(21)*C(10)-Y(25)*C(9))*C(5)+(Y(25)*C(6)- C * Y(21)*C(9))*C(8))/DET C DC(4)=(Y(20)*C(10)-Y(24)*C(9))*C(4) DC(4)=DC(4)+(Y(24)*C(6)-Y(20)*C(9))*C(7) DC(4)=DC(4)/DET DC(4)=Y(12)-DC(4) C DC(5)=(Y(21)*C(10)-Y(25)*C(9))*C(4) DC(5)=DC(5)+(Y(25)*C(6)-Y(21)*C(9))*C(7) DC(5)=DC(5)/DET DC(5)=Y(13)-DC(5) C DC(6)=(Y(20)*C(10)-Y(24)*C(9))*C(5) DC(6)=DC(6)+(Y(24)*C(6)-Y(20)*C(9))*C(8) DC(6)=DC(6)/DET DC(6)=Y(16)-DC(6) C DC(7)=(Y(21)*C(10)-Y(25)*C(9))*C(5) DC(7)=DC(7)+(Y(25)*C(6)-Y(21)*C(9))*C(8) DC(7)=DC(7)/DET DC(7)=Y(17)-DC(7) END IF DC(1)=DC(4)**2+DC(5)**2 DC(2)=DC(4)*DC(6)+DC(5)*DC(7) DC(3)=DC(6)**2+DC(7)**2 CALL AP28(3,SSDC,1,1,0.01,FFX,ISRF1, * NFFN1,IFFN1,FFN1,NFFN2,IFFN2,DC) C(1)=SSDC(1) C(2)=SSDC(2) C(3)=SSDC(3) C C Calculation of the mean arrival time at the reciever: IUF=1 K1=INT((Y(3)-O1)/D1) K2=INT((Y(4)-O2)/D2) K3=INT((Y(5)-O3)/D3) C INDX=K3*N1*N2+K2*N1+K1+1 T(1)=RAM(ISH+2*IRAM(ISH)+INDX) C IF (N1.GT.1) THEN INDX=K3*N1*N2+K2*N1+K1+2 ELSE INDX=K3*N1*N2+K2*N1+1 END IF T(2)=RAM(ISH+2*IRAM(ISH)+INDX) C IF (N2.GT.1) THEN INDX=K3*N1*N2+(K2+1)*N1+K1+1 ELSE INDX=K3*N1*N2+K1+1 END IF T(3)=RAM(ISH+2*IRAM(ISH)+INDX) C IF ((N1-1)*(N2-1).GT.0) THEN INDX=K3*N1*N2+(K2+1)*N1+K1+2 ELSE IF (N1.GT.1)THEN INDX=K3*N1*N2+K1+2 ELSE INDX=K3*N1*N2+(K2+1)*N1+1 END IF END IF T(4)=RAM(ISH+2*IRAM(ISH)+INDX) C IF (N3.GT.1) THEN INDX=(K3+1)*N1*N2+K2*N1+K1+1 ELSE INDX=K2*N1+K1+1 END IF T(5)=RAM(ISH+2*IRAM(ISH)+INDX) C IF ((N1-1)*(N3-1).GT.0) THEN INDX=(K3+1)*N1*N2+K2*N1+K1+2 ELSE IF (N1.GT.1)THEN INDX=K2*N1+K1+2 ELSE INDX=(K3+1)*N1*N2+K2*N1+1 END IF END IF T(6)=RAM(ISH+2*IRAM(ISH)+INDX) C IF ((N2-1)*(N3-1).GT.0) THEN INDX=(K3+1)*N1*N2+(K2+1)*N1+K1+1 ELSE IF (N2.GT.1)THEN INDX=(K2+1)*N1+K1+1 ELSE INDX=(K3+1)*N1*N2+K1+1 END IF END IF T(7)=RAM(ISH+2*IRAM(ISH)+INDX) C IF ((N1-1)*(N2-1)*(N3-1).GT.0) THEN INDX=(K3+1)*N1*N2+(K2+1)*N1+K1+2 ELSE IF (N1.GT.1) THEN IF (N2.GT.1)THEN INDX=(K2+1)*N1+K1+2 ELSE INDX=(K3+1)*N1*N2+K1+2 END IF ELSE INDX=(K3+1)*N1*N2+(K2+1)*N1+1 END IF END IF T(8)=RAM(ISH+2*IRAM(ISH)+INDX) C DO 70 ISM=1,8 IF (SIGN(1,T(ISM)).EQ.-1) THEN IUF=0 END IF 70 CONTINUE DO 75 ISM=3,5 IF (SIGN(1,Y(ISM)).EQ.-1) THEN IUF=0 END IF 75 CONTINUE IF (SIGN(1,(N1-1)*((N1-2)*D1-Y(3))).EQ.-1) THEN IUF=0 END IF IF (SIGN(1,(N2-1)*((N2-2)*D2-Y(4))).EQ.-1) THEN IUF=0 END IF IF (SIGN(1,(N3-1)*((N3-2)*D3-Y(5))).EQ.-1) THEN IUF=0 END IF IF (IUF.EQ.1) THEN FAK1=(Y(3)-O1)/D1-K1 FAK2=(Y(4)-O2)/D2-K2 FAK3=(Y(5)-O3)/D3-K3 TAU=T(1)+(T(2)-T(1))*FAK1+(T(3)-T(1)+ * (T(1)-T(2)-T(3)+T(4))*FAK1)*FAK2 TT= T(5)+(T(6)-T(5))*FAK1+(T(7)-T(5)+ * (T(5)-T(6)-T(7)+T(8))*FAK1)*FAK2 TAU=TAU+(TT-TAU)*FAK3 TAU=Y(1)+TAU C C Writing the output data: IPUF=IPUF+1 H11=YI(9) H21=YI(10) H31=YI(11) H13=YLI(1)*YI(6) H23=YLI(1)*YI(7) H33=YLI(1)*YI(8) H12=H23*H31-H33*H21 H22=H33*H11-H13*H31 H32=H13*H21-H23*H11 V1=YLI(4)/YLI(1)**2 V2=YLI(5)/YLI(1)**2 C E11=-H13*V1-H13*V1+H13*H13*(H13*V1+H23*V2) C E12=-H13*V2-H23*V1+H13*H23*(H13*V1+H23*V2) E22=-2*H23*V2+H23*H23*(H13*V1+H23*V2) C C CA=C11surf, CB=C22surf, CC=C21surf C H21=1. IF (H21.LT.0.1) THEN H21=0.1 END IF C11=C(1)*H21**2 C22=C(6)/H21**2 C CC11=C(4)-C(6)*E22/H21**2 C CC12=-C(6)*E22/H21**2 C21=C(4)-C(6)*E22/H21**2 C WRITE(LU7,'(1(G16.6,1X))') CA/Y(1) C WRITE(LU8,'(1(G16.6,1X))') CB/Y(1)**2 C WRITE(LU9,'(1(G16.6,1X))') CC11/Y(1) C WRITE(LU10,'(1(G16.6,1X))') CC12/Y(1) C WRITE(LU11,'(1(G16.6,1X))') CC22/Y(1) C C IF (Y0.LT.YMIN) THEN C Y0=UNDEF C END IF C IF ((C11/C22).LT.0.) THEN C Y0=UNDEF C ELSE Y0=SQRT(C11/C22) C END IF C R011=-CC11/CB C R012=-CC12/CB R0=-C21/C22 SIGMA=SQRT(Y0/C22) WRITE(LU7,'(1(G16.6,1X))') Y0 WRITE(LU8,'(1(G16.6,1X))') R0 WRITE(LU9,'(1(G16.6,1X))') SIGMA C WRITE(LU10,'(1(G16.6,1X))') R012 C WRITE(LU11,'(1(G16.6,1X))') H21 ELSE IPUF=IPUF+1 C WRITE(LU7,'(1(G16.6,1X))') UNDEF C WRITE(LU8,'(1(G16.6,1X))') UNDEF C WRITE(LU9,'(1(G16.6,1X))') UNDEF C WRITE(LU10,'(1(G16.6,1X))') UNDEF C WRITE(LU11,'(1(G16.6,1X))') UNDEF WRITE(LU7,'(I7,A)') 1,'*' WRITE(LU8,'(I7,A)') 1,'*' WRITE(LU9,'(I7,A)') 1,'*' END IF C GO TO 40 C 100 CONTINUE C DO 120 ISM=1,IRR C TIME=0. C IF (ISM.LE.11) THEN C H21=SIN(10*0.034907/2) C END IF C IF (ISM.GT.11) THEN C H21=SIN(ISM*0.017165) C END IF C IF (ISM.GE.171) THEN C H21=SIN(172*0.034907/2) C END IF C DO 110 IP=1,INR C TIME=TIME+0.1 C WRITE(LU10,'(1(G16.6,1X))') TIME C WRITE(LU11,'(1(G16.6,1X))') H21**2 C 110 CONTINUE C 120 CONTINUE C CLOSE(LU2) CLOSE(LU3) CLOSE(LU4) CLOSE(LU5) CLOSE(LU6) CLOSE(LU7) CLOSE(LU8) CLOSE(LU9) C CLOSE(LU10) C CLOSE(LU11) STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'ap.for' C ap.for INCLUDE 'forms.for' C forms.forC C C======================================================================= C