amos/zbuni.f

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