...
 
Commits (7)
...@@ -3,8 +3,10 @@ this is <CHANGELOG> ...@@ -3,8 +3,10 @@ this is <CHANGELOG>
Recent development in Seitosh (bug fixes, new features, etc) Recent development in Seitosh (bug fixes, new features, etc)
------------------------------------------------------------ ------------------------------------------------------------
11.12.2019: grepg 11.12.2019: grepg, libfapidxx, gresy
- commit 68295a7: provide option to select plot axis range - commit 68295a7: grepg: provide option to select plot axis range
- commit b498110: libfapidxx: support negative time shift
- commit a1e2be8: gresy: provide Fourier domain filters
01.11.2019: libtfxx, randomseries 01.11.2019: libtfxx, randomseries
- commit 27491f4: let RNGgaussian support selection of rng type - commit 27491f4: let RNGgaussian support selection of rng type
......
...@@ -105,7 +105,7 @@ gresynoise: gresynoise.o nyquist_check.o ...@@ -105,7 +105,7 @@ gresynoise: gresynoise.o nyquist_check.o
gresy: gresy.o nyquist_check.o gresy: gresy.o nyquist_check.o
$(FC) -o $@ $^ -lgrrefsub -lrefpar -lsffu \ $(FC) -o $@ $^ -lgrrefsub -lrefpar -lsffu \
-ltf -lsff -ltime -L$(LOCLIBDIR) -ltf -lsff -ltime -lfourier -L$(LOCLIBDIR)
rhesyg: rhesyg.o rhesyg: rhesyg.o
$(FC) -o rhesyg rhesyg.o -lgrrefsub -lrheology -lrefpar \ $(FC) -o rhesyg rhesyg.o -lgrrefsub -lrheology -lrefpar \
...@@ -119,7 +119,8 @@ afehsyg: syg.o ...@@ -119,7 +119,8 @@ afehsyg: syg.o
gresyx: gresy.o nyquist_check.o gresyx: gresy.o nyquist_check.o
$(FC) -o $@ $^ -lgrrefsub -lrefpar -lsffu -ltf \ $(FC) -o $@ $^ -lgrrefsub -lrefpar -lsffu -ltf \
-lfapidxx -ldatrwxx -lsffxx -lgsexx -ltime++ -laff -L$(LOCLIBDIR) -lfapidxx -ldatrwxx -lsffxx -lgsexx -ltime++ -laff \
-lfourier -L$(LOCLIBDIR)
#====================================================================== #======================================================================
# create package # create package
......
...@@ -76,6 +76,7 @@ Dependencies: ...@@ -76,6 +76,7 @@ Dependencies:
Seitosh libraries required to compile the code: Seitosh libraries required to compile the code:
libaff libdatrwxx libfapidxx libgrrefsub libgsexx librefpar librheology libaff libdatrwxx libfapidxx libgrrefsub libgsexx librefpar librheology
libsff libsffu libsffxx libtf libtime libtime++ libwrefsub libtfxx libsff libsffu libsffxx libtf libtime libtime++ libwrefsub libtfxx
libfourier
The libraries can be obtained from where you got the current package from The libraries can be obtained from where you got the current package from
either in a single package containing all the libraries or in individual either in a single package containing all the libraries or in individual
packages. packages.
......
...@@ -2,6 +2,7 @@ c this is <gresy.f> ...@@ -2,6 +2,7 @@ c this is <gresy.f>
c------------------------------------------------------------------------------ c------------------------------------------------------------------------------
c c
c Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart) c Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart)
c Copyright 2019 by Thomas Forbriger (BFO, KIT)
c c
c calculate seismograms from a greens function matrix file c calculate seismograms from a greens function matrix file
c c
...@@ -57,12 +58,13 @@ c 29/11/2011 V1.18 implemented radial component line source ...@@ -57,12 +58,13 @@ c 29/11/2011 V1.18 implemented radial component line source
c both new cases successfully tested against refmet c both new cases successfully tested against refmet
c notice: refmet is not yet tested against a c notice: refmet is not yet tested against a
c reference c reference
c 11/12/2019 V1.19 support Fourier domain filters
c c
c============================================================================== c==============================================================================
c c
program gresy program gresy
character*79 version character*79 version
parameter(version='GRESY V1.18 GREens function SYnthetics') parameter(version='GRESY V1.19 GREens function SYnthetics')
c dimensions c dimensions
integer maxtr, maxsamp, maxom, maxu integer maxtr, maxsamp, maxom, maxu
parameter(maxu=10000, maxtr=maxu, maxom=4100, maxsamp=maxom*2) parameter(maxu=10000, maxtr=maxu, maxom=4100, maxsamp=maxom*2)
...@@ -87,7 +89,7 @@ c functions ...@@ -87,7 +89,7 @@ c functions
real sffu_offset, sffu_tfirst real sffu_offset, sffu_tfirst
c free block c free block
integer maxfree, nfree integer maxfree, nfree
parameter(maxfree=5) parameter(maxfree=15)
character*80 free(maxfree) character*80 free(maxfree)
c SFF parameters c SFF parameters
real sffversion, tanf real sffversion, tanf
...@@ -105,14 +107,18 @@ c seismogram trace to sff file ...@@ -105,14 +107,18 @@ c seismogram trace to sff file
real dt real dt
c commandline c commandline
integer maxopt, lastarg, iargc integer maxopt, lastarg, iargc
parameter(maxopt=15) parameter(maxopt=19)
character*3 optid(maxopt) character*4 optid(maxopt)
character*80 optarg(maxopt) character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt) logical optset(maxopt), opthasarg(maxopt)
c options c options
logical debug, optlambda, optnew, optresponse, hankel1, hankel2 logical debug, optlambda, optnew, optresponse, hankel1, hankel2
logical suppress, radial, verbose, linesource, usemaster logical suppress, radial, verbose, linesource, usemaster
real lambdalim, tapfrac real lambdalim, tapfrac
c low- and high-pass frequency, time shifts (source, origin)
real lpbf,hpbf,tshiftsrc,tshiftorigin
c oredr of low- and high-pass
integer lpbo, hpbo
character*80 respfile, fileformat character*80 respfile, fileformat
c taper c taper
real tf_costap real tf_costap
...@@ -123,6 +129,8 @@ c response file ...@@ -123,6 +129,8 @@ c response file
integer nresp integer nresp
real omega(maxom) real omega(maxom)
real domega real domega
c filter function
double complex fou_eval
c any c any
integer mlu,lu, i, ierr, io, iu integer mlu,lu, i, ierr, io, iu
double precision pi2, arg, du, dom, scal, pi double precision pi2, arg, du, dom, scal, pi
...@@ -135,11 +143,12 @@ c functions ...@@ -135,11 +143,12 @@ c functions
double precision tf_dj1, tf_dy1 double precision tf_dj1, tf_dy1
c here are the keys to our commandline options c here are the keys to our commandline options
data optid/'-D','-l','-o','-t','-r','-S','-1','-2','-R','-v', data optid/'-D','-l','-o','-t','-r','-S','-1','-2','-R','-v',
& '-L','-ty','-dt', '-N', '-m'/ & '-L','-ty','-dt', '-N', '-m', '-lpb', '-hpb',
& '-tss', '-tso'/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,2*.TRUE.,6*.FALSE., data opthasarg/.FALSE.,.TRUE.,.FALSE.,2*.TRUE.,6*.FALSE.,
& 3*.true.,.false./ & 3*.true.,.false.,4*.true./
data optarg/'-','1.','-','10.','junk',6*'-','sff', data optarg/'-','1.','-','10.','junk',6*'-','sff',
& '1.e-18','1000000','-'/ & '1.e-18','1000000','-',2*'1.,0',2*'0.'/
c---------------------------------------------------------------------- c----------------------------------------------------------------------
c c
c read commandline c read commandline
...@@ -149,6 +158,7 @@ c ...@@ -149,6 +158,7 @@ c
print *,' [-m] [-l lambda] [-o]' print *,' [-m] [-l lambda] [-o]'
print *,' [-t frac] [-r respfile] [-S] [-1] [-2] [-R]' print *,' [-t frac] [-r respfile] [-S] [-1] [-2] [-R]'
print *,' [-v] [-d] [-L] [-ty f] [-dt dt] [-N n]' print *,' [-v] [-d] [-L] [-ty f] [-dt dt] [-N n]'
print *,' [-lpb f,o] [-hpb f,o] [-tss t] [-tso t]'
print *,'or gresy -help' print *,'or gresy -help'
print *,'or gresy -xhelp' print *,'or gresy -xhelp'
call tf_seitosh_reference call tf_seitosh_reference
...@@ -207,6 +217,20 @@ c ...@@ -207,6 +217,20 @@ c
print *,'-L simulate seismograms from line source (2D)' print *,'-L simulate seismograms from line source (2D)'
print *,'-R calculate radial component' print *,'-R calculate radial component'
print *,' ' print *,' '
print *,'Fourier domain filter'
print *,' '
print *,'-lpb f,o apply a Butterworth low-pass'
print *,'-hpb f,o apply a Butterworth high-pass'
print *,' f: eigenfrequency of filter (in Hz)'
print *,' o: order of filter'
print *,' '
print *,'-tss t shift source time function by t (in s)'
print *,' t>0: onset is delayed'
print *,' source time in data file is not affected'
print *,'-tso t shift origin time by t (in s)'
print *,' t<0: first sample is before source time'
print *,' header data of traces will be adjusted'
print *,' '
print *,'Definition of the Fourier transformation:' print *,'Definition of the Fourier transformation:'
print *,'The Fourier transformation used in this program and in' print *,'The Fourier transformation used in this program and in'
print *,'related programs (like gremlin, syg, and greda) is' print *,'related programs (like gremlin, syg, and greda) is'
...@@ -279,7 +303,17 @@ c ...@@ -279,7 +303,17 @@ c
if (dt.lt.1.e-20) stop 'ERROR: sampling interval too small' if (dt.lt.1.e-20) stop 'ERROR: sampling interval too small'
read(optarg(14), *, err=99) nsampwrite read(optarg(14), *, err=99) nsampwrite
usemaster=optset(15) usemaster=optset(15)
c read Fourier domain filter parameters
read(optarg(16), *, err=99) lpbf,lpbo
read(optarg(17), *, err=99) hpbf,hpbo
read(optarg(18), *, err=99) tshiftsrc
read(optarg(19), *, err=99) tshiftorigin
c
c check consistency of parameters
if ((lpbf.le.0.).and.(lpbo.gt.0))
& stop 'ERROR: invalid parameters to -lbp'
if ((hpbf.le.0.).and.(hpbo.gt.0))
& stop 'ERROR: invalid parameters to -hbp'
c c
call sff_select_output_format(fileformat, ierr) call sff_select_output_format(fileformat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting output file format' if (ierr.ne.0) stop 'ERROR: selecting output file format'
...@@ -402,6 +436,28 @@ c ...@@ -402,6 +436,28 @@ c
endif endif
c c
c---------------------------------------------------------------------- c----------------------------------------------------------------------
c apply Fourier domain filter
call foufil_clear
call foufil_frequency
if (verbose) then
if (lpbo.gt.0) print *,' apply lpb ',lpbf,',',lpbo
if (hpbo.gt.0) print *,' apply hpb ',hpbf,',',hpbo
print *,' shift source by ',tshiftsrc,'s'
print *,' shift origin by ',tshiftorigin,'s'
endif
if (lpbo.gt.0) call foufil_lpb(dble(lpbf),lpbo)
if (hpbo.gt.0) call foufil_hpb(dble(hpbf),hpbo)
do io=1,nom
response(io)=response(io)*
& cmplx(exp(ime*(tshiftorigin-tshiftsrc)*om(io)))
response(io)=response(io)*
& cmplx(fou_eval(dble(om(io))))
enddo
c
c----------------------------------------------------------------------
c c
c open seismogram file c open seismogram file
c c
...@@ -416,6 +472,26 @@ c ...@@ -416,6 +472,26 @@ c
& //rcvname(1:index(rcvname, ' ')) & //rcvname(1:index(rcvname, ' '))
free(4)=rcvtext free(4)=rcvtext
nfree=4 nfree=4
if (lpbo.gt.0) then
nfree=nfree+1
write (free(nfree), '(a,g8.2,a,i3)')
& 'apply lpb ',lpbf,',',lpbo
endif
if (hpbo.gt.0) then
nfree=nfree+1
write (free(nfree), '(a,g8.2,a,i3)')
& 'apply lpb ',hpbf,',',hpbo
endif
if (tshiftsrc.ne.0) then
nfree=nfree+1
write (free(nfree), '(a,g8.2,a)')
& 'shift source onset by ',tshiftsrc,'s'
endif
if (tshiftorigin.ne.0) then
nfree=nfree+1
write (free(nfree), '(a,g8.2,a)')
& 'shift origin time by ',tshiftorigin,'s'
endif
srctype='implicit in green file' srctype='implicit in green file'
scs='C' scs='C'
sc1=0. sc1=0.
...@@ -424,12 +500,24 @@ c ...@@ -424,12 +500,24 @@ c
date='990418' date='990418'
time='000000.000' time='000000.000'
endif endif
if (verbose) then
print *,'open output file ',seisname(1:index(seisname, ' '))
endif
if (optnew) call sff_New(lu, seisname, ierr) if (optnew) call sff_New(lu, seisname, ierr)
if (ierr.ne.0) stop 'ERROR: deleting seismogram file' if (ierr.ne.0) stop 'ERROR: deleting seismogram file'
call sff_WOpenFS(lu, seisname, free, nfree, srctype, call sff_WOpenFS(lu, seisname, free, nfree, srctype,
& scs, sc1 ,sc2 ,sc3, & scs, sc1 ,sc2 ,sc3,
& date, time, ierr) & date, time, ierr)
if (ierr.ne.0) stop 'ERROR: opening seismogram file' if (ierr.ne.0) stop 'ERROR: opening seismogram file'
if (verbose) then
print *,'information encoded in file header FREE data:'
do i=1,nfree
print *,' ',free(i)(1:76)
enddo
endif
c---------------------------------------------------------------------- c----------------------------------------------------------------------
c c
c preparatory calculations c preparatory calculations
...@@ -650,13 +738,12 @@ c---------------------------------------------------------------------- ...@@ -650,13 +738,12 @@ c----------------------------------------------------------------------
do io=nom+1,nsamp-nom+1 do io=nom+1,nsamp-nom+1
sdata(io)=(0.d0,0.d0) sdata(io)=(0.d0,0.d0)
enddo enddo
c apply response (Fourier coefficients of source wavelet) c apply response
c ------------------------------------------------------- c (Fourier coefficients of source wavelet and Fourier domain filter)
if (optresponse) then c ------------------------------------------------------------------
do io=1,nom do io=1,nom
sdata(io)=sdata(io)*response(io) sdata(io)=sdata(io)*response(io)
enddo enddo
endif
c c
c apply time shift c apply time shift
do io=1,nom do io=1,nom
...@@ -697,6 +784,11 @@ c prepare wid2line ...@@ -697,6 +784,11 @@ c prepare wid2line
cs='C' cs='C'
nstack=0 nstack=0
endif endif
c shift trave, if origin is shifted
if (tshiftorigin.ne.0.) then
call sff_modwid2shift(wid2line, 0., sngl(tshiftorigin))
endif
c
if (debug) print *,'go and write ',nsampwrite,' samples' if (debug) print *,'go and write ',nsampwrite,' samples'
call sff_WTraceI(lu, wid2line, nsampwrite, fdata, idata, last, call sff_WTraceI(lu, wid2line, nsampwrite, fdata, idata, last,
& cs, c1, c2, c3, nstack, ierr) & cs, c1, c2, c3, nstack, ierr)
......
...@@ -32,7 +32,7 @@ ...@@ -32,7 +32,7 @@
* ============================================================================ * ============================================================================
*/ */
#define TF_FAPID_SFF_MODWID2DATE_CC_VERSION \ #define TF_FAPID_SFF_MODWID2DATE_CC_VERSION \
"TF_FAPID_SFF_MODWID2DATE_CC V1.0 " "TF_FAPID_SFF_MODWID2DATE_CC V1.0"
#include <fapidxx/fapidsff.h> #include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h> #include <fapidxx/wid2container.h>
......
...@@ -32,9 +32,7 @@ ...@@ -32,9 +32,7 @@
* ============================================================================ * ============================================================================
*/ */
#define TF_FAPID_SFF_MODWID2SAMPRAT_CC_VERSION \ #define TF_FAPID_SFF_MODWID2SAMPRAT_CC_VERSION \
"TF_FAPID_SFF_MODWID2SAMPRAT_CC V1.0 " "TF_FAPID_SFF_MODWID2SAMPRAT_CC V1.0"
#define TF_FAPID_SFF_MODWID2SAMPRAT_CC_CVSID \
"$Id: $"
#include <fapidxx/fapidsff.h> #include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h> #include <fapidxx/wid2container.h>
......
...@@ -32,9 +32,7 @@ ...@@ -32,9 +32,7 @@
* ============================================================================ * ============================================================================
*/ */
#define TF_FAPID_SFF_MODWID2SAMPS_CC_VERSION \ #define TF_FAPID_SFF_MODWID2SAMPS_CC_VERSION \
"TF_FAPID_SFF_MODWID2SAMPS_CC V1.0 " "TF_FAPID_SFF_MODWID2SAMPS_CC V1.0"
#define TF_FAPID_SFF_MODWID2SAMPS_CC_CVSID \
"$Id: $"
#include <fapidxx/fapidsff.h> #include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h> #include <fapidxx/wid2container.h>
......
...@@ -28,13 +28,12 @@ ...@@ -28,13 +28,12 @@
* *
* REVISIONS and CHANGES * REVISIONS and CHANGES
* - 04/07/2014 V1.0 Thomas Forbriger * - 04/07/2014 V1.0 Thomas Forbriger
* - 11/12/2019 V1.1 support negative time shift
* *
* ============================================================================ * ============================================================================
*/ */
#define TF_FAPID_SFF_MODWID2SHIFT_CC_VERSION \ #define TF_FAPID_SFF_MODWID2SHIFT_CC_VERSION \
"TF_FAPID_SFF_MODWID2SHIFT_CC V1.0 " "TF_FAPID_SFF_MODWID2SHIFT_CC V1.1"
#define TF_FAPID_SFF_MODWID2SHIFT_CC_CVSID \
"$Id: $"
#include <fapidxx/fapidsff.h> #include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h> #include <fapidxx/wid2container.h>
...@@ -65,7 +64,11 @@ int sff_modwid2shift__(char *wid2line, real *tmin, real * ...@@ -65,7 +64,11 @@ int sff_modwid2shift__(char *wid2line, real *tmin, real *
{ {
fapidxx::WID2container wid2c(wid2line, wid2line_len); fapidxx::WID2container wid2c(wid2line, wid2line_len);
libtime::TAbsoluteTime date(wid2c.wid2.date); libtime::TAbsoluteTime date(wid2c.wid2.date);
wid2c.wid2.date+=libtime::double2time(double(*tsec)+60.*double(*tmin)); double tshift=double(*tsec)+60.*double(*tmin);
if (tshift < 0.)
{ wid2c.wid2.date-=libtime::double2time(-tshift); }
else
{ wid2c.wid2.date+=libtime::double2time(tshift); }
wid2c.encode(wid2line, wid2line_len); wid2c.encode(wid2line, wid2line_len);
return 0; return 0;
} // int sff_modwid2shift__ } // int sff_modwid2shift__
......
...@@ -32,9 +32,7 @@ ...@@ -32,9 +32,7 @@
* ============================================================================ * ============================================================================
*/ */
#define TF_FAPID_SFF_MODWID2TIME_CC_VERSION \ #define TF_FAPID_SFF_MODWID2TIME_CC_VERSION \
"TF_FAPID_SFF_MODWID2TIME_CC V1.0 " "TF_FAPID_SFF_MODWID2TIME_CC V1.0"
#define TF_FAPID_SFF_MODWID2TIME_CC_CVSID \
"$Id: $"
#include <fapidxx/fapidsff.h> #include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h> #include <fapidxx/wid2container.h>
......