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