Line data Source code
1 : MODULE m_kernel2
2 : C.............................................................kernel2
3 : c solution of 4 coupled equations (logariphmic mesh) for (l,\mu)
4 : c 1->k1, 2->k2
5 : c d/dxP1=[-k1]P1+exp(x)[e-v+2mc^2]/c^2Q1+exp(x)B/c^2\sigma_z[-1,-1]Q1
6 : c d/dxQ1=[k1]Q1-exp(x)[e-v]P1+exp(x)B[\sigma_z[1,1]P1+\sigma_z[1,2]P2]
7 : c d/dxP2=[-k2]P2+exp(x)[e-v+2mc^2]/c^2Q2+exp(x)B/c^2\sigma_z[-1,-1]Q2
8 : c d/dxQ2=[k2]Q2-exp(x)[e-v]P2+exp(x)B[\sigma_z[2,1]P1+\sigma_z[2,2]P2]
9 : c notations:
10 : c Pi=exp(x)*gi
11 : c Qi=exp(x)*cfi
12 : c \gamma=sqrt(k^2-(Z/c)^2)
13 : c v=0.5[V(1/2)+V(-1/2)]
14 : c b=0.5[V(1/2)-V(-1/2)]
15 : c Rydberg units: in charge
16 : c Hartree units: com.
17 : c NSOL= 2 - 4 equations
18 : c ------------ a. shick KFA 1996
19 : CONTAINS
20 340 : SUBROUTINE kernel2(mrad,nsol,xmj,k1,k2,xx1,xx2,e,v,b,ri,dx,
21 340 : + nmatch,nstart,dp,dq,wp,wq)
22 :
23 :
24 : USE m_constants, ONLY : c_light
25 : USE m_diff
26 : IMPLICIT NONE
27 : C ..
28 : C .. Scalar Arguments ..
29 : INTEGER, INTENT (IN) :: mrad
30 : REAL dx,e,xmj
31 : INTEGER k1,k2,nmatch,nsol,nstart
32 : C ..
33 : C .. Array Arguments ..
34 : REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
35 : + wp(2,2,mrad),wq(2,2,mrad),xx1(4),xx2(4)
36 : C ..
37 : C .. Local Scalars ..
38 : REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
39 : + dp13,dp14,dp21,dp22,dp23,dp24,dq11,dq12,dq13,dq14,dq21,dq22,
40 : + dq23,dq24,dxd8,expdxh,gk,gk1,gk1k,gkk1,gmk,gmk1,p1c,p2c,
41 : + q1c,q2c,r,vh,vmn,vmnp1,vn,vnp1,xk,xk1,xm
42 : INTEGER i,iny,ir,j,jri,n,nrk
43 : C ..
44 : C .. Local Arrays ..
45 340 : REAL bm(mrad),p1(mrad),p1s(mrad),p2(mrad),p2s(mrad),q1(mrad),
46 340 : + q1s(mrad),q2(mrad),q2s(mrad),vm(mrad)
47 : C ..
48 : C .. Intrinsic Functions ..
49 : INTRINSIC exp,real,sqrt
50 : C ..
51 : C .. Data statements ..
52 : DATA abz/0.5/
53 : C ..
54 340 : cc = c_light(2.0)
55 : c
56 229160 : DO ir = 1,mrad
57 228820 : p1(ir) = 0.0
58 228820 : q1(ir) = 0.0
59 228820 : p1s(ir) = 0.0
60 228820 : q1s(ir) = 0.0
61 : c
62 228820 : p2(ir) = 0.0
63 228820 : q2(ir) = 0.0
64 228820 : p2s(ir) = 0.0
65 228820 : q2s(ir) = 0.0
66 : c
67 686800 : DO i = 1,2
68 1601740 : DO j = 1,2
69 : c
70 915280 : wp(i,j,ir) = 0.0
71 915280 : wq(i,j,ir) = 0.0
72 : c
73 915280 : dp(i,j,ir) = 0.0
74 1372920 : dq(i,j,ir) = 0.0
75 : END DO
76 : END DO
77 : END DO
78 : C
79 340 : csq = cc*cc
80 340 : expdxh = exp(dx/2.0)
81 340 : dxd8 = dx/8.0
82 340 : jri = nmatch
83 340 : nrk = jri
84 : c-derivatives of potential
85 : CALL diff(
86 : > mrad,v,dx,jri,
87 340 : < vm)
88 : CALL diff(
89 : > mrad,b,dx,jri,
90 340 : < bm)
91 : c- begin to solve coupled Dirac equation in "magnetic" field
92 340 : xk = real(k1)
93 340 : xk1 = real(k2)
94 340 : xm = xmj
95 340 : gk = -xm/ (xk+abz)
96 340 : gk1 = -xm/ (xk1+abz)
97 340 : gkk1 = -sqrt(1.0-gk*gk)
98 340 : gk1k = gkk1
99 340 : gmk = xm/ (xk-abz)
100 340 : gmk1 = xm/ (xk1-abz)
101 : c
102 : c NB: V=R**2*POTENTIAL
103 : c NB: B=R**2*FIELD
104 : c
105 1020 : DO 30 iny = 1,nsol
106 : C WAB02140
107 680 : r = ri(nstart)
108 680 : vn = v(nstart)
109 680 : vmn = vm(nstart)
110 680 : bn = b(nstart)
111 680 : bmn = bm(nstart)
112 680 : IF (iny.EQ.1) THEN
113 340 : p1(nstart) = xx1(1)
114 340 : q1(nstart) = xx1(2)
115 340 : p2(nstart) = xx1(3)
116 340 : q2(nstart) = xx1(4)
117 : ELSE
118 340 : p1(nstart) = xx2(1)
119 340 : q1(nstart) = xx2(2)
120 340 : p2(nstart) = xx2(3)
121 340 : q2(nstart) = xx2(4)
122 : END IF
123 : c---->f.po solution of Dirac eq.
124 680 : c1 = vn/r**2 - e
125 680 : c2 = 1.0 - c1/csq
126 : c C2=2.0-C1/CSQ !H
127 680 : c3 = bn/r/r
128 680 : p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
129 : q1s(nstart) = xk*q1(nstart) + r*
130 680 : + ((c1+c3*gk)*p1(nstart)+c3*gkk1*p2(nstart))
131 : C
132 680 : p2s(nstart) = -xk1*p2(nstart) + r* (c2+c3/csq*gmk1)*q2(nstart)
133 : q2s(nstart) = xk1*q2(nstart) +
134 680 : + r* ((c1+c3*gk1)*p2(nstart)+c3*gk1k*p1(nstart))
135 : C ******************************************************************
136 : C * *
137 : C * THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS *
138 : C * *
139 : C ******************************************************************
140 : C
141 680 : n = nstart
142 : c-> y1
143 198760 : 10 p1c = p1(n)
144 198760 : q1c = q1(n)
145 198760 : p2c = p2(n)
146 198760 : q2c = q2(n)
147 : c-> k1=f(x1,y1)
148 198760 : dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
149 198760 : dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
150 198760 : dp21 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
151 198760 : dq21 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
152 : c-> y11=y1+0.5*k1
153 198760 : p1c = p1c + 0.5*dp11
154 198760 : q1c = q1c + 0.5*dq11
155 198760 : p2c = p2c + 0.5*dp21
156 198760 : q2c = q2c + 0.5*dq21
157 : c-> x1'=x1+0.5h
158 198760 : r = r*expdxh
159 : c-> interpolation of V and B for x1'
160 198760 : vnp1 = v(n+1)
161 198760 : vmnp1 = vm(n+1)
162 198760 : bnp1 = b(n+1)
163 198760 : bmnp1 = bm(n+1)
164 198760 : vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
165 198760 : bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
166 : c-> k2=f(x1',y11)
167 198760 : c1 = vh/r/r - e
168 198760 : c2 = 1.0 - c1/csq
169 : c C2=2.0-C1/CSQ !H
170 198760 : c3 = bh/r/r
171 : C WAB02850
172 198760 : dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
173 198760 : dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
174 198760 : dp22 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
175 198760 : dq22 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
176 : c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
177 198760 : p1c = p1c + 0.5* (dp12-dp11)
178 198760 : q1c = q1c + 0.5* (dq12-dq11)
179 198760 : p2c = p2c + 0.5* (dp22-dp21)
180 198760 : q2c = q2c + 0.5* (dq22-dq21)
181 : c-> k3=f(x1',y12)
182 198760 : dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
183 198760 : dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
184 198760 : dp23 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
185 198760 : dq23 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
186 : c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
187 198760 : p1c = p1c + dp13 - 0.5*dp12
188 198760 : q1c = q1c - 0.5*dq12 + dq13
189 198760 : p2c = p2c + dp23 - 0.5*dp22
190 198760 : q2c = q2c - 0.5*dq22 + dq23
191 : c-> x2=x1+h
192 198760 : n = n + 1
193 198760 : r = ri(n)
194 : c-> k4=f(x2,y13)
195 198760 : c1 = vnp1/r/r - e
196 198760 : c2 = 1.0 - c1/csq
197 : c C2=2.0-C1/CSQ !H
198 198760 : c3 = bnp1/r/r
199 : c
200 198760 : dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
201 198760 : dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
202 198760 : dp24 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
203 198760 : dq24 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
204 : c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
205 198760 : p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
206 198760 : q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
207 198760 : p2(n) = p2(n-1) + (dp21+2.0* (dp22+dp23)+dp24)/6.0
208 198760 : q2(n) = q2(n-1) + (dq21+2.0* (dq22+dq23)+dq24)/6.0
209 : c-> f(x2,y2)
210 198760 : p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
211 198760 : q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n)+c3*gkk1*p2(n))
212 198760 : p2s(n) = -xk1*p2(n) + r* (c2+c3/csq*gmk1)*q2(n)
213 198760 : q2s(n) = xk1*q2(n) + r* ((c1+c3*gk1)*p2(n)+c3*gk1k*p1(n))
214 : c-> redefinition of V ind B (2->1)
215 198760 : vn = vnp1
216 198760 : bn = bnp1
217 198760 : vmn = vmnp1
218 198760 : bmn = bmnp1
219 : C WAB03340
220 198760 : IF (n-nrk.LT.0) GOTO 10
221 : c---->Milne's method off now!
222 : c
223 : c storage of wave functions and their derivatives.
224 200120 : DO ir = 1,jri
225 199440 : wp(1,iny,ir) = p1(ir)
226 199440 : wp(2,iny,ir) = p2(ir)
227 199440 : wq(1,iny,ir) = q1(ir)
228 199440 : wq(2,iny,ir) = q2(ir)
229 : c derivatives
230 199440 : dp(1,iny,ir) = p1s(ir)
231 199440 : dp(2,iny,ir) = p2s(ir)
232 199440 : dq(1,iny,ir) = q1s(ir)
233 200120 : dq(2,iny,ir) = q2s(ir)
234 : END DO
235 : C*********************************************************** WAB05210
236 340 : 30 CONTINUE
237 : c write(oUnit,*) wp
238 :
239 340 : END SUBROUTINE kernel2
240 : END MODULE m_kernel2
|
