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

dos2unix line endings

parent 3fcb8532
This diff is collapsed.
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
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
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
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))
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
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
enddo
endif
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 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 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 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
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
......@@ -11,8 +11,8 @@ c
data PI,PI2/3.14159265358979d0,6.28318530717959d0/
data DEGTORAD/1.745329251994328d-02/
c
c A is north pole
c B is station
c A is north pole
c B is station
c C is epicentre (local cartesian coordinate origin)
c
latb=lateq*DEGTORAD
......
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
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
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
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
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
......
subroutine qpdifmat4(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
double precision alfa
double complex cldeg,clp1,cll1
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)
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
subroutine qpdifmat6(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
double complex cldeg,clp1,cll1,cup,clw
c