1*DECK ZUOIK 2 SUBROUTINE ZUOIK (ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, 3 + ELIM, ALIM) 4C***BEGIN PROLOGUE ZUOIK 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to ZBESH, ZBESI and ZBESK 7C***LIBRARY SLATEC 8C***TYPE ALL (CUOIK-A, ZUOIK-A) 9C***AUTHOR Amos, D. E., (SNL) 10C***DESCRIPTION 11C 12C ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC 13C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM 14C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW 15C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING 16C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN 17C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER 18C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE 19C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)= 20C EXP(-ELIM)/TOL 21C 22C IKFLG=1 MEANS THE I SEQUENCE IS TESTED 23C =2 MEANS THE K SEQUENCE IS TESTED 24C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE 25C =-1 MEANS AN OVERFLOW WOULD OCCUR 26C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO 27C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE 28C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO 29C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY 30C ANOTHER ROUTINE 31C 32C***SEE ALSO ZBESH, ZBESI, ZBESK 33C***ROUTINES CALLED D1MACH, ZABS, ZLOG, ZUCHK, ZUNHJ, ZUNIK 34C***REVISION HISTORY (YYMMDD) 35C 830501 DATE WRITTEN 36C 910415 Prologue converted to Version 4.0 format. (BAB) 37C 930122 Added ZLOG to EXTERNAL statement. (RWC) 38C***END PROLOGUE ZUOIK 39C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN, 40C *ZR 41 DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR, 42 * ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN, 43 * FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI, 44 * YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, 45 * ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS 46 INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW 47 DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16) 48 EXTERNAL ZABS, ZLOG 49 DATA ZEROR,ZEROI / 0.0D0, 0.0D0 / 50 DATA AIC / 1.265512123484645396D+00 / 51C***FIRST EXECUTABLE STATEMENT ZUOIK 52 NUF = 0 53 NN = N 54 ZRR = ZR 55 ZRI = ZI 56 IF (ZR.GE.0.0D0) GO TO 10 57 ZRR = -ZR 58 ZRI = -ZI 59 10 CONTINUE 60 ZBR = ZRR 61 ZBI = ZRI 62 AX = ABS(ZR)*1.7321D0 63 AY = ABS(ZI) 64 IFORM = 1 65 IF (AY.GT.AX) IFORM = 2 66 GNU = MAX(FNU,1.0D0) 67 IF (IKFLG.EQ.1) GO TO 20 68 FNN = NN 69 GNN = FNU + FNN - 1.0D0 70 GNU = MAX(GNN,FNN) 71 20 CONTINUE 72C----------------------------------------------------------------------- 73C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE 74C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET 75C THE SIGN OF THE IMAGINARY PART CORRECT. 76C----------------------------------------------------------------------- 77 IF (IFORM.EQ.2) GO TO 30 78 INIT = 0 79 CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, 80 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) 81 CZR = -ZETA1R + ZETA2R 82 CZI = -ZETA1I + ZETA2I 83 GO TO 50 84 30 CONTINUE 85 ZNR = ZRI 86 ZNI = -ZRR 87 IF (ZI.GT.0.0D0) GO TO 40 88 ZNR = -ZNR 89 40 CONTINUE 90 CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, 91 * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 92 CZR = -ZETA1R + ZETA2R 93 CZI = -ZETA1I + ZETA2I 94 AARG = ZABS(ARGR,ARGI) 95 50 CONTINUE 96 IF (KODE.EQ.1) GO TO 60 97 CZR = CZR - ZBR 98 CZI = CZI - ZBI 99 60 CONTINUE 100 IF (IKFLG.EQ.1) GO TO 70 101 CZR = -CZR 102 CZI = -CZI 103 70 CONTINUE 104 APHI = ZABS(PHIR,PHII) 105 RCZ = CZR 106C----------------------------------------------------------------------- 107C OVERFLOW TEST 108C----------------------------------------------------------------------- 109 IF (RCZ.GT.ELIM) GO TO 210 110 IF (RCZ.LT.ALIM) GO TO 80 111 RCZ = RCZ + LOG(APHI) 112 IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC 113 IF (RCZ.GT.ELIM) GO TO 210 114 GO TO 130 115 80 CONTINUE 116C----------------------------------------------------------------------- 117C UNDERFLOW TEST 118C----------------------------------------------------------------------- 119 IF (RCZ.LT.(-ELIM)) GO TO 90 120 IF (RCZ.GT.(-ALIM)) GO TO 130 121 RCZ = RCZ + LOG(APHI) 122 IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC 123 IF (RCZ.GT.(-ELIM)) GO TO 110 124 90 CONTINUE 125 DO 100 I=1,NN 126 YR(I) = ZEROR 127 YI(I) = ZEROI 128 100 CONTINUE 129 NUF = NN 130 RETURN 131 110 CONTINUE 132 ASCLE = 1.0D+3*D1MACH(1)/TOL 133 CALL ZLOG(PHIR, PHII, STR, STI, IDUM) 134 CZR = CZR + STR 135 CZI = CZI + STI 136 IF (IFORM.EQ.1) GO TO 120 137 CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) 138 CZR = CZR - 0.25D0*STR - AIC 139 CZI = CZI - 0.25D0*STI 140 120 CONTINUE 141 AX = EXP(RCZ)/TOL 142 AY = CZI 143 CZR = AX*COS(AY) 144 CZI = AX*SIN(AY) 145 CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL) 146 IF (NW.NE.0) GO TO 90 147 130 CONTINUE 148 IF (IKFLG.EQ.2) RETURN 149 IF (N.EQ.1) RETURN 150C----------------------------------------------------------------------- 151C SET UNDERFLOWS ON I SEQUENCE 152C----------------------------------------------------------------------- 153 140 CONTINUE 154 GNU = FNU + (NN-1) 155 IF (IFORM.EQ.2) GO TO 150 156 INIT = 0 157 CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, 158 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) 159 CZR = -ZETA1R + ZETA2R 160 CZI = -ZETA1I + ZETA2I 161 GO TO 160 162 150 CONTINUE 163 CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, 164 * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 165 CZR = -ZETA1R + ZETA2R 166 CZI = -ZETA1I + ZETA2I 167 AARG = ZABS(ARGR,ARGI) 168 160 CONTINUE 169 IF (KODE.EQ.1) GO TO 170 170 CZR = CZR - ZBR 171 CZI = CZI - ZBI 172 170 CONTINUE 173 APHI = ZABS(PHIR,PHII) 174 RCZ = CZR 175 IF (RCZ.LT.(-ELIM)) GO TO 180 176 IF (RCZ.GT.(-ALIM)) RETURN 177 RCZ = RCZ + LOG(APHI) 178 IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC 179 IF (RCZ.GT.(-ELIM)) GO TO 190 180 180 CONTINUE 181 YR(NN) = ZEROR 182 YI(NN) = ZEROI 183 NN = NN - 1 184 NUF = NUF + 1 185 IF (NN.EQ.0) RETURN 186 GO TO 140 187 190 CONTINUE 188 ASCLE = 1.0D+3*D1MACH(1)/TOL 189 CALL ZLOG(PHIR, PHII, STR, STI, IDUM) 190 CZR = CZR + STR 191 CZI = CZI + STI 192 IF (IFORM.EQ.1) GO TO 200 193 CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) 194 CZR = CZR - 0.25D0*STR - AIC 195 CZI = CZI - 0.25D0*STI 196 200 CONTINUE 197 AX = EXP(RCZ)/TOL 198 AY = CZI 199 CZR = AX*COS(AY) 200 CZI = AX*SIN(AY) 201 CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL) 202 IF (NW.NE.0) GO TO 180 203 RETURN 204 210 CONTINUE 205 NUF = -1 206 RETURN 207 END 208