Commit 45157473 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

initial commit

parents
Makefile.in
aclocal.m4
autom4te.cache
configure
install-sh
missing
py-compile
Makefile
config.log
config.status
*.o
src/qseis
SUBDIRS = src python
QSEIS
1) Compile and install QSEIS
> autoreconf -i # only if 'configure' script is missing
> ./configure
> make
> sudo make install
AC_INIT([qseis], [1.0])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AM_PATH_PYTHON([2.5])
AC_PROG_F77
AC_CONFIG_FILES([
Makefile
src/Makefile
python/Makefile
])
AC_OUTPUT
wrapper_PYTHON = qseis.py
wrapperdir = $(pythondir)/pyrocko/fomosto
This diff is collapsed.
bin_PROGRAMS = qseis
qseis_SOURCES = bessj0.f bessj1.f bessj.f caxcb.f cdgemp.f cmemcpy.f four1.f getdata.f qsam2ve.f qsbsj.f qsfftinv.f qsgetinp.f qshksh.f qskern.f qslayer.f qsmain.f qsmultis.f qspsv.f qsqmodel.f qssh.f qssource.f qssublay.f qsve2am.f qswavelet.f qswaveno.f qswvint.f taper.f wavelet.f qsglobal.h
double precision function bessj(n,x)
implicit none
integer n
double precision x
c
integer iacc
parameter(iacc=40)
double precision bigno,bigni
parameter(bigno=1.d+10,bigni=1.d-10)
c
integer j,jsum,m
double precision ax,bj,bjm,bjp,sum,tox,bessj0,bessj1
if(n.lt.2)pause 'bad argument n in bessj'
ax=dabs(x)
if(ax.eq.0.d0)then
bessj=0.d0
else if(ax.gt.dble(n))then
tox=2.d0/ax
bjm=bessj0(ax)
bj=bessj1(ax)
do 11 j=1,n-1
bjp=dble(j)*tox*bj-bjm
bjm=bj
bj=bjp
11 continue
bessj=bj
else
tox=2.d0/ax
m=2*((n+idint(dsqrt(dble(iacc*n))))/2)
bessj=0.d0
jsum=0
sum=0.d0
bjp=0.d0
bj=1.d0
do 12 j=m,1,-1
bjm=dble(j)*tox*bj-bjp
bjp=bj
bj=bjm
if(dabs(bj).gt.bigno)then
bj=bj*bigni
bjp=bjp*bigni
bessj=bessj*bigni
sum=sum*bigni
endif
if(jsum.ne.0)sum=sum+bj
jsum=1-jsum
if(j.eq.n)bessj=bjp
12 continue
sum=2.d0*sum-bj
bessj=bessj/sum
endif
if(x.lt.0.d0.and.mod(n,2).eq.1)bessj=-bessj
return
end
double precision function bessj0(x)
implicit none
c
c J_0(x)
c
double precision x,ax,x1,x2,theta,fct
double precision a0,a1,a2,a3,a4,a5,a6
double precision b0,b1,b2,b3,b4,b5,b6
double precision c0,c1,c2,c3,c4,c5,c6
c
data a0,a1,a2,a3,a4,a5,a6/1.00000000d0,
+ -2.24999970d0, 1.26562080d0,-0.31638660d0,
+ 0.04444790d0,-0.00394440d0, 0.00021000d0/
data b0,b1,b2,b3,b4,b5,b6/0.79788456d0,
+ -0.00000077d0,-0.00552740d0,-0.00009512d0,
+ 0.00137237d0,-0.00072805d0, 0.00014476d0/
data c0,c1,c2,c3,c4,c5,c6/-0.78539816d0,
+ -0.04166397d0,-0.00003954d0, 0.00262573d0,
+ -0.00054125d0,-0.00029333d0, 0.00013558d0/
c
ax=dabs(x)
if(ax.le.3.d0)then
x2=(ax/3.d0)**2
bessj0=a0+x2*(a1+x2*(a2+x2*(a3+x2*(a4+x2*(a5+x2*a6)))))
else
x1=3.d0/ax
fct=b0+x1*(b1+x1*(b2+x1*(b3+x1*(b4+x1*(b5+x1*b6)))))
theta=ax+c0+x1*(c1+x1*(c2+x1*(c3+x1*(c4+x1*(c5+x1*c6)))))
bessj0=fct*dcos(theta)/dsqrt(ax)
endif
c
return
end
double precision function bessj1(x)
implicit none
c
c J_1(x)
c
double precision x,ax,x1,x2,theta,fct
double precision a0,a1,a2,a3,a4,a5,a6
double precision b0,b1,b2,b3,b4,b5,b6
double precision c0,c1,c2,c3,c4,c5,c6
c
data a0,a1,a2,a3,a4,a5,a6/ 0.50000000d0,
+ -0.56249985d0, 0.21093573d0,-0.03954289d0,
+ 0.00443319d0,-0.00031761d0, 0.00001109d0/
data b0,b1,b2,b3,b4,b5,b6/ 0.79788456d0,
+ 0.00000156d0, 0.01659667d0, 0.00017105d0,
+ -0.00249511d0, 0.00113653d0,-0.00020033d0/
data c0,c1,c2,c3,c4,c5,c6/ -2.35619449d0,
+ 0.12499612d0, 0.00005650d0,-0.00637879d0,
+ 0.00074348d0, 0.00079824d0,-0.00029166d0/
c
ax=dabs(x)
if(ax.le.3.d0)then
x2=(ax/3.d0)**2
bessj1=x*(a0+x2*(a1+x2*(a2+x2*(a3+x2*(a4+x2*(a5+x2*a6))))))
else
x1=3.d0/ax
fct=b0+x1*(b1+x1*(b2+x1*(b3+x1*(b4+x1*(b5+x1*b6)))))
theta=ax+c0+x1*(c1+x1*(c2+x1*(c3+x1*(c4+x1*(c5+x1*c6)))))
bessj1=dsign(1.d0,x)*fct*dcos(theta)/dsqrt(ax)
endif
c
return
end
subroutine caxcb(a,b,n,l,m,c)
implicit none
c
c calculate c=a*b
c
integer n,l,m
double complex a(n,l),b(l,m),c(n,m)
c
integer i,j,k
c
do j=1,m
do i=1,n
c(i,j)=(0.d0,0.d0)
do k=1,l
c(i,j)=c(i,j)+a(i,k)*b(k,j)
enddo
enddo
enddo
return
end
subroutine cdgemp(a,b,n,n1,eps,key)
implicit none
c-------------------------------------------------------------------------------
c subroutine gemp to solve linear equation system i
c a: coefficient matrix(n,n); i
c b: right-hand matrix(n,n1) by input, i
c solution matrix(n,n1) by return; i
c eps: control constant; i
c key: if the main term of a column is i
c smaller than eps, key=0: anormal return, i
c else key=1: normal return. i
c-------------------------------------------------------------------------------
integer n,n1,key
double precision eps
double complex a(n,n),b(n,n1)
c
integer i,j,k,l,m
double precision p
double complex q
c
do m=1,n
p=0.d0
do i=m,n
do j=1,n
if(cdabs(a(i,j)).le.p)goto 10
p=cdabs(a(i,j))
k=i
l=j
10 continue
enddo
enddo
key=0
if(p.le.eps)return
do j=1,n
q=a(k,j)
a(k,j)=a(m,j)
a(m,j)=q
enddo
do j=1,n1
q=b(k,j)
b(k,j)=b(m,j)
b(m,j)=q
enddo
do i=1,n
q=-a(i,l)/a(m,l)
if(i.eq.m)goto 20
do j=1,n
a(i,j)=a(i,j)+a(m,j)*q
enddo
do j=1,n1
b(i,j)=b(i,j)+b(m,j)*q
enddo
20 continue
enddo
q=a(m,l)
do j=1,n
a(m,j)=a(m,j)/q
enddo
do j=1,n1
b(m,j)=b(m,j)/q
enddo
enddo
do j=1,n
do i=1,n
if(cdabs(a(i,j)).lt.0.5d0)goto 30
do k=1,n
q=a(i,k)
a(i,k)=a(j,k)
a(j,k)=q
enddo
do k=1,n1
q=b(i,k)
b(i,k)=b(j,k)
b(j,k)=q
enddo
30 continue
enddo
enddo
key=1
return
end
subroutine cmemcpy(a,b,n)
implicit none
c
c copy complex array a to b
c
integer n
double complex a(n),b(n)
c
integer i
c
do i=1,n
b(i)=a(i)
enddo
return
end
subroutine four1(data,nn,isign)
implicit none
integer nn,isign
double precision data(2*nn)
c
c fast fourier transform (fft)
c convention: f(t)=\int F(f)e^{-i2\pi ft} df
c replace data by its discrete fourier transform, if isign is input
c as 1; or replace data by nn times its inverse discrete fourier
c transform, if isign is input as -1. data is a double complex array of
c length nn or, equivalently, a double precision array of length 2*nn.
c nn must be an integer power of 2 (this is not checked!)
c
c note for convention: f(t)=\int F(f)e^{i2\pi f t} df, t-domain ==>
c f-domain, if isign=-1, and f-domain ==> t-domain, if isign=1.
c
integer i,j,n,m,mmax,istep
double precision tempr,tempi
double precision wr,wi,wpr,wpi,wtemp,theta
n=2*nn
j=1
do i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if((m.ge.2).and.(j.gt.m))then
j=j-m
m=m/2
goto 1
endif
j=j+m
enddo
mmax=2
2 if(n.gt.mmax)then
istep=2*mmax
theta=6.28318530717959d0/dble(isign*mmax)
wpr=-2.d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.d0
wi=0.d0
do m=1,mmax,2
do i=m,n,istep
j=i+mmax
tempr=wr*data(j)-wi*data(j+1)
tempi=wr*data(j+1)+wi*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
enddo
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
enddo
mmax=istep
goto 2
endif
return
end
subroutine getdata(unit,line)
implicit none
integer unit
character line*180,char*1
c
integer i
c
c this subroutine reads over all comment lines starting with "#".
c
char='#'
100 continue
if(char.eq.'#')then
read(unit,'(a)')line
i=1
char=line(1:1)
200 continue
if(char.eq.' ')then
i=i+1
char=line(i:i)
goto 200
endif
goto 100
endif
c
return
end
subroutine qsam2ve(n,ck,y,c,nj,rsite)
implicit none
c
c converting wave amplitudes to displacement-stress vectors
c by Langer block-diagonal decomposition technique for com-
c putational efficiency
c
integer n,nj
double complex ck
double complex y(4,nj),c(4,nj)
logical rsite
c
include 'qsglobal.h'
c
integer j
double complex bp,bs,swap(4)
c
if(rsite)then
bp=wbrs(n)*kprs(n)
bs=wbrs(n)*ksrs(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kprs(n)*swap(1)+ck*swap(4)
y(2,j)=wars(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ksrs(n)*swap(3)
y(4,j)=bp*swap(1)+wars(n)*swap(4)
enddo
else
bp=wb(n)*kp(n)
bs=wb(n)*ks(n)
do j=1,nj
swap(1)= c(1,j)-c(2,j)
swap(2)= c(1,j)+c(2,j)
swap(3)=-c(3,j)+c(4,j)
swap(4)= c(3,j)+c(4,j)
c
y(1,j)=kp(n)*swap(1)+ck*swap(4)
y(2,j)=wa(n)*swap(2)-bs*swap(3)
y(3,j)=ck*swap(2)-ks(n)*swap(3)
y(4,j)=bp*swap(1)+wa(n)*swap(4)
enddo
endif
c
return
end
subroutine qsbsj(dk,nk)
implicit none
integer nk
double precision dk
c
include 'qsglobal.h'
c
integer i,ir,ik
double precision k,x
double precision bessj0,bessj1,bessj
c
do ik=1,nk
k=dble(ik)*dk
do ir=1,nr
x=k*r(ir)
bsj(ik,0,ir)=bessj0(x)
bsj(ik,1,ir)=bessj1(x)
if(x.gt.2.d0)then
bsj(ik,2,ir)=bsj(ik,1,ir)*2.d0/x-bsj(ik,0,ir)
else
bsj(ik,2,ir)=bessj(2,x)
endif
if(x.gt.3.d0)then
bsj(ik,3,ir)=bsj(ik,2,ir)*4.d0/x-bsj(ik,1,ir)
else
bsj(ik,3,ir)=bessj(3,x)
endif
bsj(ik,-1,ir)=-bsj(ik,1,ir)
do i=-1,3
bsj(ik,i,ir)=bsj(ik,i,ir)*geospr(ir)
enddo
enddo
enddo
c
return
end
\ No newline at end of file
subroutine qsfftinv(icmp,istp)
implicit none
c
integer icmp,istp
c
include 'qsglobal.h'
c
integer nn,nn2,nch,nch2
parameter (nn=2*nfmax,nn2=nfmax,nch=2*nrmax,nch2=nrmax)
integer lf,mf,ir,j,it
double precision t,pi,pi2,slw0,omi
double precision f(nn),y0(nch),y(nch)
double complex s
double complex cy(nn,nch2)
c
double complex wvf(nn2)
save wvf
c
pi=4.d0*datan(1.d0)
pi2=2.d0*pi
c
if(v0.gt.0.d0)then
slw0=1.d0/v0
else
slw0=0.d0
endif
c
print *,' Output: '//outfile(icmp,istp)(1:flen(icmp,istp))
do lf=1,nf
f(lf)=dble(lf-1)*df
do ir=1,nr
cy(lf,ir)=grns(lf,icmp,ir,istp)
enddo
enddo
c
if(v0.gt.0.d0.or.dabs(tstart).gt.0.d0)then
c
c for time reduction
c
do lf=1,nf
do ir=1,nr
s=dcmplx(-fi,f(lf))
& *dcmplx(2.d0*pi*(tstart+r(ir)*slw0),0.d0)
cy(lf,ir)=cy(lf,ir)*cdexp(s)
enddo
enddo
endif
c
c seismometer filtering
c
if(asm.ne.1.d0)then
do lf=1,nf
do ir=1,nr
cy(lf,ir)=dcmplx(asm,0.d0)*cy(lf,ir)
enddo
enddo
endif
c
if(nroot+npole.gt.0)then
do lf=1,nf
s=dcmplx(-2.d0*pi*fi,2.d0*pi*f(lf))
do ir=1,nr
do j=1,nroot
cy(lf,ir)=cy(lf,ir)*(s-root(j))
enddo
do j=1,npole
cy(lf,ir)=cy(lf,ir)/(s-pole(j))
enddo
enddo
enddo
endif
c
c muliplication with wavelet spectrum
c
if(iexist.ne.1)then
call qswavelet(wvf,nf)
iexist=1
endif
do lf=1,nf
do ir=1,nr
cy(lf,ir)=cy(lf,ir)*wvf(lf)
enddo
enddo
c
mf=1
do lf=2*nf,nf+2,-1
mf=mf+1
do ir=1,nr
cy(lf,ir)=dconjg(cy(mf,ir))
enddo
enddo
do ir=1,nr
cy(nf+1,ir)=(0.d0,0.d0)
enddo
c
c convention for Fourier transform:
c f(t)=\int F(f) exp(i2\pi f t) df
c
do ir=1,nr
call four1(cy(1,ir),2*nf,+1)
enddo
c
open(20,file=outfile(icmp,istp),status='unknown')
write(20,'(a,$)')' T_sec '
do ir=1,nr-1
write(20,'(a4,a1,a1,2a4,$)')' ',