Commit f0488b8c authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

reintegrated branch gresy

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.
gresy now supports selection of sampling parameters and cloning of
sampling in a master seismogram files
there appears a problem with the mergeinfo property of the directory
it indicates a merge from trunk - ?


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 3881
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 8c8d3f6c
......@@ -101,8 +101,8 @@ gresynoise: gresynoise.o nyquist_check.o
-ltf -lgsl -lgslcblas -lsff -L$(LOCLIBDIR)
gresy: gresy.o nyquist_check.o
$(FC) -o $@ $^ -lgrrefsub -lrefread \
-ltf -lsff -L$(LOCLIBDIR)
$(FC) -o $@ $^ -lgrrefsub -lrefread -lsffu \
-ltf -lsff -ltime -L$(LOCLIBDIR)
rhesyg: rhesyg.o
$(FC) -o rhesyg rhesyg.o -lgrrefsub -lrheology -lrefread \
......@@ -114,7 +114,7 @@ afehsyg: syg.o
-ltf -L$(LOCLIBDIR)
gresyx: gresy.o nyquist_check.o
$(FC) -o $@ $^ -lgrrefsub -lrefread -ltf \
$(FC) -o $@ $^ -lgrrefsub -lrefread -lsffu -ltf \
-lfapidxx -ldatrwxx -lsffxx -lgsexx -ltime++ -laff -L$(LOCLIBDIR)
# ----- END OF Makefile -----
......@@ -50,41 +50,49 @@ c 02/03/10 V1.14 corrected (hopefully) line source integration
c 12/01/2011 V1.14b - added definition of Fourier transformation
c online help text
c 18/01/2011 V1.15 implement libfapidxx interface
c 10/04/2011 V1.16 program supports selection of sampling parameters
c 26/04/2011 master file mode is tested
c
c==============================================================================
c
program gresy
character*79 version
parameter(version='GRESY V1.15 GREens function SYnthetics')
parameter(version='GRESY V1.16 GREens function SYnthetics')
c dimensions
integer maxtr, maxsamp, maxom, maxu
parameter(maxu=10000, maxtr=maxu, maxom=4100, maxsamp=maxom*2)
integer ntr, nsamp, nom, nu
integer ntr, nsamp, nom, nu, nsampwrite
c greens function
character*78 greenname
complex green(maxom, maxu)
real om(maxom)
real slo(maxu)
c receiver specifications
real*8 rcvvred, rcvtli, rcvtre
double precision rcvvred, rcvtli, rcvtre
character*80 rcvtext
real*8 phi(maxtr)
double precision phi(maxtr)
c seismogram parameters
integer pown
character*80 seisname,rcvname
real*8 r(maxtr), rmax
double precision r(maxtr), rmax
complex*16 sdata(maxsamp)
real fdata(maxsamp)
integer idata(maxsamp)
equivalence(fdata, idata)
c functions
real sffu_offset, sffu_tfirst
c free block
integer maxfree, nfree
parameter(maxfree=5)
character*80 free(maxfree)
c SFF parameters
real sffversion, tanf
character timestamp*13, code*10
c srce line
character srctype*40, date*20, time*20, cs
real c1, c2, c3
character srctype*40, date*20, time*20, scs
real sc1, sc2, sc3
c info line
character cs
real c1, c2, c3
integer nstack
c seismogram trace to sff file
character*132 wid2line
......@@ -92,13 +100,13 @@ c seismogram trace to sff file
real dt
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=12)
parameter(maxopt=15)
character*3 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c options
logical debug, optlambda, optnew, optresponse, hankel1, hankel2
logical suppress, radial, verbose, linesource
logical suppress, radial, verbose, linesource, usemaster
real lambdalim, tapfrac
character*80 respfile, fileformat
c taper
......@@ -111,26 +119,31 @@ c response file
real omega(maxom)
real domega
c any
integer lu, i, ierr, io, iu
real*8 pi2, arg, du, dom, scal, pi
parameter(lu=20,pi=3.141592653589793115997)
integer mlu,lu, i, ierr, io, iu
double precision pi2, arg, du, dom, scal, pi
parameter(lu=20,mlu=21,pi=3.141592653589793115997)
parameter(pi2=2.d0*pi)
double complex ime
parameter(ime=(0.,1.))
c functions
real*8 tf_dj0, tf_dy0
real*8 tf_dj1, tf_dy1
double precision tf_dj0, tf_dy0
double precision tf_dj1, tf_dy1
c here are the keys to our commandline options
data optid/'-d','-l','-o','-t','-r','-S','-1','-2','-R','-v',
& '-L','-ty'/
data optid/'-D','-l','-o','-t','-r','-S','-1','-2','-R','-v',
& '-L','-ty','-dt', '-N', '-m'/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,2*.TRUE.,6*.FALSE.,
& .true./
data optarg/'-','1.','-','10.','junk',6*'-','sff'/
& 3*.true.,.false./
data optarg/'-','1.','-','10.','junk',6*'-','sff',
& '1.e-18','1000000','-'/
c----------------------------------------------------------------------
c
c read commandline
c
print *,version
print *,'Usage: gresy greenfile seisfile rcvfile [-l lambda] [-o]'
print *,'Usage: gresy greenfile seisfile rcvfile|masterfile'
print *,' [-m] [-l lambda] [-o]'
print *,' [-t frac] [-r respfile] [-S] [-1] [-2] [-R]'
print *,' [-v] [-d] [-L] [-ty f]'
print *,' [-v] [-d] [-L] [-ty f] [-dt dt] [-N n]'
print *,'or gresy -help'
print *,'or gresy -xhelp'
if (iargc().lt.1) stop 'ERROR: no arguments'
......@@ -142,14 +155,29 @@ c
print *,' '
print *,'Calculates synthetic seismograms from greens function.'
print *,' '
print *,'greenfile file containing greens function (in a format'
print *,' produced by ''syg'' or ''greda'')'
print *,'greenfile file containing Fourier-Bessel expansion'
print *,' coefficients of an impulse response'
print *,' (in a format produced by ''syg'' or'
print *,' ''greda'')'
print *,'seisfile output file to contain seismograms'
print *,' '
print *,'rcvfile receiver definition (refmet format)'
print *,'masterfile waveform data file; the traces in the output file'
print *,' ''seisfile'' to be created will have the same'
print *,' spatial and temporal coordinates as the'
print *,' ''masterfile'' '
print *,' '
print *,'rcvfile and masterfile are alternatives to each other'
print *,'see option -m and comment below'
print *,' '
print *,'-v produce verbose output'
print *,'-d produce debug output'
print *,'-D produce debug output'
print *,'-m the third parameter is understood as the name'
print *,' of a ''masterfile'' which must be a valid'
print *,' time series data file.'
print *,'-ty f select seismogram file format f (see below)'
print *,'-dt dt desired sampling interval for output'
print *,'-N n desired number of samples for output'
print *,'-l lambda limits the used slowness range by a minimum'
print *,' wavelength in meters'
print *,' (default: ',optarg(2)(1:4),')'
......@@ -190,6 +218,17 @@ c
print *,' amplitude unit: m**3/s if spectrum represents displacement'
print *,' waveform in m'
print *,' '
print *,'Output trace coordinates and time series sampling can'
print *,'be specified alternatively by'
print *,'a) rcvfile and options -N and -dt'
print *,'or'
print *,'b) masterfile'
print *,'Option -m selection alternative b), such that command '
print *,'line options -N and -dt become ineffective. In any case'
print *,'the sampling of respfile must be appropriate to match'
print *,'the desired properties of output sampling and of the'
print *,'input greenfile.'
print *,' '
print *,'$Id$'
print *,' '
call refmet_rcvinf
......@@ -207,10 +246,10 @@ c
debug=optset(1)
optlambda=optset(2)
if (optlambda) then
read(optarg(2), *) lambdalim
read(optarg(2), *, err=99) lambdalim
endif
optnew=optset(3)
read(optarg(4), *)tapfrac
read(optarg(4), *, err=99)tapfrac
tapfrac=tapfrac/100.
optresponse=optset(5)
respfile=optarg(5)
......@@ -221,8 +260,11 @@ c
verbose=optset(10)
linesource=optset(11)
fileformat=optarg(12)
c
pi2=8.d0*datan(1.d0)
read(optarg(13), *, err=99) dt
if (dt.lt.1.e-20) stop 'ERROR: sampling interval too small'
read(optarg(14), *, err=99) nsampwrite
usemaster=optset(15)
c
call sff_select_output_format(fileformat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting output file format'
......@@ -231,23 +273,40 @@ c----------------------------------------------------------------------
c
c read configuration and greens function
c
call refmet_rrcv(rcvname, rcvtext,
& rcvvred, rcvtli, rcvtre, ntr, maxtr, r, phi, 0.d0,
& 0, 1, debug)
if (.not.usemaster) then
c
c read coordinates from receiver file if no master file is used
c -------------------------------------------------------------
call refmet_rrcv(rcvname, rcvtext,
& rcvvred, rcvtli, rcvtre, ntr, maxtr, r, phi, 0.d0,
& 0, 1, debug)
c
print *,'receiver configuration was taken from:'
print *,' ',rcvtext
if (abs(rcvvred).gt.0.)
& print *,'WARNING: traveltime reduction will be ignored'
if (abs(rcvtli).gt.0.)
& print *,'WARNING: left time shift will be ignored'
if (abs(rcvtre).gt.0.)
& print *,'WARNING: right time shift will be ignored'
if (verbose) then
print *,'calculate seismograms for ',ntr,' receivers at'
do i=1,ntr
print *,r(i),' km ',phi(i),' deg'
enddo
print *,'receiver configuration was taken from:'
print *,' ',rcvtext
if (abs(rcvvred).gt.0.)
& print *,'WARNING: traveltime reduction will be ignored'
if (abs(rcvtli).gt.0.)
& print *,'WARNING: left time shift will be ignored'
if (abs(rcvtre).gt.0.)
& print *,'WARNING: right time shift will be ignored'
if (verbose) then
print *,'calculate seismograms for ',ntr,' receivers at'
do i=1,ntr
print *,r(i),' km ',phi(i),' deg'
enddo
endif
else
c open master file if selected
c ----------------------------
print *,'take spatial and temporal sampling parameters from'
print *,rcvname
call sff_select_input_format(fileformat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting input file format'
call sff_ROpenS(mlu, rcvname, sffversion, timestamp,
& code, srctype, scs, sc1, sc2, sc3,
& date, time, ierr)
if (ierr.ne.0) stop 'ERROR: opening master file'
endif
c
call greenread(greenname, debug,
......@@ -256,7 +315,8 @@ c
c
print *,'some information on green file read:'
print *,' ',nu,' slownesses from ',slo(1),'s/m to ',slo(nu),'s/m'
print *,' ',nom,' frequencies from ',om(1)/pi2,'Hz to ',om(nom)/pi2,'Hz'
print *,' ',nom,' frequencies from ',om(1)/pi2,
& 'Hz to ',om(nom)/pi2,'Hz'
print *,' green frequency interval: ',(om(2)-om(1))/pi2,' Hz'
c
if (suppress) then
......@@ -278,7 +338,8 @@ c
response(io)=(0.,0.)
enddo
c
print *,'read response from file ',respfile(1:index(respfile, ' '))
print *,'read response from file ',
& respfile(1:index(respfile, ' '))
open(lu, file=respfile, err=97)
read(lu, 50, err=96, end=95)
nresp=0
......@@ -295,7 +356,8 @@ c
if (nresp.eq.0) stop 'ERROR: no coefficient read'
print *,' response: read ',nresp,' coefficients'
print *,' response: minimum angular frequency: ',omega(1),' 1/s'
print *,' response: maximum angular frequency: ',omega(nresp),' 1/s'
print *,' response: maximum angular frequency: ',
& omega(nresp),' 1/s'
domega=(omega(nresp)-omega(1))/float(nresp-1)
if ((omega(1)/domega).gt.1.e-10)
& stop 'ERROR: first angular frequency schould be zero'
......@@ -328,34 +390,48 @@ c open seismogram file
c
free(1)=version
free(2)='green file '//greenname(1:index(greenname, ' '))
free(3)='receiver configuration from '
& //rcvname(1:index(rcvname, ' '))
free(4)=rcvtext
nfree=4
srctype='implicit in green file'
cs='C'
c1=0.
c2=0.
c3=0.
date='990418'
time='000000.000'
if (usemaster) then
free(3)='spatial and temporal sampling are taken from '
& //rcvname(1:index(rcvname, ' '))
nfree=3
else
free(3)='receiver configuration from '
& //rcvname(1:index(rcvname, ' '))
free(4)=rcvtext
nfree=4
srctype='implicit in green file'
scs='C'
sc1=0.
sc2=0.
sc3=0.
date='990418'
time='000000.000'
endif
if (optnew) call sff_New(lu, seisname, ierr)
if (ierr.ne.0) stop 'ERROR: deleting seismogram file'
call sff_WOpenFS(lu, seisname, free, nfree, srctype, cs, c1 ,c2 ,c3,
& date, time, ierr)
call sff_WOpenFS(lu, seisname, free, nfree, srctype,
& scs, sc1 ,sc2 ,sc3,
& date, time, ierr)
if (ierr.ne.0) stop 'ERROR: opening seismogram file'
c----------------------------------------------------------------------
c
c build seismic traces
c preparatory calculations
c ------------------------
c
c scale down epicentral distances to meters
do i=1,ntr
r(i)=r(i)*1000.
enddo
if (debug) then
print *,'DEBUG: receiver offsets after scaling to meters:'
print *,(r(i),i=1,ntr)
c scale down epicentral distances to meters if offsets are specified
c by receiver file
if (.not.usemaster) then
do i=1,ntr
r(i)=r(i)*1000.
enddo
if (debug) then
print *,'DEBUG: receiver offsets after scaling to meters:'
print *,(r(i),i=1,ntr)
endif
endif
c check consistency of sampling of Fourier-Bessel coefficients
c ------------------------------------------------------------
c check stepwidth
dom=om(2)-om(1)
do io=2,nom
......@@ -368,27 +444,32 @@ c
if (abs(slo(iu)-slo(iu-1)-du).gt.(0.01*du))
& print *,'WARNING: du varies more than 1%'
enddo
c check trapezoid rule stepsize
rmax=r(1)
do i=1,ntr
rmax=max(rmax,r(i))
call checkslownesssampling(sngl(r(i)),sngl(du),om(nom))
enddo
c
if (om(1).gt.(0.01*dom))
& stop 'ERROR: I assume the first frequency to be 0.Hz!'
c we have got all frequencies now and are going to fit them to
c a power of 2 number
pown=0
1 continue
pown=pown+1
nsamp=2**pown
if (nsamp.lt.maxom) goto 1
c calculate sampling interval
dt=pi2/(nsamp*dom)
print *,'traces will have ',nsamp,' samples at ',dt,
& 's sampling interval'
& stop 'ERROR: program requires the first frequency to be 0.Hz!'
c the following two checks can be done here in advance if the spatial
c sampling is defined by a reciever configuration file
c it will be done on per-trace basis later if the sampling is defined by
c a master file
c
c check trapezoid rule stepsize
if (.not.(usemaster)) then
rmax=r(1)
do i=1,ntr
rmax=max(rmax,r(i))
call checkslownesssampling(sngl(r(i)),sngl(du),om(nom))
enddo
c obtain number of samples and sampling interval
call sampling(dom, nom, dt, maxsamp, nsampwrite, nsamp,
& verbose, debug)
c
print *,'traces will have ',nsampwrite,' samples at ',dt,
& 's sampling interval'
endif
c
c report mode of calculatiomode of calculation
if (hankel1) then
print *,'I will use the first Hankel function'
elseif (hankel2) then
......@@ -398,14 +479,40 @@ c
else
print *,'I will use the Bessel function'
endif
c----------------------------------------------------------------------
c big loop over traces
c --------------------
last=.false.
do i=1,ntr
i=0
do while (.not.last)
i=i+1
c obtain sampling parameters for this trace from master file
if (usemaster) then
nsampwrite=maxsamp
call sff_RTraceI(mlu, tanf, dt,
& wid2line, nsampwrite, fdata,
& idata, code, last,
& cs, c1, c2, c3, nstack, ierr)
r(i)=sffu_offset(cs,c1,c2,c3,scs,sc1,sc2,sc3)
tanf=sffu_tfirst(wid2line, time, date)
call checkslownesssampling(sngl(r(i)),sngl(du),om(nom))
c obtain number of samples and sampling interval
call sampling(dom, nom, dt, maxsamp, nsampwrite, nsamp,
& verbose, debug)
c
print *,'trace will have ',nsampwrite,' samples at ',dt,
& 's sampling interval'
else
if (i.eq.ntr) last=.true.
tanf=0.
endif
print *,'working on trace ',i,' at ',r(i),'m'
if (i.eq.ntr) last=.true.
c calculate bessel integral
do io=1,nom
c set slowness taper
if (debug) print *,'enter frequency loop'
c if (debug) print *,'enter frequency loop'
if (optlambda) then
if (om(io).lt.dom) then
rtap=nu
......@@ -419,12 +526,12 @@ c set slowness taper
rtap=nu
ltap=int(rtap*(1.-tapfrac))
endif
if (i.eq.1) then
if ((verbose).and.(i.eq.1)) then
print *,'u-taper for ',om(io)/pi2,'Hz: ',ltap,'-',rtap,
& ' ',slo(ltap),'s/m-',slo(rtap),'s/m'
endif
sdata(io)=(0.d0,0.d0)
if (debug) print *,'enter slowness loop'
c if (debug) print *,'enter slowness loop'
if (hankel1) then
do iu=1,rtap
arg=max(slo(iu)*om(io)*r(i),1.d-100)
......@@ -484,20 +591,30 @@ c set slowness taper
do io=nom+1,nsamp-nom+1
sdata(io)=(0.d0,0.d0)
enddo
c apply response
c apply response (Fourier coefficients of source wavelet)
c -------------------------------------------------------
if (optresponse) then
do io=1,nom
sdata(io)=sdata(io)*response(io)
enddo
endif
c
c apply time shift
do io=1,nom
sdata(io)=sdata(io)*exp(ime*tanf*om(io))
enddo
c
c make Fourier coefficients complete
do io=1,nom-1
sdata(nsamp-io+1)=conjg(sdata(io+1))
enddo
c
c transform to time domain
if (debug) print *,'go fork'
call tf_dfork(nsamp, sdata, 1.d0)
if (debug) print *,'done fork'
c move trace to writing array and scale to impulse response
c
c copy trace to writing array and scale to impulse response
scal=1.d0/(sqrt(float(nsamp))*dt)
if (debug) print *,'move array'
do io=1,nsamp
......@@ -506,27 +623,37 @@ c move trace to writing array and scale to impulse response
& (abs(imag(sdata(io))).gt.(0.01*abs(real(sdata(io))))))
& print *,'WHERE does this imaginary part come from?'
enddo
c consistency check:
if (nsampwrite.gt.nsamp) stop
& 'ERROR: too many samples requested for output file'
if (.not.(usemaster)) then
c prepare wid2line
call sff_prepwid2(nsamp, 1./dt, 'NSP', 1999, 4, 18, 0, 0, 'NSP ',
& 'NSP ', 'NSP ', 0., -1., -1., -1., -1., wid2line, ierr)
if (ierr.ne.0) stop 'ERROR: preparing WID2 line'
c1=r(i)*cos(pi2*phi(i)/360.)
c2=r(i)*sin(pi2*phi(i)/360.)
c3=0.
cs='C'
nstack=0
if (debug) print *,'go and write',nsamp
call sff_WTraceI(lu, wid2line, nsamp, fdata, idata, last,
call sff_prepwid2(nsampwrite, 1./dt, 'NSP',
& 1999, 4, 18, 0, 0, 'NSP ',
& 'NSP ', 'NSP ', 0., -1., -1., -1., -1., wid2line, ierr)
if (ierr.ne.0) stop 'ERROR: preparing WID2 line'
c1=r(i)*cos(pi2*phi(i)/360.)
c2=r(i)*sin(pi2*phi(i)/360.)
c3=0.
cs='C'
nstack=0
endif
if (debug) print *,'go and write ',nsampwrite,' samples'
call sff_WTraceI(lu, wid2line, nsampwrite, fdata, idata, last,
& cs, c1, c2, c3, nstack, ierr)
if (debug) print *,'did writing'
if (ierr.ne.0) stop 'ERROR: writing trace'
enddo
c -----------------------
c end of loop over traces
c----------------------------------------------------------------------
c final check
call checkslownesssampling(sngl(rmax),sngl(du),om(nom))
c
stop
99 stop 'ERROR: reading command line argument'
98 stop 'ERROR: reading command line argument - to few distances'
97 stop 'ERROR opening response file'
96 stop 'ERROR reading response file'
95 stop 'ERROR reading response file - unexpected end'
......@@ -560,6 +687,7 @@ c
parameter(lu=20)
integer i,j
c
if (debug) print *,'in subroutine greenread'
print *,'read green file ',filename(1:index(filename, ' '))
open(lu, file=filename, form='unformatted', status='old', err=99)
......@@ -603,4 +731,83 @@ c phase slowness is given in s/m
96 stop 'ERROR: closing green file'
end
c----------------------------------------------------------------------
c
c provide sampling data
c
subroutine sampling(dom, nom, dt, Ms, Ns, Ns2, verbose, debug)
c
double precision dom
real dt
integer nom, Ms, Ns, Ns2
logical debug, verbose
c
c input:
c dom: sampling interval for angular frequency values
c nom: number of frequencies
c ms: maximum number of samples, i.e. dimension of time series array
c debug: produce debug output if true