]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MINICERN/mathlib/gen/c/wclbes.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / wclbes.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1996/04/01 15:01:57 mclareni
6* Mathlib gen
7*
8*
9#include "gen/pilot.h"
10#if defined(CERNLIB_DOUBLE)
11 SUBROUTINE WCLBES(ZZ,ETA1,ZLMIN,NL,FC,GC,FCP,GCP,SIG,KFN,MODE1,
12 1 IFAIL,IPR)
13C
14CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
15C C
16C COMPLEX COULOMB WAVEFUNCTION PROGRAM USING STEED'S METHOD C
17C Original title : COULCC C
18C C
19C A. R. Barnett Manchester March 1981 C
20C modified I.J. Thompson Daresbury, Sept. 1983 for Complex Functions C
21C C
22C The FORM (not the SUBSTANCE) of this program has been modified C
23C by K.S. KOLBIG (CERN) December 1987 C
24C C
25C original program RCWFN in CPC 8 (1974) 377-395 C
26C + RCWFF in CPC 11 (1976) 141-142 C
27C + COULFG in CPC 27 (1982) 147-166 C
28C description of real algorithm in CPC 21 (1981) 297-314 C
29C description of complex algorithm JCP 64 (1986) 490-509 C
30C this version written up in CPC 36 (1985) 363-372 C
31C C
32C WCLBES returns F,G,G',G',SIG for complex ETA, ZZ, and ZLMIN, C
33C for NL integer-spaced lambda values ZLMIN to ZLMIN+NL inclusive, C
34C thus giving complex-energy solutions to the Coulomb Schrodinger C
35C equation,to the Klein-Gordon equation and to suitable forms of C
36C the Dirac equation ,also spherical & cylindrical Bessel equations C
37C C
38C if ABS(MODE1) C
39C = 1 get F,G,F',G' for integer-spaced lambda values C
40C = 2 F,G unused arrays must be dimensioned in C
41C = 3 F, F' call to at least length (0:1) C
42C = 4 F C
43C = 11 get F,H+,F',H+' ) if KFN=0, H+ = G + i.F ) C
44C = 12 F,H+ ) >0, H+ = J + i.Y = H(1) ) in C
45C = 21 get F,H-,F',H-' ) if KFN=0, H- = G - i.F ) GC C
46C = 22 F,H- ) >0, H- = J - i.Y = H(2) ) C
47C C
48C if MODE1 < 0 then the values returned are scaled by an exponen- C
49C tial factor (dependent only on ZZ) to bring nearer C
50C unity the functions for large ABS(ZZ), small ETA , C
51C and ABS(ZL) < ABS(ZZ). C
52C Define SCALE = ( 0 if MODE1 > 0 C
53C ( IMAG(ZZ) if MODE1 < 0 & KFN < 3 C
54C ( REAL(ZZ) if MODE1 < 0 & KFN = 3 C
55C then FC = EXP(-ABS(SCALE)) * ( F, j, J, or I) C
56C and GC = EXP(-ABS(SCALE)) * ( G, y, or Y ) C
57C or EXP(SCALE) * ( H+, H(1), or K) C
58C or EXP(-SCALE) * ( H- or H(2) ) C
59C C
60C if KFN = 0,-1 complex Coulomb functions are returned F & G C
61C = 1 spherical Bessel " " " j & y C
62C = 2 cylindrical Bessel " " " J & Y C
63C = 3 modified cyl. Bessel " " " I & K C
64C C
65C and where Coulomb phase shifts put in SIG if KFN=0 (not -1) C
66C C
67C The use of MODE and KFN is independent C
68C (except that for KFN=3, H(1) & H(2) are not given) C
69C C
70C With negative orders lambda, WCLBES can still be used but with C
71C reduced accuracy as CF1 becomes unstable. The user is thus C
72C strongly advised to use reflection formulae based on C
73C H+-(ZL,,) = H+-(-ZL-1,,) * exp +-i(sig(ZL)-sig(-ZL-1)-(ZL+1/2)pi) C
74C C
75C Precision: results to within 2-3 decimals of 'machine accuracy', C
76C except in the following cases: C
77C (1) if CF1A fails because X too small or ETA too large C
78C the F solution is less accurate if it decreases with C
79C decreasing lambda (e.g. for lambda.LE.-1 & ETA.NE.0) C
80C (2) if ETA is large (e.g. >> 50) and X inside the C
81C turning point, then progressively less accuracy C
82C (3) if ZLMIN is around sqrt(ACCUR) distance from an C
83C integral order and abs(X) << 1, then errors present. C
84C RERR traces the main roundoff errors. C
85C C
86C WCLBES is coded for real*8 on IBM or equivalent ACCUR >= 10**-14 C
87C with a section of doubled REAL*16 for less roundoff errors. C
88C (If no doubled precision available, increase JMAX to eg 100)C
89C Use IMPLICIT COMPLEX*32 & REAL*16 on VS compiler ACCUR >= 10**-32 C
90C For single precision CDC (48 bits) reassign C
91C DOUBLE PRECISION=REAL etc. C
92C C
93C IPR on input = 0 : no printing of error messages C
94C ne 0 : print error messages on file 6 C
95C IFAIL in output = -2 : argument out of range C
96C = -1 : one of the continued fractions failed, C
97C or arithmetic check before final recursion C
98C = 0 : All Calculations satisfactory C
99C ge 0 : results available for orders up to & at C
100C position NL-IFAIL in the output arrays. C
101C = -3 : values at ZLMIN not found as over/underflowC
102C = -4 : roundoff errors make results meaningless C
103C C
104CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
105C C
106C Machine dependent constants : C
107C C
108C ACCU target bound on relative error (except near 0 crossings)C
109C (ACCUR should be at least 100 * ACC8) C
110C ACC8 smallest number with 1+ACC8 .ne.1 in REAL*8 arithmetic C
111C ACC16 smallest number with 1+ACC16.ne.1 in REAL*16 arithmetic C
112C FPMAX magnitude of largest floating point number * ACC8 C
113C FPMIN magnitude of smallest floating point number / ACC8 C
114C C
115C Parameters determining region of calculations : C
116C C
117C R20 estimate of (2F0 iterations)/(CF2 iterations) C
118C ASYM minimum X/(ETA**2+L) for CF1A to converge easily C
119C XNEAR minimum ABS(X) for CF2 to converge accurately C
120C LIMIT maximum no. iterations for CF1, CF2, and 1F1 series C
121C JMAX size of work arrays for Pade accelerations C
122C NDROP number of successive decrements to define instabilityC
123C C
124CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
125C
126C C309R1 = CF2, C309R2 = F11, C309R3 = F20,
127C C309R4 = CF1R, C309R5 = CF1C, C309R6 = CF1A,
128C C309R7 = RCF, C309R8 = CTIDY.
129C
130 IMPLICIT COMPLEX*16 (A-H,O-Z)
131C
132 COMMON /COULC2/ NFP,N11,NPQ1,NPQ2,N20,KAS1,KAS2
133 INTEGER NPQ(2),KAS(2)
134 EQUIVALENCE (NPQ(1),NPQ1),(NPQ(2),NPQ2)
135 EQUIVALENCE (KAS(1),KAS1),(KAS(2),KAS2)
136 DOUBLE PRECISION ZERO,ONE,TWO,HALF
137 DOUBLE PRECISION R20,ASYM,XNEAR
138 DOUBLE PRECISION FPMAX,FPMIN,FPLMAX,FPLMIN
139 DOUBLE PRECISION ACCU,ACC8,ACC16
140 DOUBLE PRECISION HPI,TLOG
141 DOUBLE PRECISION ERR,RERR,ABSC,ACCUR,ACCT,ACCH,ACCB,C309R4
142 DOUBLE PRECISION PACCQ,EPS,OFF,SCALE,SF,SFSH,TA,RK,OMEGA,ABSX
143 DOUBLE PRECISION DX1,DETA,DZLL
144
145 LOGICAL LPR,ETANE0,IFCP,RLEL,DONEM,UNSTAB,ZLNEG,AXIAL,NOCF2,NPINT
146
147 PARAMETER(ZERO = 0, ONE = 1, TWO = 2, HALF = ONE/TWO)
148 PARAMETER(CI = (0,1), CIH = HALF*CI)
149 PARAMETER(R20 = 3, ASYM = 3, XNEAR = HALF)
150 PARAMETER(LIMIT = 20000, NDROP = 5, JMAX = 50)
151
152#include "c309prec.inc"
153
154 PARAMETER(HPI = 1.57079 63267 94896 619D0)
155 PARAMETER(TLOG = 0.69314 71805 59945 309D0)
156
157 DIMENSION FC(0:*),GC(0:*),FCP(0:*),GCP(0:*),SIG(0:*)
158 DIMENSION XRCF(JMAX,4)
159
160#if defined(CERNLIB_QF2C)
161#include "defdr.inc"
162#endif
163 NINTC(W)=NINT(DREAL(W))
164 ABSC(W)=ABS(DREAL(W))+ABS(DIMAG(W))
165 NPINT(W,ACCB)=ABSC(NINTC(W)-W) .LT. ACCB .AND. DREAL(W) .LT. HALF
166C
167 MODE=MOD(ABS(MODE1),10)
168 IFCP=MOD(MODE,2) .EQ. 1
169 LPR=IPR .NE. 0
170 IFAIL=-2
171 N11=0
172 NFP=0
173 KAS(1)=0
174 KAS(2)=0
175 NPQ(1)=0
176 NPQ(2)=0
177 N20=0
178
179 ACCUR=MAX(ACCU,50*ACC8)
180 ACCT=HALF*ACCUR
181 ACCH=SQRT(ACCUR)
182 ACCB=SQRT(ACCH)
183 RERR=ACCT
184 FPLMAX=LOG(FPMAX)
185 FPLMIN=LOG(FPMIN)
186C
187 CIK=ONE
188 IF(KFN .GE. 3) CIK=CI*SIGN(ONE,FPMIN-DIMAG(ZZ))
189 X=ZZ*CIK
190 ETA=ETA1
191 IF(KFN .GT. 0) ETA=ZERO
192 ETANE0=ABSC(ETA) .GT. ACC8
193 ETAI=ETA*CI
194 DELL=ZERO
195 IF(KFN .GE. 2) DELL=HALF
196 ZM1=ZLMIN-DELL
197 SCALE=ZERO
198 IF(MODE1 .LT. 0) SCALE=DIMAG(X)
199C
200 M1=1
201 L1=M1+NL
202 RLEL=ABS(DIMAG(ETA))+ABS(DIMAG(ZM1)) .LT. ACC8
203 ABSX=ABS(X)
204 AXIAL=RLEL .AND. ABS(DIMAG(X)) .LT. ACC8*ABSX
205 IF(MODE .LE. 2 .AND. ABSX .LT. FPMIN) GO TO 310
206 XI=ONE/X
207 XLOG=LOG(X)
208
209C log with cut along the negative real axis, see also OMEGA
210
211 ID=1
212 DONEM=.FALSE.
213 UNSTAB=.FALSE.
214 LF=M1
215 IFAIL=-1
216 10 ZLM=ZM1+(LF-M1)
217 ZLL=ZM1+(L1-M1)
218C
219C *** ZLL is final lambda value, or 0.5 smaller for J,Y Bessels
220C
221 Z11=ZLL
222 IF(ID .LT. 0) Z11=ZLM
223 P11=CI*SIGN(ONE,ACC8-DIMAG(ETA))
224 LAST=L1
225C
226C *** Find phase shifts and Gamow factor at lambda = ZLL
227C
228 PK=ZLL+ONE
229 AA=PK-ETAI
230 AB=PK+ETAI
231 BB=TWO*PK
232 ZLNEG=NPINT(BB,ACCB)
233 CLGAA=WLOGAM(AA)
234 CLGAB=CLGAA
235 IF(ETANE0 .AND. .NOT.RLEL) CLGAB=WLOGAM(AB)
236 IF(ETANE0 .AND. RLEL) CLGAB=DCONJG(CLGAA)
237 SIGMA=(CLGAA-CLGAB)*CIH
238 IF(KFN .EQ. 0) SIG(L1)=SIGMA
239 IF(.NOT.ZLNEG) CLL=ZLL*TLOG-HPI*ETA-WLOGAM(BB)+(CLGAA+CLGAB)*HALF
240 THETA=X-ETA*(XLOG+TLOG)-ZLL*HPI+SIGMA
241C
242 TA=(DIMAG(AA)**2+DIMAG(AB)**2+ABS(DREAL(AA))+ABS(DREAL(AB)))*HALF
243 IF(ID .GT. 0 .AND. ABSX .LT. TA*ASYM .AND. .NOT.ZLNEG) GO TO 20
244C
245C *** use CF1 instead of CF1A, if predicted to converge faster,
246C (otherwise using CF1A as it treats negative lambda &
247C recurrence-unstable cases properly)
248C
249 RK=SIGN(ONE,DREAL(X)+ACC8)
250 P=THETA
251 IF(RK .LT. 0) P=-X+ETA*(LOG(-X)+TLOG)-ZLL*HPI-SIGMA
252 XRCF(1,1)=PK
253 F=RK*C309R6(X*RK,ETA*RK,ZLL,P,ACCT,JMAX,NFP,FEST,ERR,FPMAX,XRCF,
254 1 XRCF(1,3),XRCF(1,4))
255 FESL=LOG(FEST)+ABS(DIMAG(X))
256 NFP=-NFP
257 IF(NFP .LT. 0 .OR. UNSTAB .AND. ERR .LT. ACCB) GO TO 40
258 IF(.NOT.ZLNEG .OR. UNSTAB .AND. ERR .GT. ACCB) GO TO 20
259 IF(LPR) WRITE(6,1060) '-L',ERR
260 IF(ERR.GT.ACCB) GO TO 280
261 GO TO 40
262C
263C *** evaluate CF1 = f = F'(ZLL,ETA,X)/F(ZLL,ETA,X)
264C
265 20 IF(AXIAL) THEN
266 DX1=X
267 DETA=ETA
268 DZLL=ZLL
269 F=C309R4(DX1,DETA,DZLL,ACC8,SF,RK,ETANE0,LIMIT,ERR,NFP,
270 1 FPMIN,FPMAX,LPR)
271 FCL=SF
272 TPK1=RK
273 ELSE
274 F=C309R5(X,ETA,ZLL,ACC8,FCL,TPK1,ETANE0,LIMIT,ERR,NFP,
275 1 FPMIN,FPMAX,LPR)
276 END IF
277 IF(ERR .GT. ONE) THEN
278 IFAIL=-1
279 GO TO 290
280 END IF
281C
282C *** Make a simple check for CF1 being badly unstable:
283C
284 IF(ID .GE. 0) THEN
285 UNSTAB=DREAL((ONE-ETA*XI)*CI*DIMAG(THETA)/F) .GT. ZERO
286 1 .AND. .NOT.AXIAL .AND. ABS(DIMAG(THETA)) .GT. -LOG(ACC8)*HALF
287 2 .AND. ABSC(ETA)+ABSC(ZLL) .LT. ABSC(X)
288 IF(UNSTAB) THEN
289 ID=-1
290 LF=L1
291 L1=M1
292 RERR=ACCT
293 GO TO 10
294 END IF
295 END IF
296C
297C *** compare accumulated phase FCL with asymptotic phase for G(k+1) :
298C to determine estimate of F(ZLL) (with correct sign) to start recur
299C
300 W=X*X*(HALF/TPK1+ONE/TPK1**2)+ETA*(ETA-TWO*X)/TPK1
301 FESL=(ZLL+ONE)*XLOG+CLL-W-LOG(FCL)
302 40 FESL=FESL-ABS(SCALE)
303 RK=MAX(DREAL(FESL),FPLMIN*HALF)
304 FESL=DCMPLX(MIN(RK,FPLMAX*HALF),DIMAG(FESL))
305 FEST=EXP(FESL)
306C
307 RERR=MAX(RERR,ERR,ACC8*ABS(DREAL(THETA)))
308C
309 FCL=FEST
310 FPL=FCL*F
311 IF(IFCP) FCP(L1)=FPL
312 FC(L1)=FCL
313C
314C *** downward recurrence to lambda = ZLM. array GC,if present,stores RL
315C
316 I=MAX(-ID,0)
317 ZL=ZLL+I
318 MONO=0
319 OFF=ABS(FCL)
320 TA=ABSC(SIGMA)
321
322C
323C CORRESPONDS TO DO 70 L = L1-ID,LF,-ID
324C
325 L=L1-ID
326 LC70=(L1-LF)/ID
327 70 IF(LC70 .LE. 0) GO TO 80
328 IF(ETANE0) THEN
329 IF(RLEL) THEN
330 DSIG=ATAN2(DREAL(ETA),DREAL(ZL))
331 RL=SQRT(DREAL(ZL)**2+DREAL(ETA)**2)
332 ELSE
333 AA=ZL-ETAI
334 BB=ZL+ETAI
335 IF(ABSC(AA) .LT. ACCH .OR. ABSC(BB) .LT. ACCH) GOTO 50
336 DSIG=(LOG(AA)-LOG(BB))*CIH
337 RL=AA*EXP(CI*DSIG)
338 END IF
339 IF(ABSC(SIGMA) .LT. TA*HALF) THEN
340
341C re-calculate SIGMA because of accumulating roundoffs:
342
343 SL=(WLOGAM(ZL+I-ETAI)-WLOGAM(ZL+I+ETAI))*CIH
344 RL=(ZL-ETAI)*EXP(CI*ID*(SIGMA-SL))
345 SIGMA=SL
346 TA=ZERO
347 ELSE
348 SIGMA=SIGMA-DSIG*ID
349 END IF
350 TA=MAX(TA,ABSC(SIGMA))
351 SL=ETA+ZL*ZL*XI
352 PL=ZERO
353 IF(ABSC(ZL) .GT. ACCH) PL=(SL*SL-RL*RL)/ZL
354 FCL1=(FCL*SL+ID*ZL*FPL)/RL
355 SF=ABS(FCL1)
356 IF(SF .GT. FPMAX) GO TO 350
357 FPL=(FPL*SL+ID*PL*FCL)/RL
358 IF(MODE .LE. 1) GCP(L+ID)=PL*ID
359 ELSE
360
361C ETA = 0, including Bessels. NB RL==SL
362
363 RL=ZL*XI
364 FCL1=FCL*RL+FPL*ID
365 SF=ABS(FCL1)
366 IF(SF .GT. FPMAX) GO TO 350
367 FPL=(FCL1*RL-FCL)*ID
368 END IF
369 MONO=MONO+1
370 IF(SF .GE. OFF) MONO=0
371 FCL=FCL1
372 OFF=SF
373 FC(L)=FCL
374 IF(IFCP) FCP(L)=FPL
375 IF(KFN .EQ. 0) SIG(L)=SIGMA
376 IF(MODE .LE. 2) GC(L+ID)=RL
377 ZL=ZL-ID
378 IF(MONO .LT. NDROP .OR. AXIAL .OR.
379 1 DREAL(ZLM)*ID .GT. -NDROP .AND. .NOT.ETANE0) THEN
380 L=L-ID
381 LC70=LC70-1
382 GO TO 70
383 END IF
384 UNSTAB=.TRUE.
385C
386C *** take action if cannot or should not recur below this ZL:
387
388 50 ZLM=ZL
389 LF=L
390 IF(ID .LT. 0) GO TO 380
391 IF(.NOT.UNSTAB) LF=L+1
392 IF(L+MONO .LT. L1-2 .OR. ID .LT. 0 .OR. .NOT.UNSTAB) GO TO 80
393
394C otherwise, all L values (for stability) should be done
395C in the reverse direction:
396
397 ID=-1
398 LF=L1
399 L1=M1
400 RERR=ACCT
401 GO TO 10
402
403 80 IF(FCL .EQ. ZERO) FCL=ACC8
404 F=FPL/FCL
405C
406C *** Check, if second time around, that the 'f' values agree!
407C
408 IF(ID .GT. 0) FIRST=F
409 IF(DONEM) RERR=MAX(RERR,ABSC(F-FIRST)/ABSC(F))
410 IF(DONEM) GO TO 90
411C
412 NOCF2=.FALSE.
413 THETAM=X-ETA*(XLOG+TLOG)-ZLM*HPI+SIGMA
414C
415C *** on left x-plane, determine OMEGA by requiring cut on -x axis
416C on right x-plane, choose OMEGA (using estimate based on THETAM)
417C so H(omega) is smaller and recurs upwards accurately.
418C (x-plane boundary is shifted to give CF2(LH) a chance to converge)
419C
420 OMEGA=SIGN(ONE,DIMAG(X)+ACC8)
421 IF(DREAL(X) .GE. XNEAR) OMEGA=SIGN(ONE,DIMAG(THETAM)+ACC8)
422 SFSH=EXP(OMEGA*SCALE-ABS(SCALE))
423 OFF=EXP(MIN(TWO*MAX(ABS(DIMAG(X)),ABS(DIMAG(THETAM)),
424 1 ABS(DIMAG(ZLM))*3),FPLMAX))
425 EPS=MAX(ACC8,ACCT*HALF/OFF)
426C
427C *** Try first estimated omega, then its opposite,
428C to find the H(omega) linearly independent of F
429C i.e. maximise CF1-CF2 = 1/(F H(omega)) , to minimise H(omega)
430C
431 90 DO 100 L = 1,2
432 LH=1
433 IF(OMEGA .LT. ZERO) LH=2
434 PM=CI*OMEGA
435 ETAP=ETA*PM
436 IF(DONEM) GO TO 130
437 PQ1=ZERO
438 PACCQ=ONE
439 KASE=0
440C
441C *** Check for small X, i.e. whether to avoid CF2 :
442C
443 IF(MODE .GE. 3 .AND. ABSX .LT. ONE ) GO TO 190
444 IF(MODE .LT. 3 .AND. (NOCF2 .OR. ABSX .LT. XNEAR .AND.
445 1 ABSC(ETA)*ABSX .LT. 5 .AND. ABSC(ZLM) .LT. 4)) THEN
446 KASE=5
447 GO TO 120
448 END IF
449C
450C *** Evaluate CF2 : PQ1 = p + i.omega.q at lambda = ZLM
451C
452 PQ1=C309R1(X,ETA,ZLM,PM,EPS,LIMIT,ERR,NPQ(LH),ACC8,ACCH,
453 1 LPR,ACCUR,DELL)
454 ERR=ERR*MAX(ONE,ABSC(PQ1)/MAX(ABSC(F-PQ1),ACC8))
455C
456C *** check if impossible to get F-PQ accurately because of cancellation
457C original guess for OMEGA (based on THETAM) was wrong
458C Use KASE 5 or 6 if necessary if Re(X) < XNEAR
459 IF(ERR .LT. ACCH) GO TO 110
460 NOCF2=DREAL(X) .LT. XNEAR .AND. ABS(DIMAG(X)) .LT. -LOG(ACC8)
461 100 OMEGA=-OMEGA
462 IF(UNSTAB) GO TO 360
463 IF(LPR .AND. DREAL(X) .LT. -XNEAR) WRITE(6,1060) '-X',ERR
464 110 RERR=MAX(RERR,ERR)
465C
466C *** establish case of calculation required for irregular solution
467C
468 120 IF(KASE .GE. 5) GO TO 130
469 IF(DREAL(X) .GT. XNEAR) THEN
470
471C estimate errors if KASE 2 or 3 were to be used:
472
473 PACCQ=EPS*OFF*ABSC(PQ1)/MAX(ABS(DIMAG(PQ1)),ACC8)
474 END IF
475 IF(PACCQ .LT. ACCUR) THEN
476 KASE=2
477 IF(AXIAL) KASE=3
478 ELSE
479 KASE=1
480 IF(NPQ(1)*R20 .LT. JMAX) KASE=4
481C i.e. change to kase=4 if the 2F0 predicted to converge
482 END IF
483 130 GO TO (190,140,150,170,190,190), ABS(KASE)
484 140 IF(.NOT.DONEM)
485C
486C *** Evaluate CF2 : PQ2 = p - i.omega.q at lambda = ZLM (Kase 2)
487C
488 1 PQ2=C309R1(X,ETA,ZLM,-PM,EPS,LIMIT,ERR,NPQ(3-LH),ACC8,ACCH,
489 2 LPR,ACCUR,DELL)
490 P=(PQ2+PQ1)*HALF
491 Q=(PQ2-PQ1)*HALF*PM
492 GO TO 160
493 150 P=DREAL(PQ1)
494 Q=DIMAG(PQ1)
495C
496C *** With Kase = 3 on the real axes, P and Q are real & PQ2 = PQ1*
497C
498 PQ2=DCONJG(PQ1)
499C
500C *** solve for FCM = F at lambda = ZLM,then find norm factor W=FCM/FCL
501C
502 160 W=(PQ1-F)*(PQ2-F)
503 SF=EXP(-ABS(SCALE))
504 FCM=SQRT(Q/W)*SF
505
506C any SQRT given here is corrected by
507C using sign for FCM nearest to phase of FCL
508
509 IF(DREAL(FCM/FCL) .LT. ZERO) FCM=-FCM
510 GAM=(F-P)/Q
511 TA=ABSC(GAM+PM)
512 PACCQ=EPS*MAX(TA,ONE/TA)
513 HCL=FCM*(GAM+PM)*(SFSH/(SF*SF))
514
515C Consider a KASE = 1 Calculation
516
517 IF(PACCQ .GT. ACCUR .AND. KASE .GT. 0) THEN
518 F11V=C309R2(X,ETA,Z11,P11,ACCT,LIMIT,0,ERR,N11,FPMAX,ACC8,ACC16)
519 IF(ERR .LT. PACCQ) GO TO 200
520 END IF
521 RERR=MAX(RERR,PACCQ)
522 GO TO 230
523C
524C *** Arrive here if KASE = 4
525C to evaluate the exponentially decreasing H(LH) directly.
526C
527 170 IF(DONEM) GO TO 180
528 AA=ETAP-ZLM
529 BB=ETAP+ZLM+ONE
530 F20V=C309R3(AA,BB,-HALF*PM*XI,ACCT,JMAX,ERR,FPMAX,N20,XRCF)
531 IF(N20 .LE. 0) GO TO 190
532 RERR=MAX(RERR,ERR)
533 HCL=FPMIN
534 IF(ABS(DREAL(PM*THETAM)+OMEGA*SCALE) .GT. FPLMAX) GO TO 330
535 180 HCL=F20V*EXP(PM*THETAM+OMEGA*SCALE)
536 FCM=SFSH/((F-PQ1)*HCL)
537 GO TO 230
538C
539C *** Arrive here if KASE=1 (or if 2F0 tried mistakenly & failed)
540C
541C for small values of X, calculate F(X,SL) directly from 1F1
542C using DREAL*16 arithmetic if possible.
543C where Z11 = ZLL if ID>0, or = ZLM if ID<0
544C
545 190 F11V=C309R2(X,ETA,Z11,P11,ACCT,LIMIT,0,ERR,N11,FPMAX,ACC8,ACC16)
546 200 IF(N11 .LT. 0) THEN
547
548C F11 failed from BB = negative integer
549
550 IF(LPR) WRITE(6,1060) '-L',ONE
551 IFAIL=-1
552 GO TO 290
553 END IF
554
555C Consider a KASE 2 or 3 calculation :
556
557 IF(ERR .GT. PACCQ .AND. PACCQ .LT. ACCB) THEN
558 KASE=-2
559 IF(AXIAL) KASE=-3
560 GO TO 130
561 END IF
562 RERR=MAX(RERR,ERR)
563 IF(ERR .GT. FPMAX) GO TO 370
564 IF(ID .LT. 0) CLL=Z11*TLOG-HPI*ETA-WLOGAM(BB)
565 1 +WLOGAM(Z11+ONE+P11*ETA)-P11*SIGMA
566 EK=(Z11+ONE)*XLOG-P11*X+CLL-ABS(SCALE)
567 IF(ID .GT. 0) EK=EK-FESL+LOG(FCL)
568 IF(DREAL(EK) .GT. FPLMAX) GO TO 350
569 IF(DREAL(EK) .LT. FPLMIN) GO TO 340
570 FCM=F11V*EXP(EK)
571 IF(KASE .GE. 5) THEN
572 IF(ABSC(ZLM+ZLM-NINTC(ZLM+ZLM)) .LT. ACCH) KASE=6
573C
574C *** For abs(X) < XNEAR, then CF2 may not converge accurately, so
575C *** use an expansion for irregular soln from origin :
576C
577 SL=ZLM
578 ZLNEG=DREAL(ZLM) .LT. -ONE+ACCB
579 IF(KASE .EQ. 5 .OR. ZLNEG) SL=-ZLM-ONE
580 PK=SL+ONE
581 AA=PK-ETAP
582 AB=PK+ETAP
583 BB=TWO*PK
584 CLGAA=WLOGAM(AA)
585 CLGAB=CLGAA
586 IF(ETANE0) CLGAB=WLOGAM(AB)
587 CLGBB=WLOGAM(BB)
588 IF(KASE .EQ. 6 .AND. .NOT.ZLNEG) THEN
589 IF(NPINT(AA,ACCUR)) CLGAA=CLGAB-TWO*PM*SIGMA
590 IF(NPINT(AB,ACCUR)) CLGAB=CLGAA+TWO*PM*SIGMA
591 END IF
592 CLL=SL*TLOG-HPI*ETA-CLGBB+(CLGAA+CLGAB)*HALF
593 DSIG=(CLGAA-CLGAB)*PM*HALF
594 IF(KASE .EQ. 6) P11=-PM
595 EK=PK*XLOG-P11*X+CLL-ABS(SCALE)
596 SF=EXP(-ABS(SCALE))
597 CHI=ZERO
598 IF(.NOT.(KASE .EQ. 5 .OR. ZLNEG)) GO TO 210
599C
600C *** Use G(l) = (cos(CHI) * F(l) - F(-l-1)) / sin(CHI)
601C
602C where CHI = sig(l) - sig(-l-1) - (2l+1)*pi/2
603C
604 CHI=SIGMA-DSIG-(ZLM-SL)*HPI
605 F11V=C309R2(X,ETA,SL,P11,ACCT,LIMIT,0,ERR,NPQ(1),
606 1 FPMAX,ACC8,ACC16)
607 RERR=MAX(RERR,ERR)
608 IF(KASE .EQ. 6) GO TO 210
609 FESL=F11V*EXP(EK)
610 FCL1=EXP(PM*CHI)*FCM
611 HCL=FCL1-FESL
612 RERR=MAX(RERR,ACCT*MAX(ABSC(FCL1),ABSC(FESL))/ABSC(HCL))
613 HCL=HCL/SIN(CHI)*(SFSH/(SF*SF))
614 GO TO 220
615C
616C *** Use the logarithmic expansion for the irregular solution (KASE 6)
617C for the case that BB is integral so sin(CHI) would be zero.
618C
619 210 RL=BB-ONE
620 N=NINTC(RL)
621 ZLOG=XLOG+TLOG-PM*HPI
622 CHI=CHI+PM*THETAM+OMEGA*SCALE+AB*ZLOG
623 AA=ONE-AA
624 IF(NPINT(AA,ACCUR)) THEN
625 HCL=ZERO
626 ELSE
627 IF(ID .GT. 0 .AND. .NOT.ZLNEG) F11V=FCM*EXP(-EK)
628 HCL=EXP(CHI-CLGBB-WLOGAM(AA))*(-1)**(N+1)*(F11V*ZLOG+
629 1 C309R2(X,ETA,SL,-PM,ACCT,LIMIT,2,ERR,NPQ(2),FPMAX,ACC8,ACC16))
630 RERR=MAX(RERR,ERR)
631 END IF
632 IF(N .GT. 0) THEN
633 EK=CHI+WLOGAM(RL)-CLGAB-RL*ZLOG
634 DF=C309R2(X,ETA,-SL-ONE,-PM,ZERO,N,0,ERR,L,FPMAX,ACC8,ACC16)
635 HCL=HCL+EXP(EK)*DF
636 END IF
637 RERR=MAX(RERR,TWO*ABS(BB-NINTC(BB)))
638 220 PQ1=F-SFSH/(FCM*HCL)
639 ELSE
640 IF(MODE .LE. 2) HCL=SFSH/((F-PQ1)*FCM)
641 KASE=1
642 END IF
643C
644C *** Now have absolute normalisations for Coulomb Functions
645C FCM & HCL at lambda = ZLM
646C so determine linear transformations for Functions required :
647C
648 230 IH=ABS(MODE1)/10
649 IF(KFN .EQ. 3) IH=(3-DIMAG(CIK))/2+HALF
650 P11=ONE
651 IF(IH .EQ. 1) P11=CI
652 IF(IH .EQ. 2) P11=-CI
653 DF=-PM
654 IF(IH .GE. 1) DF=-PM+P11
655 IF(ABSC(DF) .LT. ACCH) DF=ZERO
656C
657C *** Normalisations for spherical or cylindrical Bessel functions
658C
659 IF(KFN .LE. 0) THEN
660 ALFA=ZERO
661 BETA=ONE
662 ELSE IF(KFN .EQ. 1) THEN
663 ALFA=XI
664 BETA=XI
665 ELSE
666 ALFA=XI*HALF
667 BETA=SQRT(XI/HPI)
668 IF(DREAL(BETA) .LT. ZERO) BETA=-BETA
669 END IF
670 AA=ONE
671 IF(KFN .GT. 0) AA=-P11*BETA
672
673C Calculate rescaling factors for I & K output
674
675 IF(KFN .GE. 3) THEN
676 P=EXP((ZLM+DELL)*HPI*CIK)
677 AA=BETA*HPI*P
678 BETA=BETA/P
679 Q=CIK*ID
680 END IF
681
682C Calculate rescaling factors for GC output
683
684 IF(IH .EQ. 0) THEN
685 TA=ABS(SCALE)+DIMAG(PM)*SCALE
686 RK=ZERO
687 IF(TA .LT. FPLMAX) RK=EXP(-TA)
688 ELSE
689 TA=ABS(SCALE)+DIMAG(P11)*SCALE
690 IF(ABSC(DF) .GT. ACCH .AND. TA .GT. FPLMAX) GO TO 320
691 IF(ABSC(DF) .GT. ACCH) DF=DF*EXP(TA)
692 SF=TWO*(LH-IH)*SCALE
693 RK=ZERO
694 IF(SF .GT. FPLMAX) GO TO 320
695 IF(SF .GT. FPLMIN) RK=EXP(SF)
696 END IF
697 KAS((3-ID)/2)=KASE
698 W=FCM/FCL
699 IF(LOG(ABSC(W))+LOG(ABSC(FC(LF))) .LT. FPLMIN) GO TO 340
700 IF(MODE .GE. 3) GO TO 240
701 IF(LPR .AND. ABSC(F-PQ1) .LT. ACCH*ABSC(F))
702 1 WRITE(6,1020) LH,ZLM+DELL
703 HPL=HCL*PQ1
704 IF(ABSC(HPL) .LT. FPMIN .OR. ABSC(HCL) .LT. FPMIN) GO TO 330
705C
706C *** IDward recurrence from HCL,HPL(LF) (stored GC(L) is RL if reqd)
707C *** renormalise FC,FCP at each lambda
708C *** ZL = ZLM - MIN(ID,0) here
709C
710 240 DO 270 L = LF,L1,ID
711 FCL=W*FC(L)
712 IF(ABSC(FCL) .LT. FPMIN) GO TO 340
713 IF(IFCP) FPL=W*FCP(L)
714 FC(L)=BETA*FCL
715 IF(IFCP) FCP(L)=BETA*(FPL-ALFA*FCL)*CIK
716 FC(L)=C309R8(FC(L),ACCUR)
717 IF(IFCP) FCP(L)=C309R8(FCP(L),ACCUR)
718 IF(MODE .GE. 3) GO TO 260
719 IF(L .EQ. LF) GO TO 250
720 ZL=ZL+ID
721 ZID=ZL*ID
722 RL=GC(L)
723 IF(ETANE0) THEN
724 SL=ETA+ZL*ZL*XI
725 IF(MODE .EQ. 1) THEN
726 PL=GCP(L)
727 ELSE
728 PL=ZERO
729 IF(ABSC(ZL) .GT. ACCH) PL=(SL*SL-RL*RL)/ZID
730 END IF
731 HCL1=(SL*HCL-ZID*HPL)/RL
732 HPL=(SL*HPL-PL*HCL)/RL
733 ELSE
734 HCL1=RL*HCL-HPL*ID
735 HPL=(HCL-RL*HCL1)*ID
736 END IF
737 HCL=HCL1
738 IF(ABSC(HCL) .GT. FPMAX) GO TO 320
739 250 GC(L)=AA*(RK*HCL+DF*FCL)
740 IF(MODE .EQ. 1) GCP(L)=(AA*(RK*HPL+DF*FPL)-ALFA*GC(L))*CIK
741 GC(L)=C309R8(GC(L),ACCUR)
742 IF(MODE .EQ. 1) GCP(L)=C309R8(GCP(L),ACCUR)
743 IF(KFN .GE. 3) AA=AA*Q
744 260 IF(KFN .GE. 3) BETA=-BETA*Q
745 270 LAST=MIN(LAST,(L1-L)*ID)
746 GO TO 280
747C
748C *** Come here after all soft errors to determine how many L values ok
749C
750 310 IF(LPR) WRITE(6,1000) ZZ
751 GO TO 999
752 320 IF(LPR) WRITE(6,1010) ZL+DELL,'IR',HCL,'>',FPMAX
753 GO TO 280
754 330 IF(LPR) WRITE(6,1010) ZL+DELL,'IR',HCL,'<',FPMIN
755 GO TO 280
756 340 IF(LPR) WRITE(6,1010) ZL+DELL,' ',FCL,'<',FPMIN
757 GO TO 280
758 350 IF(LPR) WRITE(6,1010) ZL+DELL,' ',FCL,'>',FPMAX
759 GO TO 280
760 360 IF(LPR) WRITE(6,1030) ZL+DELL
761 GO TO 280
762 370 IF(LPR) WRITE(6,1040) Z11,I
763 IFAIL=-1
764 GO TO 290
765 380 IF(LPR) WRITE(6,1050) ZLMIN,ZLM,ZLM+ONE,ZLMIN+NL
766 IFAIL=-1
767 GO TO 290
768 280 IF(ID .GT. 0 .OR. LAST .EQ. 0) IFAIL=LAST
769 IF(ID .LT. 0 .AND. LAST .NE. 0) IFAIL=-3
770C
771C *** Come here after ALL errors for this L range (ZLM,ZLL)
772C
773C *** so on first block, 'F' started decreasing monotonically,
774C or hit bound states for low ZL.
775C thus redo M1 to LF-1 in reverse direction, i.e. do
776C CF1A at ZLMIN & CF2 at ZLM (midway between ZLMIN & ZLMAX)
777C
778 290 IF(ID .GT. 0 .AND. LF .NE. M1) THEN
779 ID=-1
780 IF(.NOT.UNSTAB) LF=LF-1
781 DONEM=UNSTAB
782 LF=MIN(LF,L1)
783 L1=M1
784 GO TO 10
785 END IF
786 IF(IFAIL .LT. 0) GO TO 999
787 IF(LPR .AND. RERR .GT. ACCB) WRITE(6,1070) RERR
788 IF(RERR .GT. 0.1D0) IFAIL=-4
789 999 DO 998 L = 1,NL+1
790 FC(L-1)=FC(L)
791 GC(L-1)=GC(L)
792 FCP(L-1)=FCP(L)
793 GCP(L-1)=GCP(L)
794 998 SIG(L-1)=SIG(L)
795 RETURN
796C
797 1000 FORMAT(1X,'***** CERN C309 WCLBES ... ',
798 1 'CANNOT CALCULATE IRREGULAR SOLUTIONS FOR X =',
799 2 1P,2D10.2,' ABS(X) TOO SMALL')
800 1010 FORMAT(1X,'***** CERN C309 WCLBES ... ',
801 1 'AT ZL =',2F8.3,' ',A2,'REGULAR SOLUTION (',1P,2E10.1,') ',
802 2 A1,E10.1)
803 1020 FORMAT(1X,'***** CERN C309 WCLBES ... ',
804 1 'WARNING: LINEAR INDEPENDENCE BETWEEN ''F'' AND ''H(',I1,
805 2 ')'' IS LOST AT ZL = ',2F7.2/1X,'*****',22X,'(EG. COULOMB ',
806 3 'EIGENSTATE OR CF1 UNSTABLE)')
807 1030 FORMAT(1X,'***** CERN C309 WCLBES ... ',
808 2 '(ETA & L)/X TOO LARGE FOR CF1A, AND CF1 UNSTABLE AT L = ',2F8.2)
809 1040 FORMAT(1X,'***** CERN C309 WCLBES ... ',
810 1 'OVERFLOW IN 1F1 SERIES AT ZL = ',2F8.3,' AT TERM',I5)
811 1050 FORMAT(1X,'***** CERN C309 WCLBES ... ',
812 1 'BOTH BOUND-STATE POLES AND F-INSTABILITIES OCCUR OR MULTIPLE',
813 2 ' INSTABILITIES PRESENT'/
814 3 1X,'*****',22X,'TRY CALLING TWICE, 1ST FOR ZL FROM',2F8.3,
815 4 ' TO',2F8.3,' (INCL)'/1X,'*****',41X,'2ND FOR ZL FROM',2F8.3,
816 5 ' TO',2F8.3)
817 1060 FORMAT(1X,'***** CERN C309 WCLBES ... ',
818 1 'WARNING: AS ''',A2,''' REFLECTION RULES NOT USED ',
819 2 'ERRORS CAN BE UP TO',1PD12.2)
820 1070 FORMAT(1X,'***** CERN C309 WCLBES ... ',
821 1 'WARNING: OVERALL ROUNDOFF ERROR APPROXIMATELY',1PE11.1)
822 END
823#endif