-
Notifications
You must be signed in to change notification settings - Fork 1
/
qsve2am.f
51 lines (51 loc) · 1.36 KB
/
qsve2am.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
subroutine qsve2am(n,ck,y,c,nj,rsite)
implicit none
c
c converting displacement-stress vectors to wave amplitudes
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 ca,cadp,cads,swap(4)
c
if(rsite)then
ca=(0.5d0,0.d0)/accrs(n)
cadp=ca/kprs(n)
cads=ca/ksrs(n)
do j=1,nj
swap(1)=(-wars(n)*y(1,j)+ck*y(4,j))*cadp
swap(2)=(-y(2,j)+wbrs(n)*y(3,j))*ca
swap(3)=(-ck*y(2,j)+wars(n)*y(3,j))*cads
swap(4)=(wbrs(n)*y(1,j)-y(4,j))*ca
c
c(1,j)= swap(1)+swap(2)
c(2,j)=-swap(1)+swap(2)
c(3,j)=-swap(3)+swap(4)
c(4,j)= swap(3)+swap(4)
enddo
else
ca=(0.5d0,0.d0)/acc(n)
cadp=ca/kp(n)
cads=ca/ks(n)
do j=1,nj
swap(1)=(-wa(n)*y(1,j)+ck*y(4,j))*cadp
swap(2)=(-y(2,j)+wb(n)*y(3,j))*ca
swap(3)=(-ck*y(2,j)+wa(n)*y(3,j))*cads
swap(4)=(wb(n)*y(1,j)-y(4,j))*ca
c
c(1,j)= swap(1)+swap(2)
c(2,j)=-swap(1)+swap(2)
c(3,j)=-swap(3)+swap(4)
c(4,j)= swap(3)+swap(4)
enddo
endif
c
return
end