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