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

now checks stepsize

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.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 2346
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent ca85196c
......@@ -73,12 +73,12 @@ syg: syg.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
-ltf $(F2CLIB) -L$(LOCLIBDIR)
gresynoise: gresynoise.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
gresynoise: gresynoise.o nyquist_check.o
gcc -o $(LOCBINDIR)/$@ $^ -lgrrefsub -lrefread \
-ltf -lgsl -lgslcblas $(LIBSTUFF) $(F2CLIB) -L$(LOCLIBDIR)
gresy: gresy.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
gresy: gresy.o nyquist_check.o
gcc -o $(LOCBINDIR)/$@ $^ -lgrrefsub -lrefread \
-ltf $(LIBSTUFF) $(F2CLIB) -L$(LOCLIBDIR)
rhesyg: rhesyg.o
......
c this is <gresy.f>
c------------------------------------------------------------------------------
c $Id: gresy.f,v 1.15 2007-05-15 14:38:35 tforb Exp $
c $Id: gresy.f,v 1.16 2007-07-05 15:11:16 tforb Exp $
c
c 26/06/97 by Thomas Forbriger (IfG Stuttgart)
c
......@@ -22,12 +22,13 @@ c 03/04/02 V1.8 support radial component
c 09/01/03 V1.9 be not too strict when checking frequencies fread
c from response file. Check against frequency and not
c against frequency interval.
c 05/07/07 V1.10 check slowness sampling interval
c
c==============================================================================
c
program gresy
character*79 version
parameter(version='GRESY V1.9 GREens function SYnthetics')
parameter(version='GRESY V1.10 GREens function SYnthetics')
c dimensions
integer maxtr, maxsamp, maxom, maxu
parameter(maxu=4500, maxtr=maxu, maxom=4100, maxsamp=maxom*2)
......@@ -44,7 +45,7 @@ c receiver specifications
c seismogram parameters
integer pown
character*80 seisname,rcvname
real*8 r(maxtr)
real*8 r(maxtr), rmax
complex*16 sdata(maxsamp)
real fdata(maxsamp)
integer idata(maxsamp)
......@@ -132,7 +133,7 @@ c
print *,' amplitude unit: m**3/s if spectrum represents displacement'
print *,' waveform in m'
print *,' '
print *,'$Id: gresy.f,v 1.15 2007-05-15 14:38:35 tforb Exp $'
print *,'$Id: gresy.f,v 1.16 2007-07-05 15:11:16 tforb Exp $'
print *,' '
call refmet_rcvinf
stop
......@@ -256,7 +257,8 @@ 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(3)='receiver configuration from '
& //rcvname(1:index(rcvname, ' '))
free(4)=rcvtext
nfree=4
srctype='implicit in green file'
......@@ -291,6 +293,12 @@ 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),sngl(om(nom)))
enddo
c
if (om(1).gt.(0.01*dom))
& stop 'ERROR: I assume the first frequency to be 0.Hz!'
......@@ -424,6 +432,8 @@ c prepare wid2line
if (debug) print *,'did writing'
if (ierr.ne.0) stop 'ERROR: writing trace'
enddo
c final check
call checkslownesssampling(sngl(rmax),sngl(du),sngl(om(nom)))
c
stop
99 stop 'ERROR: reading command line argument'
......
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.15 2007-05-24 15:59:05 tforb Exp $)
c ($Id: gresynoise.f,v 1.16 2007-07-05 15:11:16 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -18,6 +18,7 @@ c noisymized Fourier coefficients
c 24/05/2007 V1.5 write north and east component
c corrected buffer indexing
c V1.6 return to vertical and radial component
c 05/07/2007 V1.7 check slowness sampling interval
c
c ============================================================================
c
......@@ -48,10 +49,10 @@ c
c
character*(*) version
parameter(version=
&'GRESYNOISE V1.6 calculate noise seismograms')
&'GRESYNOISE V1.7 calculate noise seismograms')
character*(*) GRESYNOISE_CVS_ID
parameter(GRESYNOISE_CVS_ID=
&'$Id: gresynoise.f,v 1.15 2007-05-24 15:59:05 tforb Exp $')
&'$Id: gresynoise.f,v 1.16 2007-07-05 15:11:16 tforb Exp $')
c
c dimensions
integer maxtr, maxsetsamp, maxom, maxu, maxset, maxsamp
......@@ -66,7 +67,7 @@ c greens function
real Zom(maxom), Rom(maxom)
real Zslo(maxu), Rslo(maxu)
c receiver specifications
real*8 rcvvred, rcvtli, rcvtre
real*8 rcvvred, rcvtli, rcvtre, rmax
character*80 rcvtext
real*8 phi(maxtr)
real scalefac(maxtr,maxset)
......@@ -352,6 +353,12 @@ c
if (abs(Zslo(iu)-Zslo(iu-1)-du).gt.(0.01*du))
& print *,'WARNING: du varies by 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),sngl(Zom(Znom)))
enddo
c
if (Zom(1).gt.(0.01*dom))
& stop 'ERROR: I assume the first frequency to be 0.Hz!'
......@@ -614,6 +621,8 @@ c prepare wid2line
& cs, c1, c2, c3, nstack, ierr)
if (debug) print *,'written'
if (ierr.ne.0) stop 'ERROR: writing trace'
c final check
call checkslownesssampling(sngl(rmax),sngl(du),sngl(Zom(Znom)))
c
stop
99 stop 'ERROR: reading command line argument'
......
c this is <nyquist_check.f>
c ----------------------------------------------------------------------------
c ($Id: nyquist_check.f,v 1.1 2007-07-05 15:11:17 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
c Check Nyquist criterion for trapezoid rule
c
c REVISIONS and CHANGES
c 05/07/2007 V1.0 Thomas Forbriger
c
c ============================================================================
c
subroutine checkslownesssampling(r,p,o)
c
real r, p, o
c
c r: receiver offset in km
c p: slowness interval in s/km
c o: maximum frequency in 1/s
c
c This function will check whether the slowness sampling is adequate for
c the oscillating Bessel functions in the slowness expansion integral when
c applying a trapezoid rule with slowness interval p.
c
c The subroutine will emit a warning, when the step is larger then pi/2.
c It will abort, when the step size is larger than 1.8*pi.
c
real d,pi2
parameter(pi2=2.*3.14159265358979)
d=o*r*p/pi2
if (d.gt.0.25) then
print *,'WARNING: The trapezoid stepsize is ',d,'*2*pi.'
print *,' Receiver offset: ',r,' km'
print *,' Maximum frequency: ',o/pi2,' Hz'
print *,' Slowness stepsize: ',p,' s/km'
if (d.gt.0.9) then
print *,'!! The stepsize for the argument of the Bessel'
print *,'!! function is larger than 1.8*pi. The'
print *,'!! trapezoid rule approximation is very likely'
print *,'!! to fail. The program aborts for this reason.'
stop 'aborted...'
else
print *,'!! The stepsize for the argument of the Bessel'
print *,'!! function is larger than pi/2. The'
print *,'!! trapezoid rule approximation is very likely'
print *,'!! to be inaccurate!'
print *,''
endif
endif
return
end
c
c ----- END OF nyquist_check.f -----
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment