Commit 86b4b6de authored by thomas.forbriger's avatar thomas.forbriger

gresy [WP]: implement Fourier domain filter

parent 1e7a37d1
......@@ -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.
......
......@@ -129,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
......@@ -301,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'
......@@ -424,6 +436,29 @@ 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))))
print *,om(io),response(io)
enddo
c
c----------------------------------------------------------------------
c
c open seismogram file
c
......@@ -672,13 +707,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
do io=1,nom
sdata(io)=sdata(io)*response(io)
enddo
endif
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
c
c apply time shift
do io=1,nom
......@@ -719,6 +753,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)
......
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