Line data Source code
1 : MODULE m_kernel1
2 : C..............................................................kernel1
3 : c solution of 2 coupled equations (logariphmic mesh) for (l,\mu:|mu|>=l)
4 : c 1->k1
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]
7 : c notations:
8 : c Pi=exp(x)*gi
9 : c Qi=exp(x)*cfi
10 : c \gamma=sqrt(k^2-(Z/c)^2)
11 : c v=0.5[V(1/2)+V(-1/2)]
12 : c b=0.5[V(1/2)-V(-1/2)]
13 : c Rydberg units: in charge
14 : c hartree units: com.
15 : c NSOL= 1 - 2 equations
16 : c ------------ a. shick KFA 1996
17 : CONTAINS
18 324 : SUBROUTINE kernel1(mrad,xmj,kap1,xx1,e,v,b,ri,dx,nmatch,nstart,
19 324 : + dp,dq,wp,wq)
20 :
21 : c
22 : USE m_constants, ONLY : c_light
23 : USE m_diff
24 : IMPLICIT NONE
25 : C ..
26 : C .. Scalar Arguments ..
27 : INTEGER, INTENT (IN) :: mrad
28 : REAL dx,e,xmj
29 : INTEGER kap1,nmatch,nstart
30 : C ..
31 : C .. Array Arguments ..
32 : REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
33 : + wp(2,2,mrad),wq(2,2,mrad),xx1(4)
34 : C ..
35 : C .. Local Scalars ..
36 : REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
37 : + dp13,dp14,dq11,dq12,dq13,dq14,dxd8,expdxh,gk,gmk,p1c,
38 : + q1c,r,vh,vmn,vmnp1,vn,vnp1,xk,xm
39 : INTEGER i,ir,j,jri,n,nrk
40 : C ..
41 : C .. Local Arrays ..
42 324 : REAL bm(mrad),p1(mrad),p1s(mrad),q1(mrad),q1s(mrad),vm(mrad)
43 : C ..
44 : C .. Intrinsic Functions ..
45 : INTRINSIC exp,real
46 : C ..
47 : C .. Data statements ..
48 : DATA abz/0.5/
49 : C ..
50 324 : cc = c_light(2.0)
51 : c
52 218376 : DO ir = 1,mrad
53 218052 : p1(ir) = 0.0
54 218052 : q1(ir) = 0.0
55 218052 : p1s(ir) = 0.0
56 218052 : q1s(ir) = 0.0
57 : c
58 654480 : DO i = 1,2
59 1526364 : DO j = 1,2
60 : c
61 872208 : wp(i,j,ir) = 0.0
62 872208 : wq(i,j,ir) = 0.0
63 : c
64 872208 : dp(i,j,ir) = 0.0
65 1308312 : dq(i,j,ir) = 0.0
66 : END DO
67 : END DO
68 : END DO
69 : C
70 324 : csq = cc*cc
71 324 : expdxh = exp(dx/2.0)
72 324 : dxd8 = dx/8.0
73 324 : jri = nmatch
74 324 : nrk = jri
75 : CALL diff(
76 : > mrad,v,dx,jri,
77 324 : < vm)
78 : CALL diff(
79 : > mrad,b,dx,jri,
80 324 : < bm)
81 324 : xk = real(kap1)
82 324 : xm = xmj
83 324 : gk = -xm/ (xk+abz)
84 324 : gmk = xm/ (xk-abz)
85 : c
86 324 : r = ri(nstart)
87 324 : vn = v(nstart)
88 324 : vmn = vm(nstart)
89 324 : bn = b(nstart)
90 324 : bmn = bm(nstart)
91 324 : p1(nstart) = xx1(1)
92 324 : q1(nstart) = xx1(2)
93 : c---->f.po solution of Dirac eq.
94 324 : c1 = vn/r**2 - e
95 324 : c2 = 1.0 - c1/csq
96 : c C2=2.0-C1/CSQ !H
97 324 : c3 = bn/r/r
98 324 : p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
99 324 : q1s(nstart) = xk*q1(nstart) + r* ((c1+c3*gk)*p1(nstart))
100 : C * THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS * WAB02520
101 324 : n = nstart
102 : c-> y1
103 91392 : 10 p1c = p1(n)
104 91392 : q1c = q1(n)
105 : c-> k1=f(x1,y1)
106 91392 : dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
107 91392 : dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
108 : c-> y11=y1+0.5*k1
109 91392 : p1c = p1c + 0.5*dp11
110 91392 : q1c = q1c + 0.5*dq11
111 : c-> x1'=x1+0.5h
112 91392 : r = r*expdxh
113 : c-> interpolation of V and B for x1'
114 91392 : vnp1 = v(n+1)
115 91392 : vmnp1 = vm(n+1)
116 91392 : bnp1 = b(n+1)
117 91392 : bmnp1 = bm(n+1)
118 91392 : vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
119 91392 : bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
120 : c-> k2=f(x1',y11)
121 91392 : c1 = vh/r/r - e
122 91392 : c2 = 1.0 - c1/csq
123 : c C2=2.0-C1/CSQ !H
124 91392 : c3 = bh/r/r
125 : C WAB02850
126 91392 : dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
127 91392 : dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
128 : c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
129 91392 : p1c = p1c + 0.5* (dp12-dp11)
130 91392 : q1c = q1c + 0.5* (dq12-dq11)
131 : c-> k3=f(x1',y12)
132 91392 : dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
133 91392 : dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
134 : c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
135 91392 : p1c = p1c + dp13 - 0.5*dp12
136 91392 : q1c = q1c - 0.5*dq12 + dq13
137 : c-> x2=x1+h
138 91392 : n = n + 1
139 91392 : r = ri(n)
140 : c-> k4=f(x2,y13)
141 91392 : c1 = vnp1/r/r - e
142 91392 : c2 = 1.0 - c1/csq
143 : c C2=2.0-C1/CSQ !H
144 91392 : c3 = bnp1/r/r
145 : c
146 91392 : dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
147 91392 : dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
148 : c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
149 91392 : p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
150 91392 : q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
151 : c-> f(x2,y2)
152 91392 : p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
153 91392 : q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n))
154 : c-> redefinition of V ind B (2->1)
155 91392 : vn = vnp1
156 91392 : bn = bnp1
157 91392 : vmn = vmnp1
158 91392 : bmn = bmnp1
159 : C WAB03340
160 91392 : IF (n-nrk.LT.0) GOTO 10
161 : c---->Milne's method off now!
162 218376 : DO ir = 1,mrad
163 218052 : wp(1,1,ir) = p1(ir)
164 218052 : wq(1,1,ir) = q1(ir)
165 : c
166 218052 : dp(1,1,ir) = p1s(ir)
167 218376 : dq(1,1,ir) = q1s(ir)
168 : END DO
169 :
170 324 : END SUBROUTINE kernel1
171 : END MODULE m_kernel1
|
