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