...
 
Commits (7)
......@@ -3,8 +3,10 @@ this is <CHANGELOG>
Recent development in Seitosh (bug fixes, new features, etc)
------------------------------------------------------------
11.12.2019: grepg
- commit 68295a7: provide option to select plot axis range
11.12.2019: grepg, libfapidxx, gresy
- 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
- commit 27491f4: let RNGgaussian support selection of rng type
......
......@@ -105,7 +105,7 @@ gresynoise: gresynoise.o nyquist_check.o
gresy: gresy.o nyquist_check.o
$(FC) -o $@ $^ -lgrrefsub -lrefpar -lsffu \
-ltf -lsff -ltime -L$(LOCLIBDIR)
-ltf -lsff -ltime -lfourier -L$(LOCLIBDIR)
rhesyg: rhesyg.o
$(FC) -o rhesyg rhesyg.o -lgrrefsub -lrheology -lrefpar \
......@@ -119,7 +119,8 @@ afehsyg: syg.o
gresyx: gresy.o nyquist_check.o
$(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
......
......@@ -76,6 +76,7 @@ Dependencies:
Seitosh libraries required to compile the code:
libaff libdatrwxx libfapidxx libgrrefsub libgsexx librefpar librheology
libsff libsffu libsffxx libtf libtime libtime++ libwrefsub libtfxx
libfourier
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
packages.
......
......@@ -2,6 +2,7 @@ c this is <gresy.f>
c------------------------------------------------------------------------------
c
c Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart)
c Copyright 2019 by Thomas Forbriger (BFO, KIT)
c
c calculate seismograms from a greens function matrix file
c
......@@ -57,12 +58,13 @@ c 29/11/2011 V1.18 implemented radial component line source
c both new cases successfully tested against refmet
c notice: refmet is not yet tested against a
c reference
c 11/12/2019 V1.19 support Fourier domain filters
c
c==============================================================================
c
program gresy
character*79 version
parameter(version='GRESY V1.18 GREens function SYnthetics')
parameter(version='GRESY V1.19 GREens function SYnthetics')
c dimensions
integer maxtr, maxsamp, maxom, maxu
parameter(maxu=10000, maxtr=maxu, maxom=4100, maxsamp=maxom*2)
......@@ -87,7 +89,7 @@ c functions
real sffu_offset, sffu_tfirst
c free block
integer maxfree, nfree
parameter(maxfree=5)
parameter(maxfree=15)
character*80 free(maxfree)
c SFF parameters
real sffversion, tanf
......@@ -105,14 +107,18 @@ c seismogram trace to sff file
real dt
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=15)
character*3 optid(maxopt)
parameter(maxopt=19)
character*4 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, usemaster
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
c taper
real tf_costap
......@@ -123,6 +129,8 @@ c response file
integer nresp
real omega(maxom)
real domega
c filter function
double complex fou_eval
c any
integer mlu,lu, i, ierr, io, iu
double precision pi2, arg, du, dom, scal, pi
......@@ -135,11 +143,12 @@ c functions
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','-dt', '-N', '-m'/
& '-L','-ty','-dt', '-N', '-m', '-lpb', '-hpb',
& '-tss', '-tso'/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,2*.TRUE.,6*.FALSE.,
& 3*.true.,.false./
& 3*.true.,.false.,4*.true./
data optarg/'-','1.','-','10.','junk',6*'-','sff',
& '1.e-18','1000000','-'/
& '1.e-18','1000000','-',2*'1.,0',2*'0.'/
c----------------------------------------------------------------------
c
c read commandline
......@@ -149,6 +158,7 @@ c
print *,' [-m] [-l lambda] [-o]'
print *,' [-t frac] [-r respfile] [-S] [-1] [-2] [-R]'
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 -xhelp'
call tf_seitosh_reference
......@@ -207,6 +217,20 @@ c
print *,'-L simulate seismograms from line source (2D)'
print *,'-R calculate radial component'
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 *,'The Fourier transformation used in this program and in'
print *,'related programs (like gremlin, syg, and greda) is'
......@@ -279,7 +303,17 @@ c
if (dt.lt.1.e-20) stop 'ERROR: sampling interval too small'
read(optarg(14), *, err=99) nsampwrite
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
call sff_select_output_format(fileformat, ierr)
if (ierr.ne.0) stop 'ERROR: selecting output file format'
......@@ -402,6 +436,28 @@ c
endif
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 open seismogram file
c
......@@ -416,6 +472,26 @@ c
& //rcvname(1:index(rcvname, ' '))
free(4)=rcvtext
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'
scs='C'
sc1=0.
......@@ -424,12 +500,24 @@ c
date='990418'
time='000000.000'
endif
if (verbose) then
print *,'open output file ',seisname(1:index(seisname, ' '))
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,
& scs, sc1 ,sc2 ,sc3,
& date, time, ierr)
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 preparatory calculations
......@@ -650,13 +738,12 @@ c----------------------------------------------------------------------
do io=nom+1,nsamp-nom+1
sdata(io)=(0.d0,0.d0)
enddo
c apply response (Fourier coefficients of source wavelet)
c -------------------------------------------------------
if (optresponse) then
c apply response
c (Fourier coefficients of source wavelet and Fourier domain filter)
c ------------------------------------------------------------------
do io=1,nom
sdata(io)=sdata(io)*response(io)
enddo
endif
c
c apply time shift
do io=1,nom
......@@ -697,6 +784,11 @@ c prepare wid2line
cs='C'
nstack=0
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'
call sff_WTraceI(lu, wid2line, nsampwrite, fdata, idata, last,
& cs, c1, c2, c3, nstack, ierr)
......
......@@ -32,7 +32,7 @@
* ============================================================================
*/
#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/wid2container.h>
......
......@@ -32,9 +32,7 @@
* ============================================================================
*/
#define TF_FAPID_SFF_MODWID2SAMPRAT_CC_VERSION \
"TF_FAPID_SFF_MODWID2SAMPRAT_CC V1.0 "
#define TF_FAPID_SFF_MODWID2SAMPRAT_CC_CVSID \
"$Id: $"
"TF_FAPID_SFF_MODWID2SAMPRAT_CC V1.0"
#include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h>
......
......@@ -32,9 +32,7 @@
* ============================================================================
*/
#define TF_FAPID_SFF_MODWID2SAMPS_CC_VERSION \
"TF_FAPID_SFF_MODWID2SAMPS_CC V1.0 "
#define TF_FAPID_SFF_MODWID2SAMPS_CC_CVSID \
"$Id: $"
"TF_FAPID_SFF_MODWID2SAMPS_CC V1.0"
#include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h>
......
......@@ -28,13 +28,12 @@
*
* REVISIONS and CHANGES
* - 04/07/2014 V1.0 Thomas Forbriger
* - 11/12/2019 V1.1 support negative time shift
*
* ============================================================================
*/
#define TF_FAPID_SFF_MODWID2SHIFT_CC_VERSION \
"TF_FAPID_SFF_MODWID2SHIFT_CC V1.0 "
#define TF_FAPID_SFF_MODWID2SHIFT_CC_CVSID \
"$Id: $"
"TF_FAPID_SFF_MODWID2SHIFT_CC V1.1"
#include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h>
......@@ -65,7 +64,11 @@ int sff_modwid2shift__(char *wid2line, real *tmin, real *
{
fapidxx::WID2container wid2c(wid2line, wid2line_len);
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);
return 0;
} // int sff_modwid2shift__
......
......@@ -32,9 +32,7 @@
* ============================================================================
*/
#define TF_FAPID_SFF_MODWID2TIME_CC_VERSION \
"TF_FAPID_SFF_MODWID2TIME_CC V1.0 "
#define TF_FAPID_SFF_MODWID2TIME_CC_CVSID \
"$Id: $"
"TF_FAPID_SFF_MODWID2TIME_CC V1.0"
#include <fapidxx/fapidsff.h>
#include <fapidxx/wid2container.h>
......