C
C Program GPANAL to calculate the amplitudes of 2-D Gaussian packets 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 SSD='string'... Name of the input file with synthetic C seismograms to be decomposed. C Default: SSD='ssd.out' C FGBR22='string'... Name of the input file with the resampled C smoothed optimized parameters of Gaussian beams Re(N0). C Default: FGBR22=' ' C FGBY22='string'... Name of the input file with the resampled C smoothed optimized parameters of Gaussian beams Im(N0). C Default: FGBY22=' ' C GPSTEP='string'... Name of the input 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=' ' 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 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 GWM=real.. Gaussian beam decay at the halfwidth is EXP(-GWM*GWM/2) C Default: GWM=3. C RELAMP=real... Relative decay of the Gaussian envelope at which C the loop over the points of the seismograms is terminated. C The relative error due to this economizing roughly C corresponds to the value of RELAMP. C Default: RELAMP=0.1 C Names of the output files: C GPAMPR='string'... Name of the output file containing C real part of Gaussian packet amplitudes. C Default: GPAMPR='gpampr.out' C GPAMPI='string'... Name of the output file containing C imaginary part of Gaussian packet amplitudes. C Default: GPAMPI='gpampi.out' C SSC='string'... Name of the output file with composed synthetic C seismograms. C Default: SSC='ssc.out' C Optional parameters specifying the form of the quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C C forms.for. C NUMLIN=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C C forms.for. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C C External functions and subroutines: EXTERNAL RARRAY,RARRAI,ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R C C Filenames and parameters: CHARACTER*80 FILSEP,SSD,SSC,FGBY22,FGBR22,GPAMPR,GPAMPI,GPSTEP INTEGER LU1 REAL PI PARAMETER (LU1=1,PI=3.141592654) C C Other variables: INTEGER IRAM(MRAM) INTEGER NX,NT,NXT,NXTP,NXTPO INTEGER IXR,ITR,IP,IOMEGA,IX,IT,IT0 INTEGER IGRD,IGRD0,IFGBY0,IFGBR0,IFGBY,IFGBR,IFGB INTEGER IGRDFR,IGRDFI,IGRDF0,IGRDF1,IGRDF2,IGRDF3 INTEGER NXR,NTR,NP,NOMEGA INTEGER IXMAX,IXMIN,ITMAX,ITMIN,IXR0,IP0,MXR,MP REAL OX,OT,DX,DT,XX,T REAL DXR,DTR,DP,DOMEGA,OXR,OTR,OP,OOMEGA,XRX,TRT,AX,BT REAL XR,TR,P,OMEGA,FR,FI,W,WR,WI,GWM,DXTPO REAL DXT,PMAX,ANIMIN,Y0MIN,PARKSI,PARNY,KAPPA REAL AK0,AK0T,AMPLR,AMPLI,DXR1,DP1,AN0R,AN0I COMPLEX CN,CNT,CN0,CAMPL,CN44,CN44T,CPT REAL THETAR,THETAI COMPLEX THETA1,THETA4,CNT1,CPT4 REAL THET2R,THET2I,THET3R,THET3I,THET4R,THET4I REAL CNT1R,CNT1I,CPT4R,CPT4I,AUX,COSWW,SINWW REAL RELAMP,TMIN,TMAX,WR0,WI0,WR3,WI3,WR4,WI4 CHARACTER*80 TEXT EQUIVALENCE (IRAM,RAM) C C....................................................................... C C Reading main input data: FILSEP=' ' WRITE(*,'(A)')'+GPANAL: Enter filename : ' READ(*,*) FILSEP IF(FILSEP.EQ.' ') THEN C GPANAL-01 CALL ERROR('GPANAL-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)')'+GPANAL: 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 input files CALL RSEP3T('SSD',SSD,'ssd.out') CALL RSEP3T('FGBY22',FGBY22,'y0m-res.out') CALL RSEP3T('FGBR22',FGBR22,'r0m-res.out') CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.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.1) CALL RSEP3R('OT',OT,0.) CALL RSEP3R('PARKSI',PARKSI,0.5) CALL RSEP3R('PARNY',PARNY,0.5) CALL RSEP3R('KAPPA',KAPPA,1.253314137) CALL RSEP3R('GWM',GWM,3.) CALL RSEP3R('RELAMP',RELAMP,0.1) RELAMP=ALOG(RELAMP) C C Reading filenames of the output files CALL RSEP3T('GPAMPR',GPAMPR,'gpampr.out') CALL RSEP3T('GPAMPI',GPAMPI,'gpampi.out') CALL RSEP3T('SSC',SSC,'ssc.out') CLOSE(LU1) C C Reading discretization steps OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED',STATUS='OLD') READ(LU1,'(A)') TEXT READ(LU1,*) NXR,DXR,OXR READ(LU1,*) NTR,DTR,OTR READ(LU1,*) NP, DP, OP READ(LU1,*) NOMEGA,DOMEGA,OOMEGA C Reading complex-valued constants READ(LU1,*) AN0R,AN0I READ(LU1,*) AK0,AK0T,Y0MIN CLOSE(LU1) C CN0=CMPLX(AN0R,AN0I) ANIMIN=Y0MIN C C Writing discretization steps to the screen WRITE(*,*) ' ' 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)') ' NXR =',NXR WRITE(*,'(A,1(G16.6,1X))') ' DXR =',DXR WRITE(*,'(A,1(G16.6,1X))') ' OXR =',OXR 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))') ' N0R =',REAL(CN0) WRITE(*,'(A,1(G16.6,1X))') ' N0I =',AIMAG(CN0) WRITE(*,'(A,1(G16.6,1X))') ' K0 =',AK0 WRITE(*,'(A,1(G16.6,1X))') ' ~K0 =',AK0T WRITE(*,'(A,1(G16.6,1X))') ' ANIMIN=',ANIMIN C C Writing constants and grid values to the screen WRITE(*,'(A,1(G16.6,1X))') ' PARKSI=',PARKSI WRITE(*,'(A,1(G16.6,1X))') ' PARNY =',PARNY WRITE(*,'(A,1(G16.6,1X))') ' KAPPA =',KAPPA WRITE(*,'(A,I6)') ' NX =',NX WRITE(*,'(A,1(G16.6,1X))') ' DX =',DX WRITE(*,'(A,1(G16.6,1X))') ' OX =',OX WRITE(*,'(A,I6)') ' NT =',NT WRITE(*,'(A,1(G16.6,1X))') ' DT =',DT WRITE(*,'(A,1(G16.6,1X))') ' OT =',OT WRITE(*,*) ' ' C NXTP=NXR*NTR*NP NXT=NX*NT C C Read synthetic seismograms (time section) to be decomposed OPEN(LU1,FILE=SSD,FORM='FORMATTED',STATUS='OLD') CALL RARRAY(LU1,' ','FORMATTED',.FALSE.,0.,NXT,RAM) CLOSE(LU1) C C Read initial parameters of Gaussian beams - imaginary part (Y0) OPEN(LU1,FILE=FGBY22,FORM='FORMATTED',STATUS='OLD') CALL RARRAY(LU1,' ','FORMATTED',.FALSE.,0.,NXTP, * RAM(NXT+1)) CLOSE(LU1) C C Read initial parameters of Gaussian beams - real part (R0) OPEN(LU1,FILE=FGBR22,FORM='FORMATTED',STATUS='OLD') CALL RARRAY(LU1,' ','FORMATTED',.FALSE.,0.,NXTP, * RAM(NXT+NXTP+1)) CLOSE(LU1) C NXTPO=NXR*NTR*NP*NOMEGA DXTPO=DXR*DTR*DP*DOMEGA DXT=DX*DT IFGBY0=NXT IFGBR0=IFGBY0+NXTP IGRDF0=NXT+2*NXTP C IF((IGRDF0+NXTPO).GT.MRAM) THEN C GPANAL-02 CALL ERROR('GPANAL-02: Too small array RAM(MRAM)') C Too small array RAM(MRAM). If possible, increase dimension MRAM C in include file C ram.inc. END IF C C----------------------------------------------------------------------- C Decomposition of the time section into 2-D Gaussian packets C WRITE(*,'(A)')'+GPANAL: Decomposition ' C BT=GWM/SQRT(AK0T) OXR=OX WRITE(*,'(A,F7.2)') ' OXR=',OXR C C Loop over coordinate of intersection of central ray C of Gaussian packet with the profile DO 500 IXR=1,NXR WRITE(*,'(2(A,I3))') '+GPANAL: Decomposition NXR=',IXR,'/',NXR XR=OXR+DXR*(IXR-1) IGRDF3=(IXR-1)*NTR C C Loop over corresponding arrival time of Gaussian packet DO 450 ITR=1,NTR TR=OTR+DTR*(ITR-1) IGRDF2=(IGRDF3+ITR-1)*NP C C Loop over component of slowness vector of Gaussian packet C along profile DO 400 IP=1,NP P=OP+DP*(IP-1) IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA IFGB=((IXR-1)*NP+IP-1)*NTR+ITR IFGBY=IFGBY0+IFGB IFGBR=IFGBR0+IFGB C C Complex valued parameter N determines initial shape C of Gaussian packets CN=CMPLX(RAM(IFGBR),RAM(IFGBY)) C C Computation of ~N according to equation (28), (Zacek, 2006) CNT=CN*CN0/(2.0*CN-CN0) C C Computation of ~p according to equation (29), (Zacek, 2006) CPT=-P*CNT/CN C C Loop over positive circular frequency of Gaussian packet DO 350 IOMEGA=1,NOMEGA C WRITE(*,'(A,I6)') ' IOMEGA =',IOMEGA OMEGA=OOMEGA+DOMEGA*(IOMEGA-1) IGRDFR=IGRDF1+IOMEGA IGRDFI=IGRDFR+NXTPO C C Computation of discretization step DXR1 according to C equation (38), (Zacek, 2003) DXR1=KAPPA*SQRT(-2.0*(1.0-PARNY)/OMEGA*AIMAG( * (CMPLX(1.0,0.0)-OOMEGA/OMEGA * *CMPLX(0.0,ANIMIN)/CN*PARKSI)/CN)) MXR=INT(DXR1/DXR) IF(MXR.LT.1) THEN MXR=1 END IF IXR0=INT(FLOAT(NXR-MXR*INT(FLOAT(NXR-1) * /FLOAT(MXR))-1)/2.0)+1 C C Computation of discretization step DP1 according to C equation (29), (Zacek, 2003) DP1=KAPPA*SQRT(AIMAG(CN0)/OMEGA) MP=INT(DP1/DP) IF(MP.LT.1) THEN MP=1 END IF IP0=INT(FLOAT(NP-MP*INT(FLOAT(NP-1) * /FLOAT(MP))-1)/2.0)+1 C IF((MOD((IXR-IXR0),MXR).EQ.0).AND. * (MOD((IP-IP0),MP).EQ.0)) THEN C FR=0.0 FI=0.0 PMAX=PI/(OMEGA*DX) C IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN C C Computation of N44 according to equation (31), C (Zacek, 2006) CN44=CN*CMPLX(0.0,AK0) * /(OMEGA*CN-CMPLX(0.0,AK0)*P**2) C C Computation of ~N44 according to equation (30), C (Zacek, 2006) CN44T=CNT*CMPLX(0.0,AK0T) * /(OMEGA*CNT-CMPLX(0.0,AK0T)*CPT**2) C C Computation of ~a according to equation (38), C (Zacek, 2006) CAMPL=0.5*CSQRT(CMPLX(0.0,-1.0)*(CN+CNT)) * *CSQRT(CMPLX(0.0,-1.0)*CN44*CN44T * *(OMEGA/CMPLX(0.0,AK0)+OMEGA/CMPLX(0.0,AK0T))) C AMPLR=REAL(CAMPL) AMPLI=AIMAG(CAMPL) C THETA1=0.5*CNT/(CNT-CMPLX(0.0,AK0T)/OMEGA*CPT**2) THETA4=THETA1*CMPLX(0.0,AK0T) THET4R=REAL (THETA4) THET4I=AIMAG(THETA4) CNT1=CNT*THETA1 CNT1R=OMEGA*REAL (CNT1) CNT1I=OMEGA*AIMAG(CNT1) CPT4=CPT*THETA4 CPT4R=REAL (CPT4) CPT4I=AIMAG(CPT4) C C Gaussian beam horizontal halfwidth AX=GWM/SQRT(OMEGA*AIMAG(CNT)) IXMIN=INT((XR-AX)/DX)+1 IF(IXMIN.LT.1) THEN IXMIN=1 END IF IXMAX=INT((XR+AX)/DX)+2 IF(IXMAX.GT.NX) THEN IXMAX=NX END IF ITMIN=INT((TR-BT)/DT)+1 IF(ITMIN.LT.1) THEN ITMIN=1 END IF ITMAX=INT((TR+BT)/DT)+2 IF(ITMAX.GT.NT) THEN ITMAX=NT END IF C C Loop over traces of wavefield (common shot gather) C to be decomposed DO 300 IX=IXMIN,IXMAX XX=OX+DX*(IX-1) XRX=XX-XR IGRD0=(IX-1)*NT AUX=XRX**2 THET2R=AUX*CNT1R-OMEGA*P*XRX THET2I=AUX*CNT1I AUX=2.0*XRX THET3R=AUX*CPT4R-OMEGA THET3I=AUX*CPT4I AUX=0.5*THET3I/THET4I TMAX=AUX*AUX-(RELAMP+THET2I)/THET4I IF(TMAX.GT.0.0) THEN TMAX=SQRT(TMAX) TMIN=(AUX-TMAX) TMAX=(AUX+TMAX) C ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),ITMIN) C ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),ITMAX) ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),1) ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),NT) ELSE ITMIN=1 ITMAX=-1 END IF C C Loop over samples of each trace of wavefield C (common shot gather) to be decomposed IF(ITMIN.LT.ITMAX) THEN IT=(ITMIN+ITMAX)/2 T=OT+DT*(IT-1) TRT=T-TR THETAR=THET2R-(THET3R-THET4R*TRT)*TRT THETAI=THET2I-(THET3I-THET4I*TRT)*TRT W=EXP(-THETAI) COSWW=W*COS(THETAR) SINWW=W*SIN(THETAR) WR=AMPLR*COSWW-AMPLI*SINWW WI=AMPLR*SINWW+AMPLI*COSWW IGRD=IGRD0+IT AUX=RAM(IGRD) FR=FR+WR*AUX FI=FI+WI*AUX IT0=IT WR0=WR WI0=WI W=EXP( (THET3I-THET4I*(2.0*TRT-DT))*DT) AUX= -(THET3R-THET4R*(2.0*TRT-DT))*DT WR3=W*COS(AUX) WI3=W*SIN(AUX) W=EXP(-2.0*THET4I*DT*DT) AUX=2.0*THET4R*DT*DT WR4=W*COS(AUX) WI4=W*SIN(AUX) C DO 251 IGRD=IGRD0+IT0+1,IGRD0+ITMAX AUX=WR3*WR4-WI3*WI4 WI3=WR3*WI4+WI3*WR4 WR3=AUX AUX=WR*WR3-WI*WI3 WI =WR*WI3+WI*WR3 WR =AUX C AUX=RAM(IGRD) FR=FR+WR*AUX FI=FI+WI*AUX 251 CONTINUE WR=WR0 WI=WI0 W=EXP(-(THET3I-THET4I*(2.0*TRT+DT))*DT) AUX= (THET3R-THET4R*(2.0*TRT+DT))*DT WR3=W*COS(AUX) WI3=W*SIN(AUX) C DO 252 IGRD=IGRD0+IT0-1,IGRD0+ITMIN,-1 AUX=WR3*WR4-WI3*WI4 WI3=WR3*WI4+WI3*WR4 WR3=AUX AUX=WR*WR3-WI*WI3 WI =WR*WI3+WI*WR3 WR =AUX C AUX=RAM(IGRD) FR=FR+WR*AUX FI=FI+WI*AUX 252 CONTINUE END IF 300 CONTINUE END IF C C Complex valued amplitudes of GPs RAM(IGRDFR)=MXR*MP*DXT*FR*OMEGA**2/(2.0*PI**3) RAM(IGRDFI)=MXR*MP*DXT*FI*OMEGA**2/(2.0*PI**3) ELSE RAM(IGRDFR)=0.0 RAM(IGRDFI)=0.0 END IF 350 CONTINUE 400 CONTINUE 450 CONTINUE 500 CONTINUE C Writing output grid values IF(GPAMPR.NE.' ') THEN OPEN(LU1,FILE=GPAMPR,FORM='FORMATTED') CALL WARRAY(LU1,' ','FORMATTED',.FALSE.,0.,.FALSE.,0., * NXTPO,RAM(IGRDF0+1)) CLOSE(LU1) END IF IF(GPAMPI.NE.' ') THEN OPEN(LU1,FILE=GPAMPI,FORM='FORMATTED') CALL WARRAY(LU1,' ','FORMATTED',.FALSE.,0.,.FALSE.,0., * NXTPO,RAM(IGRDF0+NXTPO+1)) CLOSE(LU1) END IF C C----------------------------------------------------------------------- C Composition of the synthetic time section WRITE(*,'(A)')'+GPANAL: Composition ' C BT=GWM/SQRT(AK0) C DO 600 IGRD=1,NXT RAM(IGRD)=0.0 600 CONTINUE C C Loop over coordinate of intersection of central ray C of Gaussian packet with the profile DO 1000 IXR=1,NXR WRITE(*,'(2(A,I3))') '+GPANAL: Composition NXR=',IXR,'/',NXR XR=OXR+DXR*(IXR-1) IGRDF3=(IXR-1)*NTR C C Loop over corresponding arrival time of Gaussian packet DO 950 ITR=1,NTR C WRITE(*,'(A,I6)') ' ITR =',ITR TR=OTR+DTR*(ITR-1) IGRDF2=(IGRDF3+ITR-1)*NP C C Loop over component of slowness vector of Gaussian packet C along profile DO 900 IP=1,NP P=OP+DP*(IP-1) IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA IFGB=((IXR-1)*NP+IP-1)*NTR+ITR IFGBY=IFGBY0+IFGB IFGBR=IFGBR0+IFGB C C Complex valued parameter N determines initial shape C of Gaussian packets CN=CMPLX(RAM(IFGBR),RAM(IFGBY)) C C Loop over positive circular frequency of Gaussian packet DO 850 IOMEGA=1,NOMEGA OMEGA=OOMEGA+DOMEGA*(IOMEGA-1) IGRDFR=IGRDF1+IOMEGA IGRDFI=IGRDFR+NXTPO C C IF(((RAM(IGRDFR).GT.900.0).OR.(RAM(IGRDFI).GT.900.0)).OR. C * ((RAM(IGRDFR).LT.-900.0).OR.(RAM(IGRDFI).LT.-900.0)))THEN IF(((RAM(IGRDFR).GT.0.0).OR.(RAM(IGRDFI).GT.0.0)).OR. * ((RAM(IGRDFR).LT.0.0).OR.(RAM(IGRDFI).LT.0.0)))THEN C AMPLR=RAM(IGRDFR) AMPLI=RAM(IGRDFI) PMAX=PI/(OMEGA*DX) IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN THETA1=0.5*CN/(CN-CMPLX(0.0,AK0)/OMEGA*P**2) THETA4=THETA1*CMPLX(0.0,AK0) THET4R=REAL (THETA4) THET4I=AIMAG(THETA4) CNT1=CN*THETA1 CNT1R=OMEGA*REAL (CNT1) CNT1I=OMEGA*AIMAG(CNT1) C C Gaussian beam horizontal halfwidth AX=GWM/SQRT(OMEGA*AIMAG(CN)) IXMIN=INT((XR-AX)/DX)+1 IF(IXMIN.LT.1) THEN IXMIN=1 END IF IXMAX=INT((XR+AX)/DX)+2 IF(IXMAX.GT.NX) THEN IXMAX=NX END IF ITMIN=INT((TR-BT)/DT)+1 IF(ITMIN.LT.1) THEN ITMIN=1 END IF ITMAX=INT((TR+BT)/DT)+2 IF(ITMAX.GT.NT) THEN ITMAX=NT END IF C C Loop over traces of wavefield (common shot gather) C to be decomposed DO 800 IX=IXMIN,IXMAX XX=OX+DX*(IX-1) XRX=XX-XR IGRD0=(IX-1)*NT AUX=XRX**2 THET2R=AUX*CNT1R+OMEGA*P*XRX THET2I=AUX*CNT1I AUX=2.0*P*XRX THET3R=AUX*THET4R+OMEGA THET3I=AUX*THET4I AUX=0.5*THET3I/THET4I TMAX=AUX*AUX-(RELAMP+THET2I)/THET4I IF(TMAX.GT.0.0) THEN TMAX=SQRT(TMAX) TMIN=(AUX-TMAX) TMAX=(AUX+TMAX) C ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),ITMIN) C ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),ITMAX) ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),1) ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),NT) ELSE ITMIN=1 ITMAX=-1 END IF C C Loop over samples of each trace of wavefield IF(ITMIN.LT.ITMAX) THEN IT=(ITMIN+ITMAX)/2 T=OT+DT*(IT-1) TRT=T-TR THETAR=THET2R-(THET3R-THET4R*TRT)*TRT THETAI=THET2I-(THET3I-THET4I*TRT)*TRT W=DXTPO*EXP(-THETAI) COSWW=W*COS(THETAR) SINWW=W*SIN(THETAR) WR=AMPLR*COSWW-AMPLI*SINWW WI=AMPLR*SINWW+AMPLI*COSWW C C Synthetic time section IGRD=IGRD0+IT RAM(IGRD)=RAM(IGRD)+WR IT0=IT WR0=WR WI0=WI W=EXP( (THET3I-THET4I*(2.0*TRT-DT))*DT) AUX= -(THET3R-THET4R*(2.0*TRT-DT))*DT WR3=W*COS(AUX) WI3=W*SIN(AUX) W=EXP(-2.0*THET4I*DT*DT) AUX=2.0*THET4R*DT*DT WR4=W*COS(AUX) WI4=W*SIN(AUX) C DO 751 IGRD=IGRD0+IT0+1,IGRD0+ITMAX AUX=WR3*WR4-WI3*WI4 WI3=WR3*WI4+WI3*WR4 WR3=AUX AUX=WR*WR3-WI*WI3 WI =WR*WI3+WI*WR3 WR =AUX C C Synthetic time section RAM(IGRD)=RAM(IGRD)+WR 751 CONTINUE WR=WR0 WI=WI0 W=EXP(-(THET3I-THET4I*(2.0*TRT+DT))*DT) AUX= (THET3R-THET4R*(2.0*TRT+DT))*DT WR3=W*COS(AUX) WI3=W*SIN(AUX) DO 752 IGRD=IGRD0+IT0-1,IGRD0+ITMIN,-1 AUX=WR3*WR4-WI3*WI4 WI3=WR3*WI4+WI3*WR4 WR3=AUX AUX=WR*WR3-WI*WI3 WI =WR*WI3+WI*WR3 WR =AUX C C Synthetic time section RAM(IGRD)=RAM(IGRD)+WR 752 CONTINUE END IF 800 CONTINUE END IF END IF 850 CONTINUE 900 CONTINUE 950 CONTINUE 1000 CONTINUE C Writing output grid values IF(SSC.NE.' ') THEN OPEN(LU1,FILE=SSC,FORM='FORMATTED') CALL WARRAY(LU1,' ','FORMATTED',.FALSE.,0.,.FALSE.,0., * NXT,RAM) CLOSE(LU1) END IF C WRITE(*,'(A)') '+GPANAL: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for C C======================================================================= C