C
C Program GPSTEP to calculate the discretization steps C for (a) the decomposition of the wavefield into optimized C 2-D Gaussian packets and (b) migration. C C Version: 6.20 C Date: 2008, June 5 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: C RYMM='string'... Name of the input file containing the C minimum and maximal values of the optimized and smoothed C initial parameters of Gaussian beams R0MIN,R0MAX, C Y0MIN,Y0MAX. C Default: RYMM='ryminmax.out' C Input SEP parameters: C NX=positive integer... Number of gridpoints along the profile C axis X of the decomposed wavefield. C Default: NX=1 C NT=positive integer... Number of gridpoints along the time axis T C of the decomposed wavefield. C Default: NT=1 C DX=real... Grid interval along the profile axis X of C the decomposed wavefield. C Default: DX=1. C DT=real... Grid interval along the time axis T of the decomposed C wavefield. C Default: DT=1. C OT=real, OX=real ... Coordinates of the origin of the grid. C Default: OT=0., OX=0. C ISRC=positive integer... C Default: ISRC=1 C VSURF=real... Velocity determined according to maximum velocity C in the depth of the receivers. C Default: VSURF=2000. C FMIN=real... Minimum frequency of the Gaussian packets. C Default: FMIN=1. C FMAX=real... Maximum frequency of the Gaussian packets. C Default: FMAX=100. C FREF=real... Reference frequencies of the Gaussian packets. C Default: FREF=50. C PARKSI=real... Parameter ksi with values in the interval (0,1). C Default: PARKSI=0.5 C PARNY=real... Parameters ny with values in the interval (0,1). C Default: PARNY=0.5 C KAPPA=real... Parameter kappa has the meaning of the ratio of C the step between two subsequent beams to their C characteristic half-width. C Default: KAPPA=1.253314137 C Name of the output file: C GPSTEP='string'... Name of the output file containing C discretization steps NXR,DXR,OXR,NTR,DTR,OTR,NP,DP,OP, C NOMEGA,DOMEGA,OOMEGA, complex-valued constants C N0, K0, ~K0, Im(N)min (same as Y0MIN). C Default: GPSTEP='gpstep.out' C GBGRD='string'... Name of the output file containing new C grid parameters of the Hamiltonian hypersurface (XR,TR,P) C computed for smoothed and optimized Im(N0) and Re(N0). C C======================================================================= C C External functions and subroutines: EXTERNAL LENGTH,RSEP1,RSEP3T,RSEP3I,RSEP3R,ERROR INTEGER LENGTH C C Filenames and parameters: CHARACTER*80 FILSEP,RYMM,GPSTEP,GBGRD INTEGER LU1,LU2 REAL PI PARAMETER (LU1=1,LU2=2,PI=3.141592654) C C Other variables: INTEGER N,NX,NT,NXR,NTR,NP,NOMEGA,N1,N2,N3,N1NEW,N2NEW,N3NEW,ISRC REAL DX,DT,FMAX,FMIN,FREF,OMEGAMAX,OMEGAREF,VSURF REAL DXR,OXR,DTR,OTR,DP,OP,DOMEGA,OOMEGA,D1,D2,D3,O1,O2,O3 REAL R0MIN,R0MAX,Y0MIN,Y0MAX,PARKSI,PARNY,DSRC,OXRINI REAL ANIMIN,ANIMAX,ANRMIN,ANRMAX,AK0,AK0T,KAPPA COMPLEX AN0 C C....................................................................... C C Reading main input data: FILSEP=' ' WRITE(*,'(A)')'+GPSTEP: Enter filename : ' READ(*,*) FILSEP IF(FILSEP.EQ.' ') THEN C GPSTEP-01 CALL ERROR('GPSTEP-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)')'+GPSTEP: Reading, calculating... ' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading filename of the input file CALL RSEP3T('RYMM',RYMM,'ryminmax.out') C C Reading input parameters CALL RSEP3I('NX',NX,116) CALL RSEP3R('DX',DX,25.0) CALL RSEP3R('OX',OX,0.) CALL RSEP3I('NT',NT,750) CALL RSEP3R('DT',DT,0.004) CALL RSEP3R('OT',OT,0.) C CALL RSEP3I('ISRC',ISRC,1) CALL RSEP3R('DSRC',DSRC,1.) CALL RSEP3R('OXRINI',OXRINI,1.) C Grid parameters of the Hamiltonian hypersurface (NTR,NP,NXR) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('O1',O1,1.) CALL RSEP3R('O2',O2,1.) CALL RSEP3R('O3',O3,1.) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) C CALL RSEP3R('VSURF',VSURF,2000.) CALL RSEP3R('FMAX',FMAX,1.) CALL RSEP3R('FMIN',FMIN,1.) CALL RSEP3R('FREF',FREF,1.) CALL RSEP3R('PARKSI',PARKSI,0.5) CALL RSEP3R('PARNY',PARNY,0.5) CALL RSEP3R('KAPPA',KAPPA,1.253314137) C C Reading filenames of the output files CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out') CALL RSEP3T('GBGRD',GBGRD,'gbgrd.out') C C Read maximum and minimum values of the optimized and smoothed C initial parameters of Gaussian beams R0, Y0 OPEN(LU2,FILE=RYMM,FORM='FORMATTED',STATUS='OLD') READ(LU2,*) R0MIN,R0MAX,Y0MIN,Y0MAX CLOSE(LU2) IF(Y0MIN.LT.0.) THEN C GPSTEP-02 CALL ERROR('GPSTEP-02: Y0 contains negative values') ENDIF C C N=R0+iY0 ANRMIN=R0MIN ANRMAX=R0MAX ANIMIN=Y0MIN ANIMAX=Y0MAX WRITE(*,'(A,1(G16.6,1X))') 'ANRMIN =',ANRMIN WRITE(*,'(A,1(G16.6,1X))') 'ANRMAX =',ANRMAX WRITE(*,'(A,1(G16.6,1X))') 'ANIMIN =',ANIMIN WRITE(*,'(A,1(G16.6,1X))') 'ANIMAX =',ANIMAX C OMEGAMAX=2.0*PI*FMAX OOMEGA=2.0*PI*FMIN OMEGAREF=2.0*PI*FREF IF(OOMEGA.LT.(1.0/(DT*(NT-1)))) THEN OOMEGA= 1.0/(DT*(NT-1)) END IF C C Calculation of parameters N0, K0 and ~K0, equations (33), (34) C and (35), (see Zacek, 2003) AN0=PARNY*2.0*CMPLX(0.0,ANIMIN) AK0=PARKSI*OOMEGA*ANIMIN*VSURF**2 AK0T=AK0*(1.0-PARNY)/PARNY C C Calculation of parameters N0, K0 and ~K0 according to ???? C AN0=CMPLX(ANRMIN,ANIMIN) C AK0=0.5*OOMEGA*ANIMIN*VSURF**2 C AK0T=AK0 C IF(OOMEGA.LT.OMEGAMAX) THEN C Computation of maximum discretization step DOMEGA according to C equation (30), (Zacek, 2003) DOMEGA=KAPPA*SQRT(2.0*AK0*AK0T/(AK0+AK0T)) N=NINT((OMEGAMAX-OOMEGA)/DOMEGA) NOMEGA=N+2 DOMEGA=(OMEGAMAX-OOMEGA)/(N+1) ELSE OOMEGA=OMEGAMAX DOMEGA=1.0 NOMEGA=1 END IF C C Computation of maximum discretization step DXR according to C equation (38), (Zacek, 2003) DXR=KAPPA*SQRT(-2.0*(1.0-PARNY)/OMEGAREF*AIMAG( * (CMPLX(1.0,0.0)-OOMEGA/OMEGAREF* * CMPLX(0.0,ANIMIN)/CMPLX(0.0,ANIMAX)*PARKSI)/ * CMPLX(0.0,ANIMAX))) C C Computation of DXR according to ????? C DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEGAREF*AN0))) C N=NINT((NX-1)*DX/DXR) NXR=N+2 DXR=(NX-1)*DX/(N+1) C C Computation of maximum discretization step DP according to C equation (29), (Zacek, 2003) DP=KAPPA*SQRT(AIMAG(AN0)/OMEGAREF) C N=NINT(2.0/(VSURF*DP)) NP=N+2 DP=2.0/(VSURF*(N+1)) C C Computation of maximum discretization step DTR according to C equation (40), (Zacek, 2003) DTR=KAPPA*SQRT(-2.0*PARNY*AIMAG( * (CMPLX(1.0,0.0)-CMPLX(0.0,ANIMIN)/CMPLX(0.0,ANIMIN)*PARKSI) * /CMPLX(0.0,AK0))) C C Computation of DTR according to ????? C DTR=KAPPA*SQRT(0.5/AK0) C N=NINT((NT-1)*DT/DTR) NTR=N+2 DTR=(NT-1)*DT/(N+1) C OXR=OXRINI+(ISRC-1)*DSRC C C OXR=1587.5-FLOAT(NXR-1)*DXR/2.+ISRC*25.0 C OXR=OX+FLOAT(NX-1)*DX/2.+FLOAT(ISRC-1)*DSRC OP=-1.0/VSURF OTR=0.0 C C Computing new grid parameters of the Hamiltonian hypersurface N1NEW=NTR N2NEW=NP N3NEW=NXR D1=FLOAT(N1-1)/FLOAT(N1NEW-1) D2=FLOAT(N2-1)/FLOAT(N2NEW-1) D3=FLOAT(N3-1)/FLOAT(N3NEW-1) OPEN(LU2,FILE=GBGRD,FORM='FORMATTED') CALL WSEP3I(2,'N1',N1NEW) CALL WSEP3I(2,'N2',N2NEW) CALL WSEP3I(2,'N3',N3NEW) CALL WSEP3R(2,'D1',D1) CALL WSEP3R(2,'D2',D2) CALL WSEP3R(2,'D3',D3) CALL WSEP3R(2,'O1',O1) CALL WSEP3R(2,'O2',O2) CALL WSEP3R(2,'O3',O3) CLOSE(LU2) C C Writing discretization steps to the output file GPSTEP WRITE(*,'(A)')'+GPSTEP: Writing output file.' OPEN(LU2,FILE=GPSTEP,FORM='FORMATTED') WRITE(LU2,'(3A)') '''',GPSTEP(1:LENGTH(GPSTEP)),'''' WRITE(LU2,*) NXR,DXR,OXR,' / (NXR,DXR,OXR)' WRITE(LU2,*) NTR,DTR,OTR,' / (NTR,DTR,OTR)' WRITE(LU2,*) NP, DP, OP ,' / (NP,DP,OP)' WRITE(LU2,*) NOMEGA,DOMEGA,OOMEGA,' / (NOMEGA,DOMEGA,OOMEGA)' C C Writing complex-valued constants to the output file GPSTEP WRITE(LU2,*) REAL(AN0),AIMAG(AN0),' / (Re(N0),Im(N0))' WRITE(LU2,*) AK0,AK0T,ANIMIN,' / (K0,~K0,Im(N)min)' CLOSE(LU2) C C Writing discretization steps to the screen WRITE(*,'(A,I6)') ' NXR =',NXR WRITE(*,'(A,1(G16.6,1X))') ' DXR =',DXR WRITE(*,'(A,1(G16.6,1X))') ' OXR =',OXR WRITE(*,'(A,I6)') ' NTR =',NTR WRITE(*,'(A,1(G16.6,1X))') ' DTR =',DTR WRITE(*,'(A,1(G16.6,1X))') ' OTR =',OTR WRITE(*,'(A,I6)') ' NP =',NP WRITE(*,'(A,1(G16.6,1X))') ' DP =',DP WRITE(*,'(A,1(G16.6,1X))') ' OP =',OP WRITE(*,'(A,I6)') ' NOMEGA=',NOMEGA WRITE(*,'(A,1(G16.6,1X))') ' DOMEGA=',DOMEGA WRITE(*,'(A,1(G16.6,1X))') ' OOMEGA=',OOMEGA C C Writing complex-valued constants to the screen WRITE(*,'(A,1(G16.6,1X))') ' Re(N0) =',REAL(AN0) WRITE(*,'(A,1(G16.6,1X))') ' Im(N0) =',AIMAG(AN0) WRITE(*,'(A,1(G16.6,1X))') ' K0 =',AK0 WRITE(*,'(A,1(G16.6,1X))') ' ~K0 =',AK0T WRITE(*,'(A,1(G16.6,1X))') ' Im(N)min=',ANIMIN C WRITE(*,'(A)') '+GPSTEP: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for C C======================================================================= C