Line data Source code
1 332 : SUBROUTINE inconz(e,l,xmj,kap1,kap2,vv,bb,rr,xx1,xx2)
2 :
3 : c..........................................................inconz
4 : c initial point for outcome regular solution of dirac eq.
5 : c order kap1=-L-1, kap2=L
6 : c
7 : USE m_constants, ONLY : c_light
8 : IMPLICIT NONE
9 : C .. Scalar Arguments ..
10 : REAL bb,e,rr,vv,xmj
11 : INTEGER kap1,kap2,l
12 : C ..
13 : C .. Array Arguments ..
14 : REAL xx1(4),xx2(4)
15 : C ..
16 : C .. Local Scalars ..
17 : REAL aa11,aa12,aa21,aa22,bb1,bb2,bc0,bqq,cc,cg1,cg2,cg4,cg5,cg8,
18 : + cgo,csq,det,emvpp,emvqq,rpwgpm,tz,vc0
19 : INTEGER i,j,m,mps,nsol
20 : C ..
21 : C .. Local Arrays ..
22 : REAL cgd(2),cgmd(2),gam(2),kap(2),pc(2,2,0:1),qc(2,2,0:1),wp(2,2),
23 : + wq(2,2)
24 : C ..
25 : C .. Intrinsic Functions ..
26 : INTRINSIC abs,real,int,sqrt
27 : C ..
28 332 : cc = c_light(2.0)
29 332 : csq = cc*cc
30 : c
31 : C EXPANSION COEFFICIENTS FOR THE POTENTIAL AND B-FIELD
32 : C VV=VV(1)
33 332 : tz = real(int(-vv*rr))
34 332 : vc0 = vv - (-tz)/rr
35 : C BB=BB(1)
36 332 : bc0 = bb
37 : C
38 : C CALCULATE G-COEFFICIENTS OF B-FIELD
39 : C
40 : c KAP1 = - L - 1
41 : c KAP2 = + L
42 332 : cg1 = -xmj/ (kap1+0.5)
43 332 : cg5 = -xmj/ (-kap1+0.5)
44 332 : cgd(1) = cg1
45 332 : cgmd(1) = cg5
46 332 : kap(1) = real(kap1)
47 332 : gam(1) = sqrt(kap(1)**2- (tz/cc)**2)
48 332 : IF (abs(xmj).GE.l) THEN
49 162 : cg2 = 0.00
50 162 : cg4 = 0.00
51 162 : cg8 = 0.00
52 162 : nsol = 1
53 162 : cgd(2) = 0.00
54 162 : cgo = 0.00
55 162 : cgmd(2) = 0.00
56 162 : gam(2) = 0.00
57 162 : kap(2) = 0.00
58 : ELSE
59 170 : cg2 = -sqrt(1.0- (xmj/ (kap1+0.50))**2)
60 170 : cg4 = -xmj/ (kap2+0.50)
61 170 : cg8 = -xmj/ (-kap2+0.50)
62 170 : nsol = 2
63 170 : cgd(2) = cg4
64 170 : cgo = cg2
65 170 : cgmd(2) = cg8
66 170 : kap(2) = real(kap2)
67 170 : gam(2) = sqrt(kap(2)**2- (tz/cc)**2)
68 : END IF
69 : C
70 834 : DO 10 j = 1,nsol
71 502 : i = 3 - j
72 502 : pc(j,j,0) = sqrt(abs(kap(j)-gam(j)))
73 502 : qc(j,j,0) = (kap(j)+gam(j))* (csq/tz)*pc(j,j,0)
74 502 : pc(i,j,0) = 0.0
75 502 : qc(i,j,0) = 0.0
76 332 : 10 CONTINUE
77 : C DETERMINE HIGHER EXPANSION COEFFICIENTS FOR THE WAVE FUNCTIONS
78 332 : mps = 1
79 332 : aa12 = -tz/csq
80 332 : aa21 = tz
81 332 : emvqq = (e-vc0+csq)/csq
82 332 : emvpp = -e + vc0
83 332 : bqq = bc0/csq
84 834 : DO 40 j = 1,nsol
85 1004 : DO 30 m = 1,mps
86 1344 : DO 20 i = 1,nsol
87 842 : bb1 = (emvqq+bqq*cgmd(i))*qc(i,j,m-1)
88 : bb2 = (emvpp+bc0*cgd(i))*pc(i,j,m-1) +
89 842 : + bc0*cgo*pc(3-i,j,m-1)
90 842 : aa11 = gam(j) + m + kap(i)
91 842 : aa22 = gam(j) + m - kap(i)
92 842 : det = aa11*aa22 - aa12*aa21
93 842 : pc(i,j,m) = (bb1*aa22-aa12*bb2)/det
94 842 : qc(i,j,m) = (aa11*bb2-bb1*aa21)/det
95 502 : 20 CONTINUE
96 502 : 30 CONTINUE
97 332 : 40 CONTINUE
98 : C
99 : C PERFORM SUMMATION OVER WAVE FUNCTION - EXPANSION COEFFICIENTS
100 : C FOR THE FIRST - MESH - POINT
101 : c RR= RC(1)
102 834 : DO 80 j = 1,nsol
103 502 : rpwgpm = rr** (gam(j))
104 1344 : DO 50 i = 1,nsol
105 842 : wp(i,j) = pc(i,j,0)*rpwgpm
106 842 : wq(i,j) = qc(i,j,0)*rpwgpm
107 502 : 50 CONTINUE
108 1004 : DO 70 m = 1,mps
109 502 : rpwgpm = rpwgpm*rr
110 1344 : DO 60 i = 1,nsol
111 842 : wp(i,j) = wp(i,j) + pc(i,j,m)*rpwgpm
112 842 : wq(i,j) = wq(i,j) + qc(i,j,m)*rpwgpm
113 502 : 60 CONTINUE
114 502 : 70 CONTINUE
115 332 : 80 CONTINUE
116 : C---> First point solutions construction
117 332 : IF (nsol.EQ.2) THEN
118 170 : xx1(1) = wp(1,1)
119 170 : xx1(2) = wq(1,1)
120 170 : xx1(3) = wp(2,1)
121 170 : xx1(4) = wq(2,1)
122 : c
123 170 : xx2(1) = wp(1,2)
124 170 : xx2(2) = wq(1,2)
125 170 : xx2(3) = wp(2,2)
126 170 : xx2(4) = wq(2,2)
127 : ELSE
128 162 : xx1(1) = wp(1,1)
129 162 : xx1(2) = wq(1,1)
130 : c not needed
131 162 : xx1(3) = 0.0
132 162 : xx1(4) = 0.0
133 162 : xx2(1) = 0.0
134 162 : xx2(2) = 0.0
135 162 : xx2(3) = 0.0
136 162 : xx2(4) = 0.0
137 : END IF
138 332 : RETURN
139 : END
|
