1*DECK ZSERI 2 SUBROUTINE ZSERI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, 3 + ALIM) 4C***BEGIN PROLOGUE ZSERI 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to ZBESI and ZBESK 7C***LIBRARY SLATEC 8C***TYPE ALL (CSERI-A, ZSERI-A) 9C***AUTHOR Amos, D. E., (SNL) 10C***DESCRIPTION 11C 12C ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY 13C MEANS OF THE POWER SERIES FOR LARGE ABS(Z) IN THE 14C REGION ABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN. 15C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO 16C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE 17C CONDITION ABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE 18C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). 19C 20C***SEE ALSO ZBESI, ZBESK 21C***ROUTINES CALLED D1MACH, DGAMLN, ZABS, ZDIV, ZLOG, ZMLT, ZUCHK 22C***REVISION HISTORY (YYMMDD) 23C 830501 DATE WRITTEN 24C 910415 Prologue converted to Version 4.0 format. (BAB) 25C 930122 Added ZLOG to EXTERNAL statement. (RWC) 26C***END PROLOGUE ZSERI 27C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z 28 DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL, 29 * AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, 30 * ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, 31 * STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI, 32 * ZR, DGAMLN, D1MACH, ZABS 33 INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW 34 DIMENSION YR(N), YI(N), WR(2), WI(2) 35 EXTERNAL ZABS, ZLOG 36 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / 37C***FIRST EXECUTABLE STATEMENT ZSERI 38 NZ = 0 39 AZ = ZABS(ZR,ZI) 40 IF (AZ.EQ.0.0D0) GO TO 160 41 ARM = 1.0D+3*D1MACH(1) 42 RTR1 = SQRT(ARM) 43 CRSCR = 1.0D0 44 IFLAG = 0 45 IF (AZ.LT.ARM) GO TO 150 46 HZR = 0.5D0*ZR 47 HZI = 0.5D0*ZI 48 CZR = ZEROR 49 CZI = ZEROI 50 IF (AZ.LE.RTR1) GO TO 10 51 CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI) 52 10 CONTINUE 53 ACZ = ZABS(CZR,CZI) 54 NN = N 55 CALL ZLOG(HZR, HZI, CKR, CKI, IDUM) 56 20 CONTINUE 57 DFNU = FNU + (NN-1) 58 FNUP = DFNU + 1.0D0 59C----------------------------------------------------------------------- 60C UNDERFLOW TEST 61C----------------------------------------------------------------------- 62 AK1R = CKR*DFNU 63 AK1I = CKI*DFNU 64 AK = DGAMLN(FNUP,IDUM) 65 AK1R = AK1R - AK 66 IF (KODE.EQ.2) AK1R = AK1R - ZR 67 IF (AK1R.GT.(-ELIM)) GO TO 40 68 30 CONTINUE 69 NZ = NZ + 1 70 YR(NN) = ZEROR 71 YI(NN) = ZEROI 72 IF (ACZ.GT.DFNU) GO TO 190 73 NN = NN - 1 74 IF (NN.EQ.0) RETURN 75 GO TO 20 76 40 CONTINUE 77 IF (AK1R.GT.(-ALIM)) GO TO 50 78 IFLAG = 1 79 SS = 1.0D0/TOL 80 CRSCR = TOL 81 ASCLE = ARM*SS 82 50 CONTINUE 83 AA = EXP(AK1R) 84 IF (IFLAG.EQ.1) AA = AA*SS 85 COEFR = AA*COS(AK1I) 86 COEFI = AA*SIN(AK1I) 87 ATOL = TOL*ACZ/FNUP 88 IL = MIN(2,NN) 89 DO 90 I=1,IL 90 DFNU = FNU + (NN-I) 91 FNUP = DFNU + 1.0D0 92 S1R = CONER 93 S1I = CONEI 94 IF (ACZ.LT.TOL*FNUP) GO TO 70 95 AK1R = CONER 96 AK1I = CONEI 97 AK = FNUP + 2.0D0 98 S = FNUP 99 AA = 2.0D0 100 60 CONTINUE 101 RS = 1.0D0/S 102 STR = AK1R*CZR - AK1I*CZI 103 STI = AK1R*CZI + AK1I*CZR 104 AK1R = STR*RS 105 AK1I = STI*RS 106 S1R = S1R + AK1R 107 S1I = S1I + AK1I 108 S = S + AK 109 AK = AK + 2.0D0 110 AA = AA*ACZ*RS 111 IF (AA.GT.ATOL) GO TO 60 112 70 CONTINUE 113 S2R = S1R*COEFR - S1I*COEFI 114 S2I = S1R*COEFI + S1I*COEFR 115 WR(I) = S2R 116 WI(I) = S2I 117 IF (IFLAG.EQ.0) GO TO 80 118 CALL ZUCHK(S2R, S2I, NW, ASCLE, TOL) 119 IF (NW.NE.0) GO TO 30 120 80 CONTINUE 121 M = NN - I + 1 122 YR(M) = S2R*CRSCR 123 YI(M) = S2I*CRSCR 124 IF (I.EQ.IL) GO TO 90 125 CALL ZDIV(COEFR, COEFI, HZR, HZI, STR, STI) 126 COEFR = STR*DFNU 127 COEFI = STI*DFNU 128 90 CONTINUE 129 IF (NN.LE.2) RETURN 130 K = NN - 2 131 AK = K 132 RAZ = 1.0D0/AZ 133 STR = ZR*RAZ 134 STI = -ZI*RAZ 135 RZR = (STR+STR)*RAZ 136 RZI = (STI+STI)*RAZ 137 IF (IFLAG.EQ.1) GO TO 120 138 IB = 3 139 100 CONTINUE 140 DO 110 I=IB,NN 141 YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) 142 YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) 143 AK = AK - 1.0D0 144 K = K - 1 145 110 CONTINUE 146 RETURN 147C----------------------------------------------------------------------- 148C RECUR BACKWARD WITH SCALED VALUES 149C----------------------------------------------------------------------- 150 120 CONTINUE 151C----------------------------------------------------------------------- 152C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE 153C UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3 154C----------------------------------------------------------------------- 155 S1R = WR(1) 156 S1I = WI(1) 157 S2R = WR(2) 158 S2I = WI(2) 159 DO 130 L=3,NN 160 CKR = S2R 161 CKI = S2I 162 S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI) 163 S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR) 164 S1R = CKR 165 S1I = CKI 166 CKR = S2R*CRSCR 167 CKI = S2I*CRSCR 168 YR(K) = CKR 169 YI(K) = CKI 170 AK = AK - 1.0D0 171 K = K - 1 172 IF (ZABS(CKR,CKI).GT.ASCLE) GO TO 140 173 130 CONTINUE 174 RETURN 175 140 CONTINUE 176 IB = L + 1 177 IF (IB.GT.NN) RETURN 178 GO TO 100 179 150 CONTINUE 180 NZ = N 181 IF (FNU.EQ.0.0D0) NZ = NZ - 1 182 160 CONTINUE 183 YR(1) = ZEROR 184 YI(1) = ZEROI 185 IF (FNU.NE.0.0D0) GO TO 170 186 YR(1) = CONER 187 YI(1) = CONEI 188 170 CONTINUE 189 IF (N.EQ.1) RETURN 190 DO 180 I=2,N 191 YR(I) = ZEROR 192 YI(I) = ZEROI 193 180 CONTINUE 194 RETURN 195C----------------------------------------------------------------------- 196C RETURN WITH NZ.LT.0 IF ABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE 197C THE CALCULATION IN CBINU WITH N=N-ABS(NZ) 198C----------------------------------------------------------------------- 199 190 CONTINUE 200 NZ = -NZ 201 RETURN 202 END 203