Commit 3787e51f authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

reintegrate merge (from branches/greda): greda now supports

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.
transformation of radial component seismograms


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 4370
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent b3251288
......@@ -106,6 +106,17 @@ c - use GSL instead of numerical recipes
c 30.12.2010 V3.21 - implemented libfapidxx
c 12/01/2011 V3.21b - added definition of Fourier transformation
c online help text
c 10/01/2012 V3.22 - implemented radial component
c tested with Fourier-Bessel transformation and
c HOP inversion both with Bessel kernel and
c Hankel(1) kernel; both work well or at least as
c well as reconstruction tests for vertical
c components (HOP apparently looses some signal
c energy due to damping); K-damping, exp-damping, and
c boxcar damping all producde nuemrical errors in
c the test; slant stack analysis does not
c distinguish between vertical and radial
c component
c
c==============================================================================
c
......@@ -115,7 +126,7 @@ c first we declare some general variables
c
character*79 version,CVSID
parameter(version=
&'GREDA V3.21b Calculate Fourier-Bessel expansion coefficients')
&'GREDA V3.22 Calculate Fourier-Bessel expansion coefficients')
parameter(CVSID=
& '$Id$')
c
......@@ -167,7 +178,7 @@ c file
c
c selected method
character*70 method
character*2 kernel
character*4 kernel
c
c----------------------------------------------------------------------
c here we go with command line parameters
......@@ -187,7 +198,7 @@ c parameters
logical pwofile, pwoautofile, pwosel
c command line
integer maxopt, lastarg, iargc
parameter(maxopt=36)
parameter(maxopt=37)
character*4 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
......@@ -195,20 +206,20 @@ c here are the keys to our commandline options
data optid/2h-d,2h-t,2h-s,2h-n,2h-f,2h-o,2h-1,2h-2,2h-v,2h-T,2h-S,
& 2h-M,2h-L,2h-q,2h-E,2h-K,2h-D,2h-a,2h-b,2h-g,2h-H,2h-O,2h-Q,2h-B,
& 2h-P,2h-r,2h-F,2h-N,2h-X,2h-W,4h-tap,4h-spo,3h-pw,4h-pwf,4h-pwa,
& '-ty'/
& '-ty','-R'/
data opthasarg/.FALSE.,4*.TRUE.,4*.FALSE.,.TRUE.,3*.FALSE.,2*.TRUE.,
& 6*.FALSE.,2*.TRUE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,.TRUE.,.FALSE.,
& 7*.TRUE./
& 7*.TRUE.,.false./
data optarg/1h-,3h10.,2h8.,1h5,4h100.,4*1h-,2h0.,3*1h-,5h1.,1.,
& 2h1.,6*1h-,4h0.01,5h1.,1.,4*1h-,3h10.,1h-,6h1.,10.,
& 11h0.,0.,0.,0.,4*3hnil,'sff'/
& 11h0.,0.,0.,0.,4*3hnil,'sff','-'/
c
c======================================================================
c
c go on with executable code
c
print *,version
print *,'Usage: greda datafile coeffile [-ty format]'
print *,'Usage: greda datafile coeffile [-ty format] [-R]'
print *,' [-L] [-K] [-P] [-D] [-M] [-H]'
print *,' [-n nslo] [-s smax] [-f fmax]'
print *,' [-t frac] [-T frac] [-E edge]'
......@@ -247,6 +258,12 @@ c
print *,'coeffile Name of file to contain results.'
print *,'-ty format input file format (see list below)'
print *,' '
print *,'-R radial component seismograms of a vertical'
print *,' single force or an explosion are to be'
print *,' transformed.'
print *,' (tested for HOP-inversion and'
print *,' Fourier-Bessel transformation)'
print *,' '
print *,'Available alternatives to the Fourier-Bessel transform:'
print *,'-L Calculate by linear inversion (exp-damping).'
print *,'-K Calculate by linear inversion (K_0-damping).'
......@@ -524,6 +541,7 @@ c print *,'DEBUG: returned from tf_cmdline'
pwoautofile=optset(35)
pwafilename=optarg(35)
informat=optarg(36)
radial=optset(37)
if (debug) print *,'DEBUG: read options'
c
if ((.not.(planewave)).and.(smax.lt.0.))
......@@ -537,16 +555,30 @@ c print *,'DEBUG: returned from pwo_init'
c----------------------------------------------------------------------
c report
print *,' '
if (hankel1) then
print *,'I will use the hankel function H^(1)_0 of order zero'
kernel='H1'
elseif (hankel2) then
print *,'I will use the hankel function H^(2)_0 of order zero'
kernel='H2'
if (radial) then
if (hankel1) then
print *,'I will use the Hankel function H^(1)_1 of order one'
kernel='H1_1'
elseif (hankel2) then
print *,'I will use the Hankel function H^(2)_1 of order one'
kernel='H2_1'
else
print *,'I will use the Bessel function of the first kind',
& ' J_1 of order one'
kernel='J_1'
endif
else
print *,'I will use the Bessel function of the first kind J_0',
& ' of order zero'
kernel='J'
if (hankel1) then
print *,'I will use the Hankel function H^(1)_0 of order zero'
kernel='H1_0'
elseif (hankel2) then
print *,'I will use the Hankel function H^(2)_0 of order zero'
kernel='H2_0'
else
print *,'I will use the Bessel function of the first kind J_0',
& ' of order zero'
kernel='J_0'
endif
endif
print *,' '
c
......@@ -1624,6 +1656,11 @@ c
c
c expand the model from coefficients alpha using
c J_0 or H^1_0 or H^2_0 as representers
c or
c J_1 or H^1_1 or H^2_1 as representers, respectively
c
c This subroutine is used by almost all methods, including the
c Fourier-Bessel transformation
c
c get common block
include 'greda_dim.inc'
......@@ -1636,6 +1673,7 @@ c
c
integer it, iu
double precision tf_dj0, arg, tf_dy0
double precision tf_dj1, tf_dy1
real thisom, slowness
double complex sum, add
c
......@@ -1646,28 +1684,53 @@ c evaluate expansion for every slowness
if (applywltaper) then
call setwltaper(omega,slowness,r)
endif
do it=1,ntr
arg=slowness*omega*r(it)
add=(0.d0,0.d0)
if ((uzerospecial).and.(iu.eq.1)) then
if (radial) then
do it=1,ntr
arg=slowness*omega*r(it)
add=(0.d0,0.d0)
if ((uzerospecial).and.(iu.eq.1)) then
c use J_1
arg=alpha(it)*tf_dj1(arg)*wltaperfac(it)
elseif (hankel1) then
c use H^(1)_1
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj1(arg)+(0.d0,1.d0)*tf_dy1(arg))*
& wltaperfac(it)*0.5
elseif (hankel2) then
c use H^(2)_1
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj1(arg)-(0.d0,1.d0)*tf_dy1(arg))*
& wltaperfac(it)*0.5
else
c use J_1
add=alpha(it)*tf_dj1(arg)*wltaperfac(it)
endif
sum=sum+add
enddo
else
do it=1,ntr
arg=slowness*omega*r(it)
add=(0.d0,0.d0)
if ((uzerospecial).and.(iu.eq.1)) then
c use J_0
arg=alpha(it)*tf_dj0(arg)*wltaperfac(it)
elseif (hankel1) then
arg=alpha(it)*tf_dj0(arg)*wltaperfac(it)
elseif (hankel1) then
c use H^(1)_0
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj0(arg)+(0.d0,1.d0)*tf_dy0(arg))*
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj0(arg)+(0.d0,1.d0)*tf_dy0(arg))*
& wltaperfac(it)*0.5
elseif (hankel2) then
elseif (hankel2) then
c use H^(2)_0
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj0(arg)-(0.d0,1.d0)*tf_dy0(arg))*
arg=max(arg,1.d-100)
add=alpha(it)*(tf_dj0(arg)-(0.d0,1.d0)*tf_dy0(arg))*
& wltaperfac(it)*0.5
else
else
c use J_0
add=alpha(it)*tf_dj0(arg)*wltaperfac(it)
endif
sum=sum+add
enddo
add=alpha(it)*tf_dj0(arg)*wltaperfac(it)
endif
sum=sum+add
enddo
endif
ogreen(iu)=sum
enddo
c
......@@ -1712,6 +1775,11 @@ c
c----------------------------------------------------------------------
c
c calculate model expansion coefficients alpha from gram matrix
c
c This subroutine solves the system of linear equations set up with the
c Gram matrix. All methods of linear inversion using the Gram matrix
c make use of this subroutine, except the HOP-method, which uses an
c analytical expression for the Gram matrix.
c
subroutine modexp(iom, s)
c
......@@ -1768,9 +1836,13 @@ c here we go with the various strategies
c
c======================================================================
c
c first we build code for a simple Fourier Bessel transform
c Method: straight foward Fourier-Bessel transformation
c -----------------------------------------------------
c
c This subsroutine prepares values to be used in subroutine expmodel to
c calculate Fourier-Bessel expansion coefficients by a straight forward
c Fourier-Bessel transformation.
c
c----------------------------------------------------------------------
c
subroutine backcoeff(verbose, r, omega, spect, chain, firstrec,
& numbers)
......@@ -1852,8 +1924,9 @@ c
end
c
c======================================================================
c
c do a linear inversion with exp(-rohq*uq) damping factor
c
c Method: Linear inversion using a Gram matrix formulation and
c exponential damping with exp(-rohq*uq) damping factor
c
c----------------------------------------------------------------------
c
......@@ -1870,6 +1943,7 @@ c
c
integer i,j
double precision tf_gsl_sf_bessel_I0
double precision tf_gsl_sf_bessel_I1
double precision mwr, wr2, er2, earg, ibarg
c
er2=0.5d0/rhoq
......@@ -1877,13 +1951,23 @@ c
wr2=omq*er2
c
c calculate upper half and diagonal
do i=1,ntr
do j=i,ntr
earg=mwr*(r(i)**2+r(j)**2)
ibarg=wr2*r(i)*r(j)
gram(i,j)=er2*exp(earg)*tf_gsl_sf_bessel_I0(ibarg)
if (radial) then
do i=1,ntr
do j=i,ntr
earg=mwr*(r(i)**2+r(j)**2)
ibarg=wr2*r(i)*r(j)
gram(i,j)=er2*exp(earg)*tf_gsl_sf_bessel_I1(ibarg)
enddo
enddo
enddo
else
do i=1,ntr
do j=i,ntr
earg=mwr*(r(i)**2+r(j)**2)
ibarg=wr2*r(i)*r(j)
gram(i,j)=er2*exp(earg)*tf_gsl_sf_bessel_I0(ibarg)
enddo
enddo
endif
c
c matrix is symmetric
do j=1,ntr-1
......@@ -1932,9 +2016,12 @@ c
c
c======================================================================
c
c Method: slant stack
c -------------------
c
c assume plane waves and do it by stacking seismograms
c
c----------------------------------------------------------------------
c This method cannot distinguish between vertical and radial component.
c
subroutine planestack(omega, slo, r, green, iomega, applywltaper)
c
......@@ -1986,10 +2073,11 @@ c
c
c======================================================================
c
c do a direct matrix inversion
c Method: direct discrete matrix inversion
c ----------------------------------------
c
c This is a DIRTY APPROACH and will most probably fail
c
c----------------------------------------------------------------------
c
c calculate forward matrix
c
subroutine forwardmat(omega, slo, r, hankel1, hankel2, uzerospecial)
......@@ -2041,8 +2129,9 @@ c
end
c
c======================================================================
c
c do a linear inversion with K_0 damping factor
c
c Method: Linear inversion using a Gram matrix formulation and
c damping with K_0 damping factor
c
c----------------------------------------------------------------------
c
......@@ -2090,39 +2179,62 @@ c
c
c----------------------------------------------------------------------
c
subroutine scalkmet(nslo, rho, ogreen, slo, ntr)
subroutine scalkmet(nslowness, rho, ogreen, slo, ntraces)
c
c scale coefficient matrix correct for K_0 damping
c
include 'greda_dim.inc'
include 'greda.inc'
include 'greda_pwo.inc'
c
integer nslo, ntr
integer nslowness, ntraces
real rho
complex ogreen(maxslo)
real slo(maxslo)
c
integer iu, it
double precision tf_gsl_sf_bessel_k1
double precision tf_gsl_sf_bessel_k0, factor
c
c scale to reflectivity inner product
do iu=1,nslo
ogreen(iu)=ogreen(iu)*tf_gsl_sf_bessel_k0(dble(slo(iu)*rho))
enddo
if (radial) then
do iu=1,nslowness
ogreen(iu)=ogreen(iu)*tf_gsl_sf_bessel_k1(dble(slo(iu)*rho))
enddo
else
do iu=1,nslowness
ogreen(iu)=ogreen(iu)*tf_gsl_sf_bessel_k0(dble(slo(iu)*rho))
enddo
endif
c
c scale phasor walkout
if (pwo_nlu.gt.0) then
do iu=1,pwo_nlu
factor=tf_gsl_sf_bessel_k0(dble(pwo_p(pwo_lu(iu))*rho))
do it=1,ntr
pwo_phasor(iu,it)=pwo_phasor(iu,it)*factor
if (radial) then
if (pwo_nlu.gt.0) then
do iu=1,pwo_nlu
factor=tf_gsl_sf_bessel_k1(dble(pwo_p(pwo_lu(iu))*rho))
do it=1,ntraces
pwo_phasor(iu,it)=pwo_phasor(iu,it)*factor
enddo
enddo
enddo
endif
else
if (pwo_nlu.gt.0) then
do iu=1,pwo_nlu
factor=tf_gsl_sf_bessel_k0(dble(pwo_p(pwo_lu(iu))*rho))
do it=1,ntraces
pwo_phasor(iu,it)=pwo_phasor(iu,it)*factor
enddo
enddo
endif
endif
cc
c
return
end
c
c======================================================================
c
c method independent subroutines to handle wavelength specific tapers
c
c----------------------------------------------------------------------
c
subroutine initwltaper(length,fraction)
......@@ -2182,7 +2294,11 @@ c
c
c======================================================================
c
c do a linear inversion with discrete calculated gram matrix
c
c Method: Linear inversion using a Gram matrix formulation and
c boxcar damping
c
c This requires a discretized integration to calculate the Gram matrix
c
c----------------------------------------------------------------------
c
......@@ -2204,7 +2320,7 @@ c
c calculate upper half and diagonal
do i=1,ntr
do j=i,ntr
gram(i,j)=gramdisint(maxslo, nslo, slo, r(i), r(j), om)
gram(i,j)=gramdisint(maxslo, nslo, slo, r(i), r(j), om, radial)
enddo
enddo
c
......@@ -2223,7 +2339,9 @@ c
c calculate discrete gram integral
c
double precision function gramdisint(maxslo, nslo,
& slo, rj, rk, om)
& slo, rj, rk, om, radialcomp)
c
logical radialcomp
c
integer maxslo, nslo
real slo(maxslo), om, rj, rk
......@@ -2231,18 +2349,29 @@ c
integer iu
double precision du, result, pdu, argj, argk, pargj, pargk
double precision tf_dj0
double precision tf_dj1
c
pdu=slo(2)-slo(1)
pargj=om*rj
pargk=om*rk
result=0.d0
do iu=1,nslo
du=pdu
if ((iu.eq.1).or.(iu.eq.nslo)) du=pdu*0.5d0
argj=pargj*slo(iu)
argk=pargk*slo(iu)
result=result+tf_dj0(argj)*tf_dj0(argk)*slo(iu)*du
enddo
if (radialcomp) then
do iu=1,nslo
du=pdu
if ((iu.eq.1).or.(iu.eq.nslo)) du=pdu*0.5d0
argj=pargj*slo(iu)
argk=pargk*slo(iu)
result=result+tf_dj1(argj)*tf_dj1(argk)*slo(iu)*du
enddo
else
do iu=1,nslo
du=pdu
if ((iu.eq.1).or.(iu.eq.nslo)) du=pdu*0.5d0
argj=pargj*slo(iu)
argk=pargk*slo(iu)
result=result+tf_dj0(argj)*tf_dj0(argk)*slo(iu)*du
enddo
endif
gramdisint=result
c
return
......@@ -2250,7 +2379,10 @@ c
c
c======================================================================
c
c do a linear inversion with Parkers hint
c Method: Linear inversion using a Gram matrix formulation and
c Lorentz damping
c
c This is the HOP-method published by Henry, Orcutt and Parker.
c
c----------------------------------------------------------------------
c
......@@ -2274,6 +2406,8 @@ c
integer i,j
double precision tf_gsl_sf_bessel_k0
double precision tf_gsl_sf_bessel_i0
double precision tf_gsl_sf_bessel_k1
double precision tf_gsl_sf_bessel_i1
double precision arg
c
double precision q0,q1,pi,maxval,alpa,alpaphi
......@@ -2283,14 +2417,25 @@ c
c
i=firstrec
j=1
do while (i.gt.0)
arg=r(i)*rho
p(j)=tf_gsl_sf_bessel_i0(arg)
q(j)=tf_gsl_sf_bessel_k0(arg)
tidx(j)=i
i=chain(i)
j=j+1
enddo
if (radial) then
do while (i.gt.0)
arg=r(i)*rho
p(j)=tf_gsl_sf_bessel_i1(arg)
q(j)=tf_gsl_sf_bessel_k1(arg)
tidx(j)=i
i=chain(i)
j=j+1
enddo
else
do while (i.gt.0)
arg=r(i)*rho
p(j)=tf_gsl_sf_bessel_i0(arg)
q(j)=tf_gsl_sf_bessel_k0(arg)
tidx(j)=i
i=chain(i)
j=j+1
enddo
endif
c
c the array gram is used here to hold the inverse of the gram matrix!
if (ntr.lt.3) stop 'ERROR less than 3 traces is senseless'
......
......@@ -25,6 +25,7 @@ c REVISIONS and CHANGES
c 23/10/97 V1.0 Thomas Forbriger
c 12/06/02 V1.1 supply wavelength specific taper factors
c 17/09/02 V1.2 move dimensions to greda_dim.inc
c 10/01/12 V1.3 provide component selection
c
c==============================================================================
c
......@@ -50,6 +51,9 @@ c wavelength taper length in units of wavelength
c
c wavelength taper falling edge start in units of the taper length
real wltaperfrac
c
c true, if radial compoment has to be transformed
logical radial
c
c----------------------------------------------------------------------
c
......@@ -57,5 +61,6 @@ c common blocks
common /hilbert/alpha, gram
common /seismo/spectra, nsamp, nom, ntr, nslo
common /wltaper/wltaperfac, wltaperlen, wltaperfrac
common /component/ radial
c
c ----- END OF greda.inc -----
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