github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zbuni.f (about) 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,ZABS,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, ZABS, 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 = ZABS(CMPLX(CYR(1),CYI(1),kind=KIND(1.0D0))) 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/ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 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