Commit a2140999 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
SUBDIRS = src python
EXTRA_DIST = doc
QSSP
1) Compile and install QSSP
> autoreconf -i # only if 'configure' script is missing
> ./configure
> make
> sudo make install
AC_INIT([qssp], [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
This diff is collapsed.
wrapper_PYTHON = qssp.py
wrapperdir = $(pythondir)/pyrocko/fomosto
This diff is collapsed.
bin_PROGRAMS = qssp
qssp_SOURCES = brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f four1.f getdata.f moments.f pemain.f qpdifmat2.f qpdifmat4.f qpdifmat6.f qpdlegendre.f qpfftinv.f qpgetinp.f qpgrnspec.f qplegendre.f qpmain.f qppsvkern.f qppsvkerng.f qpqmodel.f qpshkern.f qpsmat0.f qpsmatc.f qpsmat.f qpsprop0.f qpsprop1.f qpsprop.f qpspropg0.f qpspropg1.f qpspropg.f qpstart0.f qpstart4.f qpstart6.f qpsublayer.f qptmat.f qptprop1.f qptprop.f qpwvint.f ruku.f spbdphj.f spbdphy.f spbjh.f spbphj.f spbphy.f spbpsj.f spbpsy.f taper.f qpglobal.h
subroutine brunewvlet(tau,df,mm2,wvf)
implicit none
integer mm2
double precision tau,df,fi
double complex wvf(mm2)
c
integer l
double precision f
c
double precision pi2
double complex c1
data pi2/6.28318530717959d0/
data c1/(1.d0,0.d0)/
c
c for wavelet: Brune's source time function
c
do l=1,mm2
f=df*dble(l-1)
wvf(l)=c1/dcmplx(1.d0,pi2*f*tau)**2
enddo
return
end
\ No newline at end of file
subroutine butterworth(n,fc,df,nf,lpf)
implicit none
integer n,nf
double precision fc,df
double complex lpf(nf)
c
integer l,k
double complex s,cs
c
double precision PI
data PI/3.14159265358979d0/
c
if(n.le.0.or.fc.le.0.d0)then
do l=1,nf
lpf(l)=(1.d0,0.d0)
enddo
else if(n.eq.2*(n/2))then
do l=1,nf
s=dcmplx(0.d0,df*dble(l-1)/fc)
lpf(l)=(1.d0,0.d0)
do k=1,n/2
cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
enddo
enddo
else
do l=1,nf
s=dcmplx(0.d0,df*dble(l-1)/fc)
lpf(l)=(1.d0,0.d0)/(s+(1.d0,0.d0))
do k=1,(n-1)/2
cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
enddo
enddo
endif
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)=dcmplx(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 cdsvd500(ca,cb,n,m,eps,key)
implicit none
c-------------------------------------------------------------------------------
c Solve complex linear equation system by single-value decomposiztion i
c method (modified from ludcmp and lubksb in the <Numerical Recipies> i
c ca: coefficient matrix(n,n); i
c cb: right-hand matrix(n,m) by input, i
c solution matrix(n,m) by return; i
c cunit: unit of the culomn vectors 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 i
c Note: n <= 500 will be NOT CHECKED! i
c-------------------------------------------------------------------------------
integer n,m,key
double precision eps
double complex ca(n,n),cb(n,m)
c
integer NMAX
parameter (NMAX=500)
integer i,ii,imax,j,k,ll
integer indx(NMAX)
double precision aamax,dum
double precision vv(NMAX)
double complex cdum,csum
c
do i=1,n
aamax=0.d0
do j=1,n
aamax=dmax1(aamax,cdabs(ca(i,j)))
enddo
if(aamax.le.eps)then
key=0
return
endif
vv(i)=1.d0/aamax
enddo
do j=1,n
do i=1,j-1
csum=ca(i,j)
do k=1,i-1
csum=csum-ca(i,k)*ca(k,j)
enddo
ca(i,j)=csum
enddo
aamax=0.d0
do i=j,n
csum=ca(i,j)
do k=1,j-1
csum=csum-ca(i,k)*ca(k,j)
enddo
ca(i,j)=csum
dum=vv(i)*cdabs(csum)
if(dum.ge.aamax) then
imax=i
aamax=dum
endif
enddo
if(j.ne.imax) then
do k=1,n
cdum=ca(imax,k)
ca(imax,k)=ca(j,k)
ca(j,k)=cdum
enddo
vv(imax)=vv(j)
endif
indx(j)=imax
if(cdabs(ca(j,j)).le.eps)then
key=0
return
endif
if(j.ne.n) then
cdum=(1.d0,0.d0)/ca(j,j)
do i=j+1,n
ca(i,j)=ca(i,j)*cdum
enddo
endif
enddo
c
do k=1,m
ii=0
do i=1,n
ll=indx(i)
csum=cb(ll,k)
cb(ll,k)=cb(i,k)
if(ii.ne.0) then
do j=ii,i-1
csum=csum-ca(i,j)*cb(j,k)
enddo
else if(cdabs(csum).ne.0.d0)then
ii=i
endif
cb(i,k)=csum
enddo
do i=n,1,-1
csum=cb(i,k)
do j=i+1,n
csum=csum-ca(i,j)*cb(j,k)
enddo
cb(i,k)=csum/ca(i,i)
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 disazi(rearth,lateq,loneq,latst,lonst,xnorth,yeast)
implicit none
c
double precision rearth,lateq,loneq,latst,lonst,xnorth,yeast
c
integer iangle
double precision latb,lonb,latc,lonc,angleb,anglec
double precision dis,a,b,c,s,aa
c
double precision PI,PI2,DEGTORAD
data PI,PI2/3.14159265358979d0,6.28318530717959d0/
data DEGTORAD/1.745329251994328d-02/
c
c A is north pole
c B is station
c C is epicentre (local cartesian coordinate origin)
c
latb=lateq*DEGTORAD
lonb=loneq*DEGTORAD
latc=latst*DEGTORAD
lonc=lonst*DEGTORAD
c
if(lonb.lt.0.d0)lonb=lonb+PI2
if(lonc.lt.0.d0)lonc=lonc+PI2
c
b=0.5d0*PI-latb
c=0.5d0*PI-latc
if(lonc.gt.lonb)then
aa=lonc-lonb
if(aa.le.PI)then
iangle=1
else
aa=PI2-aa
iangle=-1
endif
else
aa=lonb-lonc
if(aa.le.PI)then
iangle=-1
else
aa=PI2-aa
iangle=1
endif
endif
s=dcos(b)*dcos(c)+dsin(b)*dsin(c)*dcos(aa)
a=dacos(dsign(dmin1(dabs(s),1.d0),s))
dis=a*rearth
if(a*b*c.eq.0.d0)then
angleb=0.d0
anglec=0.d0
else
s=0.5d0*(a+b+c)
a=dmin1(a,s)
b=dmin1(b,s)
c=dmin1(c,s)
anglec=2.d0*dasin(dmin1(1.d0,
& dsqrt(dsin(s-a)*dsin(s-b)/(dsin(a)*dsin(b)))))
angleb=2.d0*dasin(dmin1(1.d0,
& dsqrt(dsin(s-a)*dsin(s-c)/(dsin(a)*dsin(c)))))
if(iangle.eq.1)then
angleb=PI2-angleb
else
anglec=PI2-anglec
endif
endif
c
c cartesian coordinates of the station
c
xnorth=dis*dcos(anglec)
yeast=dis*dsin(anglec)
c
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 moments(munit,strike,dip,rake,
& mtt,mpp,mrr,mtp,mpr,mrt)
implicit none
double precision munit,strike,dip,rake,
& mtt,mpp,mrr,mtp,mpr,mrt
c
double precision DEG2RAD
parameter(DEG2RAD=1.745329251994328d-02)
c
double precision st,di,ra
double precision sss,sss2,ss2s,css,css2,cs2s
double precision ssd,ss2d,csd,cs2d,ssr,csr
c
st=strike*DEG2RAD
di=dip*DEG2RAD
ra=rake*DEG2RAD
c
sss=dsin(st)
sss2=sss*sss
ss2s=dsin(2.d0*st)
css=dcos(st)
css2=css*css
cs2s=dcos(2.d0*st)
c
ssd=dsin(di)
ss2d=dsin(2.d0*di)
csd=dcos(di)
cs2d=dcos(2.d0*di)
c
ssr=dsin(ra)
csr=dcos(ra)
c
mtt=munit*(-ssd*csr*ss2s-ss2d*ssr*sss2)
mpp=munit*(ssd*csr*ss2s-ss2d*ssr*css2)
mrr=-(mtt+mpp)
mtp=munit*(-ssd*csr*cs2s-0.5d0*ss2d*ssr*ss2s)
mpr=munit*(csd*csr*sss-cs2d*ssr*css)
mrt=munit*(-csd*csr*css-cs2d*ssr*sss)
c
return
end
\ No newline at end of file
subroutine qpdifmat2(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
double precision mass,ro1,dro
double complex cup,clw,crr,drr,crorr,cxirr,clarr,cmurr,cgrrr
c
double complex c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
crr=dcmplx(rr,0.d0)
cup=(crr-rrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
crorr=cup*croup(ly)+clw*crolw(ly)
clarr=cup*claup(ly)+clw*clalw(ly)
cmurr=cup*cmuup(ly)+clw*cmulw(ly)
cxirr=clarr+c2*cmurr
c
dro=(roup(ly)-rolw(ly))/(rrup(ly)-rrlw(ly))
ro1=rolw(ly)-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))
& *((4.d0/3.d0)*ro1*(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)+rr*rrlw(ly)**2+rrlw(ly)**3))
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/cxirr/crr
c
mat(2,1)=c4*(cmurr*(c1+c2*clarr/cxirr)/crr-crorr*cgrrr)
& -crorr*comi2*crr
mat(2,2)=c2*clarr/cxirr/crr
c
return
end
\ No newline at end of file
subroutine qpdifmat4(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
double precision alfa
double complex cldeg,clp1,cll1
double complex crr,crorr,clarr,cgrrr,cgarr
c
double complex c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
cldeg=dcmplx(dble(ldeg),0.d0)
clp1=cldeg+c1
cll1=cldeg*clp1
c
crr=dcmplx(rr,0.d0)
c
alfa=(rr-rrlw(ly))/(rrup(ly)-rrlw(ly))
crorr=dcmplx(rolw(ly)
& *dexp((dlog(roup(ly))-dlog(rolw(ly)))*alfa),0.d0)
clarr=clalw(ly)+(claup(ly)-clalw(ly))*dcmplx(alfa,0.d0)
cgarr=cgalw(ly)+(cgaup(ly)-cgalw(ly))*dcmplx(alfa,0.d0)
cgrrr=cgrlw(ly)+(cgrup(ly)-cgrlw(ly))*dcmplx(alfa,0.d0)
c
mat(1,1)=crorr*cgrrr/clarr-c1/crr
mat(1,2)=cll1/crr-crorr*comi2*crr/clarr
mat(1,3)=-crorr*crr/clarr
mat(1,4)=(0.d0,0.d0)
c
mat(2,1)=c1/crr
mat(2,2)=(0.d0,0.d0)
mat(2,3)=(0.d0,0.d0)
mat(2,4)=(0.d0,0.d0)
c
mat(3,1)=cgarr/crr
mat(3,2)=(0.d0,0.d0)
mat(3,3)=-clp1/crr
mat(3,4)=c1/crr
c
mat(4,1)=cgarr*clp1/crr
mat(4,2)=-cgarr*cll1/crr
mat(4,3)=(0.d0,0.d0)
mat(4,4)=cldeg/crr
return
end
\ No newline at end of file
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