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

proceeding

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: 2287
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 050077d2
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.8 2007-05-21 10:08:44 tforb Exp $)
c ($Id: gresynoise.f,v 1.9 2007-05-22 11:51:59 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -37,10 +37,10 @@ c nset selected through command line
c must be: nset <= maxset
c nom number of Fourier coefficients to store in each set
c selected through Green's coefficient file
c nsamp number of samples to obtain from each set
c largest values with nsamp <= 2*maxom
c noutsamp number of output samples
c (nset+1)*nsamp/2
c nsetsamp number of samples to obtain from each set
c largest values with nsetsamp <= 2*maxom
c nsamp number of output samples
c (nset+1)*nsetsamp/2
c
program gresynoise
c
......@@ -49,13 +49,14 @@ c
&'GRESYNOISE V1.4i calculate noise seismograms')
character*(*) GRESYNOISE_CVS_ID
parameter(GRESYNOISE_CVS_ID=
&'$Id: gresynoise.f,v 1.8 2007-05-21 10:08:44 tforb Exp $')
&'$Id: gresynoise.f,v 1.9 2007-05-22 11:51:59 tforb Exp $')
c
c dimensions
integer maxtr, maxsamp, maxom, maxu
integer maxtr, maxsetsamp, maxom, maxu, maxset, maxsamp
c parameter(maxu=1000, maxtr=10000, maxom=4100, maxsamp=maxom*2)
parameter(maxu=5000, maxtr=10000, maxom=4100, maxsamp=100000)
integer ntr, nsamp, Znom, Znu, Rnom, Rnu
parameter(maxu=5000, maxtr=10000, maxom=4100)
parameter(maxset=20, maxsamp=(maxset+1)*maxom, maxsetsamp=2*maxom)
integer ntr, nsamp, Znom, Znu, Rnom, Rnu, nset, nsetsamp
c greens function
character*78 greenname, Rgreenname
complex Zgreen(maxom, maxu)
......@@ -66,17 +67,19 @@ c receiver specifications
real*8 rcvvred, rcvtli, rcvtre
character*80 rcvtext
real*8 phi(maxtr)
real scalefac(maxtr)
real scalefac(maxtr,maxset)
c seismogram parameters
integer pown
character*80 Zseisname,rcvname,Rseisname
real*8 r(maxtr)
complex*16 Zsdata(maxsamp)
complex*16 Ztrans(maxsetsamp)
complex*16 Zsets(maxom,maxset)
complex*16 Zfourier(maxom)
real Zfdata(maxsamp)
integer Zidata(maxsamp)
equivalence(Zfdata, Zidata)
complex*16 Rsdata(maxsamp)
complex*16 Rtrans(maxsetsamp)
complex*16 Rsets(maxom,maxset)
complex*16 Rfourier(maxom)
real Rfdata(maxsamp)
integer Ridata(maxsamp)
......@@ -102,7 +105,7 @@ c source response
complex ime
parameter(ime=(0.,1.))
c any
integer lu, i, iargc, ierr, io, iu
integer lu, i, iargc, ierr, io, iu, is
real*8 pi2, arg, du, dom, scal, pi
parameter(lu=20,pi=3.141592653589793115997)
c functions
......@@ -110,9 +113,9 @@ c functions
real*8 tf_dj1, tf_dy1
c options
logical debug, optlambda, optnew, optresponse, hankel1, hankel2
logical suppress, optrndscale, optnsamp
logical suppress, optrndscale
real lambdalim, tapfrac, scalexp
integer argnsamp
integer argnsets
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=11)
......@@ -124,7 +127,7 @@ c here are the keys to our commandline options
& 2h-r,2h-n/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,.TRUE.,4*.FALSE.,
& .TRUE.,.FALSE.,.TRUE./
data optarg/1h-,2h1.,1h-,3h10.,4*1h-,2h0.,2*1h-/
data optarg/1h-,2h1.,1h-,3h10.,4*1h-,2h0.,2*1h1/
c------------------------------------------------------------------------------
c basic information
c
......@@ -168,8 +171,9 @@ c
print *,'-S suppress zero frequency and zero slowness'
print *,'-1 use Hankel 1 instead of Bessel'
print *,'-2 use Hankel 2 instead of Bessel'
print *,'-n n set number of output samples to largest'
print *,' possible value, if possible larger than n'
print *,'-n n number of sets to use'
print *,' time series will become larger with each'
print *,' additional set of Fourier coefficients'
print *,' '
print *,'Input file units:'
print *,' '
......@@ -205,8 +209,11 @@ c
hankel2=optset(8)
read (optarg(9), *, err=99) scalexp
optrndscale=optset(10)
optnsamp=optset(11)
if (optnsamp) read (optarg(11), *, err=99) argnsamp
read (optarg(11), *, err=99) argnsets
argnsets=min(maxset, argnsets)
argnsets=max(1,argnsets)
print *,'number of sets to be used: ',argnsets
nset=argnsets
c
pi2=2.d0*pi
c
......@@ -276,14 +283,14 @@ c noise signal
enddo
c
c initialize arrays
do io=1,maxom
Zfourier(io)=(0.,0.)
Rfourier(io)=(0.,0.)
do is=1,maxset
do io=1,maxom
Zsets(io,is)=(0.,0.)
Rsets(io,is)=(0.,0.)
enddo
enddo
do io=1,maxsamp
Zsdata(io)=(0.,0.)
Zfdata(io)=0.
Rsdata(io)=(0.,0.)
Rfdata(io)=0.
enddo
c
......@@ -308,9 +315,11 @@ c
c scale down epicentral distances to meters
do i=1,ntr
r(i)=r(i)*1000.
scalefac(i)=1.
do is=1,nset
scalefac(i,is)=1.
enddo
enddo
if (optrndscale) call tf_gsl_rng_ugaussian(scalefac, ntr)
if (optrndscale) call tf_gsl_rng_ugaussian(scalefac, ntr*nset)
c check stepwidth
dom=Zom(2)-Zom(1)
do io=2,Znom
......@@ -328,17 +337,19 @@ c
& 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
if (.not.optnsamp) argnsamp=2*maxom
pown=0
1 continue
pown=pown+1
nsamp=2**pown
if (nsamp.lt.argnsamp) goto 1
if (nsamp.lt.(2*Znom)) stop 'ERROR: cannot use enough samples'
nsetsamp=2**pown
if (nsetsamp.lt.maxsetsamp) goto 1
if (nsetsamp.lt.(2*Znom)) stop 'ERROR: cannot use enough samples'
c calculate sampling interval
dt=pi2/(nsamp*dom)
print *,'traces will have ',nsamp,' samples at ',dt,
dt=pi2/(nsetsamp*dom)
nsamp=(nset+1)*nsetsamp/2
print *,'sets will have ',nsetsamp,' samples at ',dt,
& 's sampling interval'
print *,'final output series will have ',nsamp,
& ' samples'
c
if (hankel1) then
print *,'I will use the first Hankel function'
......@@ -427,21 +438,24 @@ c end of frequency loop
c----------------------------------------------------------------------
if (debug) print *,'left loops'
c apply response
if (optresponse) then
call tf_gsl_rng_uniform(phase, Znom)
do is=1,nset
if (optresponse) then
call tf_gsl_rng_uniform(phase, Znom)
do io=1,Znom
Zfourier(io)=Zfourier(io)*cdexp(ime*pi2*phase(io))
Rfourier(io)=Rfourier(io)*cdexp(ime*pi2*phase(io))
cc print *,i,io,Zfourier(io),Rfourier(io),phase(io)
enddo
endif
c scale trace and stack
scal=1.d0/(sqrt(float(nsamp))*dt)*(r(i)**scalexp)
scal=scal*scalefac(i,is)
cc print *,'rec # ', i, ' scal: ',scal
do io=1,Znom
Zfourier(io)=Zfourier(io)*cdexp(ime*pi2*phase(io))
Rfourier(io)=Rfourier(io)*cdexp(ime*pi2*phase(io))
cc print *,i,io,Zfourier(io),Rfourier(io),phase(io)
print *,i,io,is,Zfourier(io),Rfourier(io),scal
Zsets(io,is)=Zsets(io,is)+Zfourier(io)*scal
Rsets(io,is)=Zsets(io,is)+Rfourier(io)*scal
enddo
endif
c scale trace and stack
scal=1.d0/(sqrt(float(nsamp))*dt)*(r(i)**scalexp)
scal=scal*scalefac(i)
cc print *,'rec # ', i, ' scal: ',scal
do io=1,Znom
Zsdata(io)=Zsdata(io)+Zfourier(io)*scal
Rsdata(io)=Zsdata(io)+Rfourier(io)*scal
enddo
enddo
c receiver loop ends here
......@@ -451,34 +465,51 @@ c print *,Zsdata(io),Rsdata(io)
c enddo
c
print *,'make Fourier coefficients complete...'
print *,'transform sets and stack...'
do io=1,nsamp
Zfdata(io)=0.
Rfdata(io)=0.
enddo
c
do is=1,nset
c make Fourier coefficients complete
do io=Znom+1,nsamp-Znom+1
Zsdata(io)=(0.d0,0.d0)
Rsdata(io)=(0.d0,0.d0)
enddo
do io=1,Znom-1
Zsdata(nsamp-io+1)=conjg(Zsdata(io+1))
Rsdata(nsamp-io+1)=conjg(Rsdata(io+1))
enddo
do io=1,maxom
Ztrans(io)=(0.d0,0.d0)
Rtrans(io)=(0.d0,0.d0)
enddo
do io=1,Znom
Ztrans(io)=Zsets(io,is)
Rtrans(io)=Rsets(io,is)
enddo
do io=1,Znom-1
Ztrans(nsamp-io+1)=conjg(Ztrans(io+1))
Rtrans(nsamp-io+1)=conjg(Rtrans(io+1))
enddo
c
print *,'transform to time domain...'
print *,'transform set ',is,' to time domain...'
c transform to time domain
if (debug) print *,'go fork'
call tf_dfork(nsamp, Zsdata, 1.d0)
call tf_dfork(nsamp, Rsdata, 1.d0)
if (debug) print *,'done fork'
if (debug) print *,'go fork'
call tf_dfork(nsetsamp, Ztrans, 1.d0)
call tf_dfork(nsetsamp, Rtrans, 1.d0)
if (debug) print *,'done fork'
c extract time series
do io=1,nsamp
Zfdata(io)=Zfdata(io)+real(Zsdata(io))
if ((debug).and.
& (abs(imag(Zsdata(io))).gt.(0.01*abs(real(Zsdata(io))))))
& print *,'WHERE does this imaginary part come from?'
Rfdata(io)=Rfdata(io)+real(Rsdata(io))
if ((debug).and.
& (abs(imag(Rsdata(io))).gt.(0.01*abs(real(Rsdata(io))))))
& print *,'WHERE does this imaginary part come from?'
c print *,Zfdata(io),Rfdata(io)
ltap=nsetsamp/2
rtap=ltap+1
if (is.eq.1) ltap=1
if (is.eq.nset) rtap=nsetsamp
do io=1,nsetsamp
scal=tf_costap(io,1,ltap,rtap,nsetsamp)
i=io+is*(nsetsamp/2)
Zfdata(i)=Zfdata(i)+real(Ztrans(io))*scal
if ((debug).and.
& (abs(imag(Ztrans(io))).gt.(0.01*abs(real(Ztrans(io))))))
& print *,'WHERE does this imaginary part come from?'
Rfdata(i)=Rfdata(i)+real(Rtrans(io))*scal
if ((debug).and.
& (abs(imag(Rtrans(io))).gt.(0.01*abs(real(Rtrans(io))))))
& print *,'WHERE does this imaginary part come from?'
c print *,Zfdata(io),Rfdata(io)
enddo
enddo
c----------------------------------------------------------------------
c
......
Supports Markdown
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