MINICERN/mathlib/gen/c/wclbes.F

   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)
  13 C
  14 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  15 C                                                                      C
  16 C  COMPLEX COULOMB WAVEFUNCTION PROGRAM USING STEED'S METHOD           C
  17 C  Original title : COULCC                                             C
  18 C                                                                      C
  19 C  A. R. Barnett           Manchester  March   1981                    C
  20 C  modified I.J. Thompson  Daresbury, Sept. 1983 for Complex Functions C
  21 C                                                                      C
  22 C  The FORM (not the SUBSTANCE) of this program has been modified      C
  23 C   by K.S. KOLBIG (CERN)    December 1987                             C
  24 C                                                                      C
  25 C  original program  RCWFN       in    CPC  8 (1974) 377-395           C
  26 C                 +  RCWFF       in    CPC 11 (1976) 141-142           C
  27 C                 +  COULFG      in    CPC 27 (1982) 147-166           C
  28 C  description of real algorithm in    CPC 21 (1981) 297-314           C
  29 C  description of complex algorithm    JCP 64 (1986) 490-509           C
  30 C  this version written up       in    CPC 36 (1985) 363-372           C
  31 C                                                                      C
  32 C  WCLBES returns F,G,G',G',SIG for complex ETA, ZZ, and ZLMIN,        C
  33 C   for NL integer-spaced lambda values ZLMIN to ZLMIN+NL inclusive,   C
  34 C   thus giving  complex-energy solutions to the Coulomb Schrodinger   C
  35 C   equation,to the Klein-Gordon equation and to suitable forms of     C
  36 C   the Dirac equation ,also spherical & cylindrical Bessel equations  C
  37 C                                                                      C
  38 C  if ABS(MODE1)                                                       C
  39 C            = 1  get F,G,F',G'   for integer-spaced lambda values     C
  40 C            = 2      F,G      unused arrays must be dimensioned in    C
  41 C            = 3      F,  F'          call to at least length (0:1)    C
  42 C            = 4      F                                                C
  43 C            = 11 get F,H+,F',H+' ) if KFN=0, H+ = G + i.F        )    C
  44 C            = 12     F,H+        )       >0, H+ = J + i.Y = H(1) ) in C
  45 C            = 21 get F,H-,F',H-' ) if KFN=0, H- = G - i.F        ) GC C
  46 C            = 22     F,H-        )       >0, H- = J - i.Y = H(2) )    C
  47 C                                                                      C
  48 C     if MODE1 < 0 then the values returned are scaled by an exponen-  C
  49 C                  tial factor (dependent only on ZZ) to bring nearer  C
  50 C                  unity the functions for large ABS(ZZ), small ETA ,  C
  51 C                  and ABS(ZL) < ABS(ZZ).                              C
  52 C        Define SCALE = (  0        if MODE1 > 0                       C
  53 C                       (  IMAG(ZZ) if MODE1 < 0  &  KFN < 3           C
  54 C                       (  REAL(ZZ) if MODE1 < 0  &  KFN = 3           C
  55 C        then FC = EXP(-ABS(SCALE)) * ( F, j, J, or I)                 C
  56 C         and GC = EXP(-ABS(SCALE)) * ( G, y, or Y )                   C
  57 C               or EXP(SCALE)       * ( H+, H(1), or K)                C
  58 C               or EXP(-SCALE)      * ( H- or H(2) )                   C
  59 C                                                                      C
  60 C  if  KFN  =  0,-1  complex Coulomb functions are returned   F & G    C
  61 C           =  1   spherical Bessel      "      "     "       j & y    C
  62 C           =  2 cylindrical Bessel      "      "     "       J & Y    C
  63 C           =  3 modified cyl. Bessel    "      "     "       I & K    C
  64 C                                                                      C
  65 C          and where Coulomb phase shifts put in SIG if KFN=0 (not -1) C
  66 C                                                                      C
  67 C  The use of MODE and KFN is independent                              C
  68 C    (except that for KFN=3,  H(1) & H(2) are not given)               C
  69 C                                                                      C
  70 C  With negative orders lambda, WCLBES can still be used but with      C
  71 C    reduced accuracy as CF1 becomes unstable. The user is thus        C
  72 C    strongly advised to use reflection formulae based on              C
  73 C    H+-(ZL,,) = H+-(-ZL-1,,) * exp +-i(sig(ZL)-sig(-ZL-1)-(ZL+1/2)pi) C
  74 C                                                                      C
  75 C  Precision:  results to within 2-3 decimals of 'machine accuracy',   C
  76 C              except in the following cases:                          C
  77 C              (1) if CF1A fails because X too small or ETA too large  C
  78 C               the F solution  is less accurate if it decreases with  C
  79 C               decreasing lambda (e.g. for lambda.LE.-1 & ETA.NE.0)   C
  80 C              (2) if ETA is large (e.g. >> 50) and X inside the       C
  81 C                turning point, then progressively less accuracy       C
  82 C              (3) if ZLMIN is around sqrt(ACCUR) distance from an     C
  83 C               integral order and abs(X) << 1, then errors present.   C
  84 C               RERR traces the main roundoff errors.                  C
  85 C                                                                      C
  86 C   WCLBES is coded for real*8 on IBM or equivalent  ACCUR >= 10**-14  C
  87 C          with a section of doubled REAL*16 for less roundoff errors. C
  88 C          (If no doubled precision available, increase JMAX to eg 100)C
  89 C   Use IMPLICIT COMPLEX*32 & REAL*16 on VS compiler ACCUR >= 10**-32  C
  90 C   For single precision CDC (48 bits) reassign                        C
  91 C        DOUBLE PRECISION=REAL etc.                                    C
  92 C                                                                      C
  93 C   IPR    on input   = 0 : no printing of error messages              C
  94 C                    ne 0 : print error messages on file 6             C
  95 C   IFAIL  in output = -2 : argument out of range                      C
  96 C                    = -1 : one of the continued fractions failed,     C
  97 C                           or arithmetic check before final recursion C
  98 C                    =  0 : All Calculations satisfactory              C
  99 C                    ge 0 : results available for orders up to & at    C
 100 C                             position NL-IFAIL in the output arrays.  C
 101 C                    = -3 : values at ZLMIN not found as over/underflowC
 102 C                    = -4 : roundoff errors make results meaningless   C
 103 C                                                                      C
 104 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 105 C                                                                      C
 106 C     Machine dependent constants :                                    C
 107 C                                                                      C
 108 C     ACCU     target bound on relative error (except near 0 crossings)C
 109 C               (ACCUR should be at least 100 * ACC8)                  C
 110 C     ACC8    smallest number with 1+ACC8 .ne.1 in REAL*8  arithmetic  C
 111 C     ACC16    smallest number with 1+ACC16.ne.1 in REAL*16 arithmetic C
 112 C     FPMAX    magnitude of largest floating point number * ACC8       C
 113 C     FPMIN    magnitude of smallest floating point number / ACC8      C
 114 C                                                                      C
 115 C     Parameters determining region of calculations :                  C
 116 C                                                                      C
 117 C        R20      estimate of (2F0 iterations)/(CF2 iterations)        C
 118 C        ASYM     minimum X/(ETA**2+L) for CF1A to converge easily     C
 119 C        XNEAR    minimum ABS(X) for CF2 to converge accurately        C
 120 C        LIMIT    maximum no. iterations for CF1, CF2, and 1F1 series  C
 121 C        JMAX     size of work arrays for Pade accelerations           C
 122 C        NDROP    number of successive decrements to define instabilityC
 123 C                                                                      C
 124 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 125 C
 126 C     C309R1 = CF2,   C309R2 = F11,    C309R3 = F20,
 127 C     C309R4 = CF1R,  C309R5 = CF1C,   C309R6 = CF1A,
 128 C     C309R7 = RCF,   C309R8 = CTIDY.
 129 C
 130       IMPLICIT COMPLEX*16 (A-H,O-Z)
 131 C
 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
 166 C
 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)
 186 C
 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)
 199 C
 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
 209 C       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)
 218 C
 219 C ***       ZLL  is final lambda value, or 0.5 smaller for J,Y Bessels
 220 C
 221       Z11=ZLL
 222       IF(ID .LT. 0) Z11=ZLM
 223       P11=CI*SIGN(ONE,ACC8-DIMAG(ETA))
 224       LAST=L1
 225 C
 226 C ***       Find phase shifts and Gamow factor at lambda = ZLL
 227 C
 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
 241 C
 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
 244 C
 245 C ***         use CF1 instead of CF1A, if predicted to converge faster,
 246 C                 (otherwise using CF1A as it treats negative lambda &
 247 C                  recurrence-unstable cases properly)
 248 C
 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
 262 C
 263 C ***    evaluate CF1  =  f   =  F'(ZLL,ETA,X)/F(ZLL,ETA,X)
 264 C
 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
 281 C
 282 C ***  Make a simple check for CF1 being badly unstable:
 283 C
 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
 296 C
 297 C *** compare accumulated phase FCL with asymptotic phase for G(k+1) :
 298 C     to determine estimate of F(ZLL) (with correct sign) to start recur
 299 C
 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)
 306 C
 307       RERR=MAX(RERR,ERR,ACC8*ABS(DREAL(THETA)))
 308 C
 309       FCL=FEST
 310       FPL=FCL*F
 311       IF(IFCP) FCP(L1)=FPL
 312       FC(L1)=FCL
 313 C
 314 C *** downward recurrence to lambda = ZLM. array GC,if present,stores RL
 315 C
 316       I=MAX(-ID,0)
 317       ZL=ZLL+I
 318       MONO=0
 319       OFF=ABS(FCL)
 320       TA=ABSC(SIGMA)
 321
 322 C
 323 C     CORRESPONDS TO   DO 70 L = L1-ID,LF,-ID
 324 C
 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
 341 C               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
 361 C                      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.
 385 C
 386 C ***    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
 394 C             otherwise, all L values (for stability) should be done
 395 C                        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
 405 C
 406 C *** Check, if second time around, that the 'f' values agree!
 407 C
 408       IF(ID .GT. 0) FIRST=F
 409       IF(DONEM) RERR=MAX(RERR,ABSC(F-FIRST)/ABSC(F))
 410       IF(DONEM) GO TO 90
 411 C
 412       NOCF2=.FALSE.
 413       THETAM=X-ETA*(XLOG+TLOG)-ZLM*HPI+SIGMA
 414 C
 415 C *** on left x-plane, determine OMEGA by requiring cut on -x axis
 416 C     on right x-plane, choose OMEGA (using estimate based on THETAM)
 417 C       so H(omega) is smaller and recurs upwards accurately.
 418 C     (x-plane boundary is shifted to give CF2(LH) a chance to converge)
 419 C
 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)
 426 C
 427 C ***    Try first estimated omega, then its opposite,
 428 C        to find the H(omega) linearly independent of F
 429 C        i.e. maximise  CF1-CF2 = 1/(F H(omega)) , to minimise H(omega)
 430 C
 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
 440 C
 441 C ***            Check for small X, i.e. whether to avoid CF2 :
 442 C
 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
 449 C
 450 C ***  Evaluate   CF2 : PQ1 = p + i.omega.q  at lambda = ZLM
 451 C
 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))
 455 C
 456 C *** check if impossible to get F-PQ accurately because of cancellation
 457 C                original guess for OMEGA (based on THETAM) was wrong
 458 C                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)
 465 C
 466 C ***  establish case of calculation required for irregular solution
 467 C
 468   120 IF(KASE .GE. 5) GO TO 130
 469       IF(DREAL(X) .GT. XNEAR) THEN
 470
 471 C          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
 481 C             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)
 485 C
 486 C ***  Evaluate   CF2 : PQ2 = p - i.omega.q  at lambda = ZLM   (Kase 2)
 487 C
 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)
 495 C
 496 C ***   With Kase = 3 on the real axes, P and Q are real & PQ2 = PQ1*
 497 C
 498       PQ2=DCONJG(PQ1)
 499 C
 500 C *** solve for FCM = F at lambda = ZLM,then find norm factor W=FCM/FCL
 501 C
 502   160 W=(PQ1-F)*(PQ2-F)
 503       SF=EXP(-ABS(SCALE))
 504       FCM=SQRT(Q/W)*SF
 505
 506 C                  any SQRT given here is corrected by
 507 C                  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
 515 C                            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
 523 C
 524 C *** Arrive here if KASE = 4
 525 C     to evaluate the exponentially decreasing H(LH) directly.
 526 C
 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
 538 C
 539 C *** Arrive here if KASE=1   (or if 2F0 tried mistakenly & failed)
 540 C
 541 C           for small values of X, calculate F(X,SL) directly from 1F1
 542 C               using DREAL*16 arithmetic if possible.
 543 C           where Z11 = ZLL if ID>0, or = ZLM if ID<0
 544 C
 545   190 F11V=C309R2(X,ETA,Z11,P11,ACCT,LIMIT,0,ERR,N11,FPMAX,ACC8,ACC16)
 546   200 IF(N11 .LT. 0) THEN
 547
 548 C                               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
 555 C                      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
 573 C
 574 C ***  For abs(X) < XNEAR, then CF2 may not converge accurately, so
 575 C ***      use an expansion for irregular soln from origin :
 576 C
 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
 599 C
 600 C ***  Use  G(l)  =  (cos(CHI) * F(l) - F(-l-1)) /  sin(CHI)
 601 C
 602 C      where CHI = sig(l) - sig(-l-1) - (2l+1)*pi/2
 603 C
 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
 615 C
 616 C *** Use the logarithmic expansion for the irregular solution (KASE 6)
 617 C        for the case that BB is integral so sin(CHI) would be zero.
 618 C
 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
 643 C
 644 C ***  Now have absolute normalisations for Coulomb Functions
 645 C          FCM & HCL  at lambda = ZLM
 646 C      so determine linear transformations for Functions required :
 647 C
 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
 656 C
 657 C *** Normalisations for spherical or cylindrical Bessel functions
 658 C
 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
 673 C                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
 682 C                  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
 705 C
 706 C *** IDward recurrence from HCL,HPL(LF) (stored GC(L) is RL if reqd)
 707 C *** renormalise FC,FCP at each lambda
 708 C ***    ZL   = ZLM - MIN(ID,0) here
 709 C
 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
 747 C
 748 C *** Come here after all soft errors to determine how many L values ok
 749 C
 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
 770 C
 771 C *** Come here after ALL errors for this L range (ZLM,ZLL)
 772 C
 773 C *** so on first block, 'F' started decreasing monotonically,
 774 C                        or hit bound states for low ZL.
 775 C     thus redo M1 to LF-1 in reverse direction, i.e. do
 776 C      CF1A at ZLMIN & CF2 at ZLM (midway between ZLMIN & ZLMAX)
 777 C
 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
 796 C
 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