Commit 49b14a9b authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

correction

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: 2283
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 85850918
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.5 2007-05-15 14:59:32 tforb Exp $)
c ($Id: gresynoise.f,v 1.6 2007-05-15 16:05:16 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -12,6 +12,7 @@ c REVISIONS and CHANGES
c 09/05/2007 V1.0 Thomas Forbriger
c 10/05/2007 V1.1 introduced randomization
c 11/05/2007 V1.2 allow selection of number of output samples
c 15/05/2007 V1.3 reordered loops; phase is of correct type now
c
c ============================================================================
c
......@@ -19,10 +20,10 @@ c
c
character*(*) version
parameter(version=
&'GRESYNOISE V1.2 calculate noise seismograms')
&'GRESYNOISE V1.3 calculate noise seismograms')
character*(*) GRESYNOISE_CVS_ID
parameter(GRESYNOISE_CVS_ID=
&'$Id: gresynoise.f,v 1.5 2007-05-15 14:59:32 tforb Exp $')
&'$Id: gresynoise.f,v 1.6 2007-05-15 16:05:16 tforb Exp $')
c
c dimensions
integer maxtr, maxsamp, maxom, maxu
......@@ -45,10 +46,12 @@ c seismogram parameters
character*80 Zseisname,rcvname,Rseisname
real*8 r(maxtr)
complex*16 Zsdata(maxsamp)
complex*16 Zfourier(maxom)
real Zfdata(maxsamp)
integer Zidata(maxsamp)
equivalence(Zfdata, Zidata)
complex*16 Rsdata(maxsamp)
complex*16 Rfourier(maxom)
real Rfdata(maxsamp)
integer Ridata(maxsamp)
equivalence(Rfdata, Ridata)
......@@ -69,7 +72,7 @@ c taper
real tf_costap
integer ltap, rtap
c source response
real phase(maxom)
double precision phase(maxom)
complex ime
parameter(ime=(0.,1.))
c any
......@@ -246,15 +249,24 @@ c noise signal
& stop 'ERROR: inconsistent slowness'
enddo
c
c initialize arrays
do io=1,maxom
Zfourier(io)=(0.,0.)
Rfourier(io)=(0.,0.)
enddo
do io=1,maxsamp
Zsdata(io)=(0.,0.)
Zfdata(io)=0.
Rsdata(io)=(0.,0.)
Rfdata(io)=0.
enddo
c
c suppress zero frequency and slowness if requested
if (suppress) then
print *,'suppress zero frequency and zero slowness coefficients'
do io=1,maxom
Zgreen(io,1)=(0.,0.)
Rgreen(io,1)=(0.,0.)
Zsdata(io)=(0.,0.)
Zfdata(io)=0.
Rsdata(io)=(0.,0.)
Rfdata(io)=0.
enddo
do iu=1,maxu
Zgreen(1,iu)=(0.,0.)
......@@ -310,6 +322,9 @@ c
print *,'I will use the Bessel function'
endif
last=.false.
c
c start receiver loop
c -------------------
do i=1,ntr
print *,'working on trace ',i,' at ',r(i),'m'
if (i.eq.ntr) last=.true.
......@@ -338,87 +353,111 @@ c set slowness taper
print *,'u-taper for ',Zom(io)/pi2,'Hz: ',ltap,'-',rtap,
& ' ',Zslo(ltap),'s/m-',Zslo(rtap),'s/m'
endif
Zsdata(io)=(0.d0,0.d0)
Rsdata(io)=(0.d0,0.d0)
Zfourier(io)=(0.d0,0.d0)
Rfourier(io)=(0.d0,0.d0)
c
c start slowness loop (conditional: select base function)
c -------------------------------------------------------
if (debug) print *,'enter slowness loop'
if (hankel1) then
do iu=1,rtap
arg=max(Zslo(iu)*Zom(io)*r(i),1.d-100)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
Rfourier(io)=Rfourier(io)+Rgreen(io,iu)*
& (tf_dj1(arg)+(0.d0,1.d0)*tf_dy1(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
Zfourier(io)=Zfourier(io)+Zgreen(io,iu)*
& (tf_dj0(arg)+(0.d0,1.d0)*tf_dy0(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
Zsdata(io)=0.5d0*Zsdata(io)
Rsdata(io)=0.5d0*Rsdata(io)
Zfourier(io)=0.5d0*Zfourier(io)
Rfourier(io)=0.5d0*Rfourier(io)
elseif (hankel2) then
do iu=1,rtap
arg=max(Zslo(iu)*Zom(io)*r(i),1.d-100)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
Rfourier(io)=Rfourier(io)+Rgreen(io,iu)*
& (tf_dj1(arg)-(0.d0,1.d0)*tf_dy1(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
Zfourier(io)=Zfourier(io)+Zgreen(io,iu)*
& (tf_dj0(arg)-(0.d0,1.d0)*tf_dy0(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
Zsdata(io)=0.5d0*Zsdata(io)
Rsdata(io)=0.5d0*Rsdata(io)
Zfourier(io)=0.5d0*Zfourier(io)
Rfourier(io)=0.5d0*Rfourier(io)
else
do iu=1,rtap
arg=Zslo(iu)*Zom(io)*r(i)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
Rfourier(io)=Rfourier(io)+Rgreen(io,iu)*
& tf_dj1(arg)*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
Zfourier(io)=Zfourier(io)+Zgreen(io,iu)*
& tf_dj0(arg)*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
endif
c end slowness loop
c -----------------
enddo
c end of frequency loop for Z-component
c end of frequency loop
c----------------------------------------------------------------------
if (debug) print *,'left loops'
do io=Znom+1,nsamp-Znom+1
Zsdata(io)=(0.d0,0.d0)
Rsdata(io)=(0.d0,0.d0)
enddo
c apply response
if (optresponse) then
call tf_gsl_rng_uniform(phase, Znom)
do io=1,Znom
Zsdata(io)=Zsdata(io)*exp(ime*pi2*phase(io))
Rsdata(io)=Rsdata(io)*exp(ime*pi2*phase(io))
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
do io=1,Znom-1
Zsdata(nsamp-io+1)=conjg(Zsdata(io+1))
Rsdata(nsamp-io+1)=conjg(Rsdata(io+1))
enddo
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'
c move trace to writing array and scale to impulse response
c scale trace and stack
scal=1.d0/(sqrt(float(nsamp))*dt)*(r(i)**scalexp)
scal=scal*scalefac(i)
if (debug) print *,'stack traces'
do io=1,nsamp
Zfdata(io)=Zfdata(io)+real(Zsdata(io))*scal
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))*scal
if ((debug).and.
& (abs(imag(Rsdata(io))).gt.(0.01*abs(real(Rsdata(io))))))
& print *,'WHERE does this imaginary part come from?'
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
c -----------------------
c do io=1,Znom
c print *,Zsdata(io),Rsdata(io)
c enddo
c
print *,'make Fourier coefficients complete...'
c
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
c
print *,'transform 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'
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)
enddo
c----------------------------------------------------------------------
c
c open seismogram file
print *,'write seismograms...'
c
free(1)=version
free(2)='green file '//greenname(1:index(greenname, ' '))
......@@ -553,6 +592,5 @@ c phase slowness is given in s/m
96 stop 'ERROR: closing green file'
end
c ----- END OF gresy.f -----
c
c ----- END OF gresynoise.f -----
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