www.pudn.com > Cshot.rar > cshot2.f
* Copyright (c) Colorado School of Mines, 2007.
* All rights reserved.
* Copyright (c) Colorado School of Mines, 2007.
* All rights reserved.
ccccccccc
c CSHOT2 builds traces given the ray data calculated by CSHOT1.
c There are 2 modes of operation:
c (1) Generate shot gathers
c (2) Generate receiver gathers
c For case (1), you specify the record numbers of the shots to build.
c Also, the traces to generate within each shot are given. You can
c create a common offset data set by choosing to build only 1 trace/shot
c and by being careful with the geometry cards in CSHOT1 (ie, you have to make
c sure that the trace in question is the same distance from the shot in all
c records).
c In case (2) you specify a range of shot records to search for the common
c receiver data. You also specify the range of receivers for which gathers
c are to be built.
ccccccccc
cc Variable List:
c AMP Wavelet is caled by this amplitude factor (calculated
c by CSHOT1)
c CAUSTC TRUE if the ray has passed through a caustic (wavelet
c will be phase shifted)
c DATIN Unit number of input data (generated by CSHOT1)
c DATOUT Unit number of output data (traces)
c DELTAT Sample rate used in calculating wavelet
c EOF TRUE if the end of file of the input data has been reached
c FLL,FHL,FLR,FHR Amplitude spectrum of wavelet
c FNYQ Nyquist frequency
c FUDGE Constant scale factor applied to all amplitudes
c HEAD TRUE if this is a head wave event
c I Counter
c IBACK Number of records to backspace in input data file to get
c back to start of current shot (used when paneling the output)
c IDUMMY Dummy integer variable
c INFIL Name of the input data file (created by CSHOT1)
c IRECD Record number
c IRECD1 First (shot) record in sort
c IRECDN Last (shot) record in sort
c IRECV Station number of the receiver
c IRPREV Previous record number
c IROUT Loop variable over receiver sorts
c IROUT1 First receiver sort to output
c IROUTN Last receiver sort to output
c ITR1 First trace to build
c ITRN Last trace to build
c ISOUT Current output record number
c ITRACE Trace number
c ITRMAX Maximum trace to output in this panel of the shot record
c ITRMIN Minimum trace to output in this panel of the shot record
c J,K,L Counters
c JOB Job descriptor (s, r, or o)
c KEY Either the recievr number or the trace number, depending
c on the sort option
c MAXTPT MAximum number of samples/trace
c MAXTR Maximum number of traces/panel - for shot records
c MAXWAV MAximum number of samples in the wavelet (at input sample rate)
c NMOD2 Number samples in wavelet modulo 2
c NPTBY2 NPTWAV / 2
c NPTWAV Number of time samples in wavelet (at input sample rate)
c NREC Number of traces/shot in input data set
c NRSWAV Number of points in the resampled wavelet
c NSHOTS Number of shots in the input data set
c NTPTS Number of samples/trace
c NTRACE Number of traces per output panel
c NWAVMX Maximum number of points in resampled wavelet
c OFFSRT TRUE for common offset (ie, common trace) sort
c OUTFIL Name of the output file containing the traces
c PARIN Unit number of input parameter file PARAM2
c PHASE Phase due to post critical reflections
c PI 3.1415927E+000
c RECSRT TRUE for a receiver sort
c RESAMP Factor by which to resample wavelet (should be even)
c SHTSRT TRUE for shot sort
c SRATE Sample rate in seconds
c TARGET Receiver number or trace number for sort
c TIME Arrival time of wavelet
c TMAX Trace end time
c TRACE(,) Array containing traces
c WAVLEN Length in seconds of the wavelet
c WHEAD() Head wave wavelet
c WIMAG() Imaginary part of (zero phase) wavelet
c WREAL() Real part of (zero phase) wavelet 10240 5104
c XDUMMY Dummy float variable
c ZDUMMY Dummy float variable
INTEGER DATIN, DATOUT, parin, STDERR
PARAMETER ( STDERR = 0,
: DATIN = 12,
: DATOUT = 13,
: parin = 14)
INTEGER MAXTR, MAXTPT, MAXWAV
PARAMETER ( MAXTR = 500,
: MAXTPT =1001,
c : MAXWAV = 513)
: MAXWAV = 257)
integer resamp, nwavmx
c resamp should be an EVEN integer OR 1 for no resampling
parameter ( resamp = 64)
parameter ( nwavmx = maxwav*resamp )
REAL WREAL(2*nwavmx), WIMAG(nwavmx),
: whead(nwavmx), TRACE(MAXTPT,maxtr)
REAL AMP, deltat, FUDGE,
: FLL, FHL, FLR, FHR, FNYQ,
: PHASE, PI, SRATE, TIME, TMAX,
: WAVLEN, xdummy, zdummy
INTEGER I, J, K, L, iback,
: isout, irecd1, irecdn,
: irecd, irprev, irecv, key,
: itrace, itrmin, itrmax, ntrace,
: nmod2, NPTWAV, nptby2, NREC, nrswav,
: nshots, NTPTS, itr1, itrn,
: irout, irout1, iroutn, idummy
CHARACTER INFIL*20, OUTFIL*20, job*1
logical shtsrt, recsrt, caustc, head,
: eof, offsrt
PARAMETER ( PI = 3.141592653589,
: FUDGE = 1.)
OPEN(UNIT=parin,FILE='param2',STATUS='old',ERR=5)
GO TO 10
5 WRITE(STDERR,'(A)') 'Can''t open file : ''param2'''
STOP
10 REWIND parin
c read job description
READ(parin,'(A)',err=15) job
go to 16
15 write(stderr,'(1x,a)')
:'MAIN: Can''t read job option, file ''param2'''
stop
16 continue
c read first and last records to sort on
read(parin,*,err=18) irecd1, irecdn
go to 19
18 write(stderr,'(1x,a,a)')
:'MAIN: Error reading first, last records for sort,',
:' file ''param2'''
stop
19 continue
if(job.eq.'s'.or.job.eq.'S') THEN
shtsrt = .true.
c read first and last traces to build
read(parin,*,err=20) itr1, itrn
go to 21
20 write(stderr,'(1x,a,a)')
: 'MAIN: Error reading first, last traces to build,',
: ' file ''param2'''
stop
21 continue
else if(job.eq.'r'.or.job.eq.'R') THEN
recsrt = .true.
c read receiver numbers for sort
read(parin,*,err=23) irout1, iroutn
go to 24
23 write(stderr,'(1x,a)')
: 'MAIN: Error reading receiver sort numbers, file ''param2'''
stop
24 continue
else
write(stderr,'(/1x,a,/33x,a/)')
: 'MAIN: unrecognized option - use s for shot records,',
: 'r for receiver sort.'
stop
end if
c get corner frequencies of wavelet
READ(parin,*,err=30) FLL,FHL,FLR,FHR
go to 31
30 write(stderr,'(1x,a)')
:'MAIN: Error reading filter frequencies, file ''param2'''
stop
31 continue
IF(FLL.LT.0.) THEN
WRITE(STDERR,'(1x,A)')'MAIN: negative frequency in spectrum'
WRITE(STDERR,'(1x,A)')'Specify positive half of spectrum only'
STOP
END IF
c fll<=fhl<=flr<=fhr
if(fhl.lt.fll) fhl = fll
if(fhr.lt.flr) fhr = flr
if(flr.lt.fhl) then
write(stderr,'(1x,a)')
: 'MAIN: Error in frequency specification of wavelet.'
stop
end if
c length of wavelet
READ(parin,*,err=35) WAVLEN
go to 36
35 write(stderr,'(1x,a)')
:'MAIN: Error reading wavelet length, file ''param2'''
stop
36 continue
c sample rate in seconds
READ(parin,*,err=40) SRATE
go to 41
40 write(stderr,'(1x,a)')
:'MAIN: Error reading sample rate, file ''param2'''
stop
41 continue
c resampling rate for wavelet
deltat = srate / resamp
c Check frequencies less than Nyquist.
FNYQ = 1. / ( 2. * SRATE)
IF(FHR.GT.FNYQ) THEN
WRITE(STDERR,'(1x,A)') 'MAIN : frequency beyond Nyquist.'
STOP
END IF
c number of points in wavelet
nptwav = wavlen / srate + 1
nmod2 = mod(nptwav,2)
if(nmod2.eq.0) then
c make number of points in wavelet an odd number
nptwav = nptwav + 1
end if
if(resamp.eq.1) then
nrswav = nptwav
else
nrswav = nptwav * resamp + 1
end if
nptby2 = ( nptwav - 1 ) / 2
IF(nrswav.GT.nwavmx) THEN
WRITE(STDERR,'(1x,A)') 'MAIN : too many points in wavelet.'
STOP
END IF
c trace length
READ(parin,*,err=50) TMAX
go to 55
50 write(stderr,'(1x,a)')
:'MAIN: Error reading trace length.'
stop
55 continue
NTPTS = 1 + ( TMAX + SRATE / 2. ) / SRATE
IF(NTPTS.GT.MAXTPT) THEN
WRITE(STDERR,'(A)') 'MAIN : too many points / trace.'
STOP
END IF
IF(NPTWAV.GT.NTPTS) THEN
WRITE(STDERR,'(A)') 'MAIN : wavelet length > trace length.'
STOP
END IF
c Read the names of the input and output data files
READ(parin,'(A)',err=70) INFIL
go to 75
70 write(stderr,'(1x,a)')
:'MAIN: Error reading input file name in file ''param2'''
stop
75 continue
READ(parin,'(A)',err=80) OUTFIL
go to 85
80 write(stderr,'(1x,a)')
:'MAIN: Error reading output file name in file ''param2'''
stop
85 continue
OPEN(UNIT=DATIN,FILE=infil,STATUS='old',
:form='UNFORMATTED',ERR=140)
c :ERR=140)
GO TO 145
140 WRITE(STDERR,'(1x,A)') 'Error opening input data file.'
STOP
145 REWIND DATIN
c OPEN(UNIT=DATOUT,FILE=OUTFIL,ERR=150)
OPEN(UNIT=DATOUT,FORM='UNFORMATTED',FILE=OUTFIL,
:ERR=150)
GO TO 155
150 WRITE(STDERR,'(1x,A)') 'Error creating output trace file.'
STOP
155 REWIND DATOUT
c Constructing the wavelet.
CALL WAVLET(FLL,FHL,FLR,FHR,NrsWAV,deltat,
:nwavmx,PI,WREAL,WIMAG,whead,stderr)
if(shtsrt) then
c SHOT GATHERS
c initialize shot counter
isout = irecd1
iback = 0
nrec = itrn - itr1 + 1
c number of traces/record = number of receivers/shot
900 itrmin = itr1
if(nrec.gt.maxtr) then
c will have to panel the output
itrmax = itrmin + maxtr - 1
ntrace = maxtr
else
itrmax = itrn
ntrace = nrec
end if
c zero the trace array for this shot
DO 950 i = 1, ntrace
do 925 j = 1, ntpts
trace(j,i) = 0.
925 continue
950 continue
eof = .false.
1000 read(datin,end=1050) irecd,itrace,irecv,time,
: xdummy,idummy,zdummy,
: amp,phase,caustc,head
if(irecd.lt.isout) go to 1000
if(irecd.gt.isout) go to 1100
iback = iback + 1
if(itrace.lt.itrmin.or.itrace.gt.itrmax) go to 1000
c trace number within this panel
itrace = itrace - itrmin + 1
amp = amp * fudge
c put the wavelet on the trace
call addwav(phase,caustc,nptwav,nptby2,time,ntpts,trace,
: amp,wimag,wreal,whead,pi,itrace,srate,
: deltat,head,resamp,maxtr,maxtpt)
go to 1000
1050 eof = .true.
if(irecd.lt.irecd1) then
c at eof and did not find any records within range
c close file, stop
write(stderr,'(1x,a)') 'No shots within specified range.'
write(stderr,'(1x,a,i4)') 'Max record is ',irecd
close(unit=datout)
stop
end if
1100 continue
c output the shot panel
write(stderr,'(1x,a,i4,1x,a,i4,a,i4)')
: 'Outputting shot ',isout,'traces',itrmin,' - ',itrmax
DO 1300 I = 1, ntrace
write(datout)(trace(j,i),j=1,ntpts)
1300 CONTINUE
if(itrmax.lt.itrn) then
c more traces for this shot
itrmin = itrmin + maxtr
itrmax = itrmax + maxtr
if(itrmax.gt.itrn) itrmax = itrn
ntrace = itrmax - itrmin + 1
c backspace to beginning of data for this shot
do 1305 l = 1, iback + 1
backspace datin
1305 continue
iback = 0
c zero the trace array for the next panel
DO 1320 i = 1, ntrace
do 1310 j = 1, ntpts
trace(j,i) = 0.
1310 continue
1320 continue
go to 1000
end if
if(eof.or.isout.eq.irecdn) then
c done reading events
c close file, stop
close(unit=datout)
stop
else
c next shot
isout = isout + 1
c backspace once to pick up first record of new shot
backspace datin
iback = 0
go to 900
end if
end if
if(recsrt) then
c RECEIVER SORT
c The number of traces is not known in advance when performing a
c receiver sort - it depends on the recording geometry. The
c number of traces will be the number of different shots for which
c the receiver was active (and recorded some events).
c Output one trace at a time
ntrace = 1
do 2000 irout = irout1, iroutn
c initialize the current record number to a non-existent record
irprev = -1
c irprev will become positive when a target event has been read
1500 read(datin,end=1750) irecd,itrace,irecv,time,
: xdummy,idummy,zdummy,
: amp,phase,caustc,head
if(irecd.lt.irecd1) go to 1500
if(irecd.gt.irecdn) go to 1750
if(irecv.ne.irout) go to 1500
if(irecd.ne.irprev) then
c this is a new shot (may be the first)
if(irprev.gt.0) then
c output the trace from the previous shot
c write(stderr,'(1x,a,i4)')'receiver sort-record ',irprev
write(datout)(trace(j,ntrace),j=1,ntpts)
c zero trace for new shot
do 1710 j = 1, ntpts
trace(j,ntrace) = 0.
1710 continue
end if
irprev = irecd
end if
amp = amp*fudge
call addwav(phase,caustc,nptwav,nptby2,time,ntpts,trace,
: amp,wimag,wreal,whead,pi,ntrace,srate,deltat,head,resamp,
: maxtr,maxtpt)
go to 1500
1750 continue
c output last trace for this receiver sort
if(irprev.gt.0) then
c write(stderr,'(1x,a,i4)')'receiver sort-record ',irprev
write(stderr,'(1x,a,i4)')'Done receiver ',irout
write(datout)(trace(j,ntrace),j=1,ntpts)
end if
c done outputting sort for this receiver
c zero trace for new receiver
do 1810 j = 1, ntpts
trace(j,ntrace) = 0.
1810 continue
rewind datin
2000 continue
c close file, stop
close(unit=datout)
stop
end if
STOP
END
c-----------------------------------------------------