src/numerical/slatec/fortran/zbuni.f

   1 *DECK ZBUNI
   2       SUBROUTINE ZBUNI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST,
   3      +   FNUL, TOL, ELIM, ALIM)
   4 C***BEGIN PROLOGUE  ZBUNI
   5 C***SUBSIDIARY
   6 C***PURPOSE  Subsidiary to ZBESI and ZBESK
   7 C***LIBRARY   SLATEC
   8 C***TYPE      ALL (CBUNI-A, ZBUNI-A)
   9 C***AUTHOR  Amos, D. E., (SNL)
  10 C***DESCRIPTION
  11 C
  12 C     ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE ABS(Z).GT.
  13 C     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
  14 C     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
  15 C     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
  16 C     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
  17 C
  18 C***SEE ALSO  ZBESI, ZBESK
  19 C***ROUTINES CALLED  D1MACH, ZABS, ZUNI1, ZUNI2
  20 C***REVISION HISTORY  (YYMMDD)
  21 C   830501  DATE WRITTEN
  22 C   910415  Prologue converted to Version 4.0 format.  (BAB)
  23 C***END PROLOGUE  ZBUNI
  24 C     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z
  25       DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU,
  26      * ELIM, FNU, FNUI, FNUL, GNU, RAZ, RZI, RZR, STI, STR, S1I, S1R,
  27      * S2I, S2R, TOL, YI, YR, ZI, ZR, ZABS, ASCLE, BRY, C1R, C1I, C1M,
  28      * D1MACH
  29       INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
  30       DIMENSION YR(N), YI(N), CYR(2), CYI(2), BRY(3)
  31       EXTERNAL ZABS
  32 C***FIRST EXECUTABLE STATEMENT  ZBUNI
  33       NZ = 0
  34       AX = ABS(ZR)*1.7321D0
  35       AY = ABS(ZI)
  36       IFORM = 1
  37       IF (AY.GT.AX) IFORM = 2
  38       IF (NUI.EQ.0) GO TO 60
  39       FNUI = NUI
  40       DFNU = FNU + (N-1)
  41       GNU = DFNU + FNUI
  42       IF (IFORM.EQ.2) GO TO 10
  43 C-----------------------------------------------------------------------
  44 C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
  45 C     -PI/3.LE.ARG(Z).LE.PI/3
  46 C-----------------------------------------------------------------------
  47       CALL ZUNI1(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
  48      * ELIM, ALIM)
  49       GO TO 20
  50    10 CONTINUE
  51 C-----------------------------------------------------------------------
  52 C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
  53 C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
  54 C     AND HPI=PI/2
  55 C-----------------------------------------------------------------------
  56       CALL ZUNI2(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
  57      * ELIM, ALIM)
  58    20 CONTINUE
  59       IF (NW.LT.0) GO TO 50
  60       IF (NW.NE.0) GO TO 90
  61       STR = ZABS(CYR(1),CYI(1))
  62 C----------------------------------------------------------------------
  63 C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
  64 C----------------------------------------------------------------------
  65       BRY(1)=1.0D+3*D1MACH(1)/TOL
  66       BRY(2) = 1.0D0/BRY(1)
  67       BRY(3) = BRY(2)
  68       IFLAG = 2
  69       ASCLE = BRY(2)
  70       CSCLR = 1.0D0
  71       IF (STR.GT.BRY(1)) GO TO 21
  72       IFLAG = 1
  73       ASCLE = BRY(1)
  74       CSCLR = 1.0D0/TOL
  75       GO TO 25
  76    21 CONTINUE
  77       IF (STR.LT.BRY(2)) GO TO 25
  78       IFLAG = 3
  79       ASCLE=BRY(3)
  80       CSCLR = TOL
  81    25 CONTINUE
  82       CSCRR = 1.0D0/CSCLR
  83       S1R = CYR(2)*CSCLR
  84       S1I = CYI(2)*CSCLR
  85       S2R = CYR(1)*CSCLR
  86       S2I = CYI(1)*CSCLR
  87       RAZ = 1.0D0/ZABS(ZR,ZI)
  88       STR = ZR*RAZ
  89       STI = -ZI*RAZ
  90       RZR = (STR+STR)*RAZ
  91       RZI = (STI+STI)*RAZ
  92       DO 30 I=1,NUI
  93         STR = S2R
  94         STI = S2I
  95         S2R = (DFNU+FNUI)*(RZR*STR-RZI*STI) + S1R
  96         S2I = (DFNU+FNUI)*(RZR*STI+RZI*STR) + S1I
  97         S1R = STR
  98         S1I = STI
  99         FNUI = FNUI - 1.0D0
 100         IF (IFLAG.GE.3) GO TO 30
 101         STR = S2R*CSCRR
 102         STI = S2I*CSCRR
 103         C1R = ABS(STR)
 104         C1I = ABS(STI)
 105         C1M = MAX(C1R,C1I)
 106         IF (C1M.LE.ASCLE) GO TO 30
 107         IFLAG = IFLAG+1
 108         ASCLE = BRY(IFLAG)
 109         S1R = S1R*CSCRR
 110         S1I = S1I*CSCRR
 111         S2R = STR
 112         S2I = STI
 113         CSCLR = CSCLR*TOL
 114         CSCRR = 1.0D0/CSCLR
 115         S1R = S1R*CSCLR
 116         S1I = S1I*CSCLR
 117         S2R = S2R*CSCLR
 118         S2I = S2I*CSCLR
 119    30 CONTINUE
 120       YR(N) = S2R*CSCRR
 121       YI(N) = S2I*CSCRR
 122       IF (N.EQ.1) RETURN
 123       NL = N - 1
 124       FNUI = NL
 125       K = NL
 126       DO 40 I=1,NL
 127         STR = S2R
 128         STI = S2I
 129         S2R = (FNU+FNUI)*(RZR*STR-RZI*STI) + S1R
 130         S2I = (FNU+FNUI)*(RZR*STI+RZI*STR) + S1I
 131         S1R = STR
 132         S1I = STI
 133         STR = S2R*CSCRR
 134         STI = S2I*CSCRR
 135         YR(K) = STR
 136         YI(K) = STI
 137         FNUI = FNUI - 1.0D0
 138         K = K - 1
 139         IF (IFLAG.GE.3) GO TO 40
 140         C1R = ABS(STR)
 141         C1I = ABS(STI)
 142         C1M = MAX(C1R,C1I)
 143         IF (C1M.LE.ASCLE) GO TO 40
 144         IFLAG = IFLAG+1
 145         ASCLE = BRY(IFLAG)
 146         S1R = S1R*CSCRR
 147         S1I = S1I*CSCRR
 148         S2R = STR
 149         S2I = STI
 150         CSCLR = CSCLR*TOL
 151         CSCRR = 1.0D0/CSCLR
 152         S1R = S1R*CSCLR
 153         S1I = S1I*CSCLR
 154         S2R = S2R*CSCLR
 155         S2I = S2I*CSCLR
 156    40 CONTINUE
 157       RETURN
 158    50 CONTINUE
 159       NZ = -1
 160       IF(NW.EQ.(-2)) NZ=-2
 161       RETURN
 162    60 CONTINUE
 163       IF (IFORM.EQ.2) GO TO 70
 164 C-----------------------------------------------------------------------
 165 C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
 166 C     -PI/3.LE.ARG(Z).LE.PI/3
 167 C-----------------------------------------------------------------------
 168       CALL ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
 169      * ELIM, ALIM)
 170       GO TO 80
 171    70 CONTINUE
 172 C-----------------------------------------------------------------------
 173 C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
 174 C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
 175 C     AND HPI=PI/2
 176 C-----------------------------------------------------------------------
 177       CALL ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
 178      * ELIM, ALIM)
 179    80 CONTINUE
 180       IF (NW.LT.0) GO TO 50
 181       NZ = NW
 182       RETURN
 183    90 CONTINUE
 184       NLAST = N
 185       RETURN
 186       END