1*DECK ZUNI2 2 SUBROUTINE ZUNI2 (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, 3 + TOL, ELIM, ALIM) 4C***BEGIN PROLOGUE ZUNI2 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to ZBESI and ZBESK 7C***LIBRARY SLATEC 8C***TYPE ALL (CUNI2-A, ZUNI2-A) 9C***AUTHOR Amos, D. E., (SNL) 10C***DESCRIPTION 11C 12C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF 13C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I 14C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. 15C 16C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC 17C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. 18C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER 19C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. 20C Y(I)=CZERO FOR I=NLAST+1,N 21C 22C***SEE ALSO ZBESI, ZBESK 23C***ROUTINES CALLED D1MACH, ZABS, ZAIRY, ZUCHK, ZUNHJ, ZUOIK 24C***REVISION HISTORY (YYMMDD) 25C 830501 DATE WRITTEN 26C 910415 Prologue converted to Version 4.0 format. (BAB) 27C***END PROLOGUE ZUNI2 28C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, 29C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN 30 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI, 31 * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR, 32 * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII, 33 * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, 34 * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, 35 * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, 36 * CYI, D1MACH, ZABS, CAR, SAR 37 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, 38 * NN, NUF, NW, NZ, IDUM 39 DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), 40 * CSRR(3), CYR(2), CYI(2) 41 EXTERNAL ZABS 42 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / 43 DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), 44 * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/ 45 DATA HPI, AIC / 46 1 1.57079632679489662D+00, 1.265512123484645396D+00/ 47C***FIRST EXECUTABLE STATEMENT ZUNI2 48 NZ = 0 49 ND = N 50 NLAST = 0 51C----------------------------------------------------------------------- 52C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- 53C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, 54C EXP(ALIM)=EXP(ELIM)*TOL 55C----------------------------------------------------------------------- 56 CSCL = 1.0D0/TOL 57 CRSC = TOL 58 CSSR(1) = CSCL 59 CSSR(2) = CONER 60 CSSR(3) = CRSC 61 CSRR(1) = CRSC 62 CSRR(2) = CONER 63 CSRR(3) = CSCL 64 BRY(1) = 1.0D+3*D1MACH(1)/TOL 65C----------------------------------------------------------------------- 66C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI 67C----------------------------------------------------------------------- 68 ZNR = ZI 69 ZNI = -ZR 70 ZBR = ZR 71 ZBI = ZI 72 CIDI = -CONER 73 INU = FNU 74 ANG = HPI*(FNU-INU) 75 C2R = COS(ANG) 76 C2I = SIN(ANG) 77 CAR = C2R 78 SAR = C2I 79 IN = INU + N - 1 80 IN = MOD(IN,4) + 1 81 STR = C2R*CIPR(IN) - C2I*CIPI(IN) 82 C2I = C2R*CIPI(IN) + C2I*CIPR(IN) 83 C2R = STR 84 IF (ZI.GT.0.0D0) GO TO 10 85 ZNR = -ZNR 86 ZBI = -ZBI 87 CIDI = -CIDI 88 C2I = -C2I 89 10 CONTINUE 90C----------------------------------------------------------------------- 91C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER 92C----------------------------------------------------------------------- 93 FN = MAX(FNU,1.0D0) 94 CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, 95 * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 96 IF (KODE.EQ.1) GO TO 20 97 STR = ZBR + ZETA2R 98 STI = ZBI + ZETA2I 99 RAST = FN/ZABS(STR,STI) 100 STR = STR*RAST*RAST 101 STI = -STI*RAST*RAST 102 S1R = -ZETA1R + STR 103 S1I = -ZETA1I + STI 104 GO TO 30 105 20 CONTINUE 106 S1R = -ZETA1R + ZETA2R 107 S1I = -ZETA1I + ZETA2I 108 30 CONTINUE 109 RS1 = S1R 110 IF (ABS(RS1).GT.ELIM) GO TO 150 111 40 CONTINUE 112 NN = MIN(2,ND) 113 DO 90 I=1,NN 114 FN = FNU + (ND-I) 115 CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, 116 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 117 IF (KODE.EQ.1) GO TO 50 118 STR = ZBR + ZETA2R 119 STI = ZBI + ZETA2I 120 RAST = FN/ZABS(STR,STI) 121 STR = STR*RAST*RAST 122 STI = -STI*RAST*RAST 123 S1R = -ZETA1R + STR 124 S1I = -ZETA1I + STI + ABS(ZI) 125 GO TO 60 126 50 CONTINUE 127 S1R = -ZETA1R + ZETA2R 128 S1I = -ZETA1I + ZETA2I 129 60 CONTINUE 130C----------------------------------------------------------------------- 131C TEST FOR UNDERFLOW AND OVERFLOW 132C----------------------------------------------------------------------- 133 RS1 = S1R 134 IF (ABS(RS1).GT.ELIM) GO TO 120 135 IF (I.EQ.1) IFLAG = 2 136 IF (ABS(RS1).LT.ALIM) GO TO 70 137C----------------------------------------------------------------------- 138C REFINE TEST AND SCALE 139C----------------------------------------------------------------------- 140C----------------------------------------------------------------------- 141 APHI = ZABS(PHIR,PHII) 142 AARG = ZABS(ARGR,ARGI) 143 RS1 = RS1 + LOG(APHI) - 0.25D0*LOG(AARG) - AIC 144 IF (ABS(RS1).GT.ELIM) GO TO 120 145 IF (I.EQ.1) IFLAG = 1 146 IF (RS1.LT.0.0D0) GO TO 70 147 IF (I.EQ.1) IFLAG = 3 148 70 CONTINUE 149C----------------------------------------------------------------------- 150C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 151C EXPONENT EXTREMES 152C----------------------------------------------------------------------- 153 CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM) 154 CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM) 155 STR = DAIR*BSUMR - DAII*BSUMI 156 STI = DAIR*BSUMI + DAII*BSUMR 157 STR = STR + (AIR*ASUMR-AII*ASUMI) 158 STI = STI + (AIR*ASUMI+AII*ASUMR) 159 S2R = PHIR*STR - PHII*STI 160 S2I = PHIR*STI + PHII*STR 161 STR = EXP(S1R)*CSSR(IFLAG) 162 S1R = STR*COS(S1I) 163 S1I = STR*SIN(S1I) 164 STR = S2R*S1R - S2I*S1I 165 S2I = S2R*S1I + S2I*S1R 166 S2R = STR 167 IF (IFLAG.NE.1) GO TO 80 168 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 169 IF (NW.NE.0) GO TO 120 170 80 CONTINUE 171 IF (ZI.LE.0.0D0) S2I = -S2I 172 STR = S2R*C2R - S2I*C2I 173 S2I = S2R*C2I + S2I*C2R 174 S2R = STR 175 CYR(I) = S2R 176 CYI(I) = S2I 177 J = ND - I + 1 178 YR(J) = S2R*CSRR(IFLAG) 179 YI(J) = S2I*CSRR(IFLAG) 180 STR = -C2I*CIDI 181 C2I = C2R*CIDI 182 C2R = STR 183 90 CONTINUE 184 IF (ND.LE.2) GO TO 110 185 RAZ = 1.0D0/ZABS(ZR,ZI) 186 STR = ZR*RAZ 187 STI = -ZI*RAZ 188 RZR = (STR+STR)*RAZ 189 RZI = (STI+STI)*RAZ 190 BRY(2) = 1.0D0/BRY(1) 191 BRY(3) = D1MACH(2) 192 S1R = CYR(1) 193 S1I = CYI(1) 194 S2R = CYR(2) 195 S2I = CYI(2) 196 C1R = CSRR(IFLAG) 197 ASCLE = BRY(IFLAG) 198 K = ND - 2 199 FN = K 200 DO 100 I=3,ND 201 C2R = S2R 202 C2I = S2I 203 S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) 204 S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) 205 S1R = C2R 206 S1I = C2I 207 C2R = S2R*C1R 208 C2I = S2I*C1R 209 YR(K) = C2R 210 YI(K) = C2I 211 K = K - 1 212 FN = FN - 1.0D0 213 IF (IFLAG.GE.3) GO TO 100 214 STR = ABS(C2R) 215 STI = ABS(C2I) 216 C2M = MAX(STR,STI) 217 IF (C2M.LE.ASCLE) GO TO 100 218 IFLAG = IFLAG + 1 219 ASCLE = BRY(IFLAG) 220 S1R = S1R*C1R 221 S1I = S1I*C1R 222 S2R = C2R 223 S2I = C2I 224 S1R = S1R*CSSR(IFLAG) 225 S1I = S1I*CSSR(IFLAG) 226 S2R = S2R*CSSR(IFLAG) 227 S2I = S2I*CSSR(IFLAG) 228 C1R = CSRR(IFLAG) 229 100 CONTINUE 230 110 CONTINUE 231 RETURN 232 120 CONTINUE 233 IF (RS1.GT.0.0D0) GO TO 140 234C----------------------------------------------------------------------- 235C SET UNDERFLOW AND UPDATE PARAMETERS 236C----------------------------------------------------------------------- 237 YR(ND) = ZEROR 238 YI(ND) = ZEROI 239 NZ = NZ + 1 240 ND = ND - 1 241 IF (ND.EQ.0) GO TO 110 242 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) 243 IF (NUF.LT.0) GO TO 140 244 ND = ND - NUF 245 NZ = NZ + NUF 246 IF (ND.EQ.0) GO TO 110 247 FN = FNU + (ND-1) 248 IF (FN.LT.FNUL) GO TO 130 249C FN = CIDI 250C J = NUF + 1 251C K = MOD(J,4) + 1 252C S1R = CIPR(K) 253C S1I = CIPI(K) 254C IF (FN.LT.0.0D0) S1I = -S1I 255C STR = C2R*S1R - C2I*S1I 256C C2I = C2R*S1I + C2I*S1R 257C C2R = STR 258 IN = INU + ND - 1 259 IN = MOD(IN,4) + 1 260 C2R = CAR*CIPR(IN) - SAR*CIPI(IN) 261 C2I = CAR*CIPI(IN) + SAR*CIPR(IN) 262 IF (ZI.LE.0.0D0) C2I = -C2I 263 GO TO 40 264 130 CONTINUE 265 NLAST = ND 266 RETURN 267 140 CONTINUE 268 NZ = -1 269 RETURN 270 150 CONTINUE 271 IF (RS1.GT.0.0D0) GO TO 140 272 NZ = N 273 DO 160 I=1,N 274 YR(I) = ZEROR 275 YI(I) = ZEROI 276 160 CONTINUE 277 RETURN 278 END 279