1 SUBROUTINE CUNI2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, 2 * ALIM) 3C***BEGIN PROLOGUE CUNI2 4C***REFER TO CBESI,CBESK 5C 6C CUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF 7C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I 8C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. 9C 10C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC 11C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. 12C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER 13C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. 14C Y(I)=CZERO FOR I=NLAST+1,N 15C 16C***ROUTINES CALLED CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH 17C***END PROLOGUE CUNI2 18 COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL, 19 * CSR, CSS, CY, CZERO, C1, C2, DAI, PHI, RZ, S1, S2, Y, Z, ZB, 20 * ZETA1, ZETA2, ZN, ZAR 21 REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M, 22 * C2R, ELIM, FN, FNU, FNUL, HPI, RS1, SAR, TOL, YY, R1MACH 23 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, 24 * NN, NUF, NW, NZ, IDUM 25 DIMENSION BRY(3), Y(N), CIP(4), CSS(3), CSR(3), CY(2) 26 DATA CZERO,CONE,CI/(0.0E0,0.0E0),(1.0E0,0.0E0),(0.0E0,1.0E0)/ 27 DATA CIP(1),CIP(2),CIP(3),CIP(4)/ 28 1 (1.0E0,0.0E0), (0.0E0,1.0E0), (-1.0E0,0.0E0), (0.0E0,-1.0E0)/ 29 DATA HPI, AIC / 30 1 1.57079632679489662E+00, 1.265512123484645396E+00/ 31C 32 NZ = 0 33 ND = N 34 NLAST = 0 35C----------------------------------------------------------------------- 36C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- 37C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, 38C EXP(ALIM)=EXP(ELIM)*TOL 39C----------------------------------------------------------------------- 40 CSCL = CMPLX(1.0E0/TOL,0.0E0) 41 CRSC = CMPLX(TOL,0.0E0) 42 CSS(1) = CSCL 43 CSS(2) = CONE 44 CSS(3) = CRSC 45 CSR(1) = CRSC 46 CSR(2) = CONE 47 CSR(3) = CSCL 48 BRY(1) = 1.0E+3*R1MACH(1)/TOL 49 YY = AIMAG(Z) 50C----------------------------------------------------------------------- 51C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI 52C----------------------------------------------------------------------- 53 ZN = -Z*CI 54 ZB = Z 55 CID = -CI 56 INU = INT(FNU) 57 ANG = HPI*(FNU-FLOAT(INU)) 58 CAR = COS(ANG) 59 SAR = SIN(ANG) 60 C2 = CMPLX(CAR,SAR) 61 ZAR = C2 62 IN = INU + N - 1 63 IN = MOD(IN,4) 64 C2 = C2*CIP(IN+1) 65 IF (YY.GT.0.0E0) GO TO 10 66 ZN = CONJG(-ZN) 67 ZB = CONJG(ZB) 68 CID = -CID 69 C2 = CONJG(C2) 70 10 CONTINUE 71C----------------------------------------------------------------------- 72C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER 73C----------------------------------------------------------------------- 74 FN = AMAX1(FNU,1.0E0) 75 CALL CUNHJ(ZN, FN, 1, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM) 76 IF (KODE.EQ.1) GO TO 20 77 CFN = CMPLX(FNU,0.0E0) 78 S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2)) 79 GO TO 30 80 20 CONTINUE 81 S1 = -ZETA1 + ZETA2 82 30 CONTINUE 83 RS1 = REAL(S1) 84 IF (ABS(RS1).GT.ELIM) GO TO 150 85 40 CONTINUE 86 NN = MIN0(2,ND) 87 DO 90 I=1,NN 88 FN = FNU + FLOAT(ND-I) 89 CALL CUNHJ(ZN, FN, 0, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM) 90 IF (KODE.EQ.1) GO TO 50 91 CFN = CMPLX(FN,0.0E0) 92 AY = ABS(YY) 93 S1 = -ZETA1 + CFN*(CFN/(ZB+ZETA2)) + CMPLX(0.0E0,AY) 94 GO TO 60 95 50 CONTINUE 96 S1 = -ZETA1 + ZETA2 97 60 CONTINUE 98C----------------------------------------------------------------------- 99C TEST FOR UNDERFLOW AND OVERFLOW 100C----------------------------------------------------------------------- 101 RS1 = REAL(S1) 102 IF (ABS(RS1).GT.ELIM) GO TO 120 103 IF (I.EQ.1) IFLAG = 2 104 IF (ABS(RS1).LT.ALIM) GO TO 70 105C----------------------------------------------------------------------- 106C REFINE TEST AND SCALE 107C----------------------------------------------------------------------- 108C----------------------------------------------------------------------- 109 APHI = CABS(PHI) 110 AARG = CABS(ARG) 111 RS1 = RS1 + ALOG(APHI) - 0.25E0*ALOG(AARG) - AIC 112 IF (ABS(RS1).GT.ELIM) GO TO 120 113 IF (I.EQ.1) IFLAG = 1 114 IF (RS1.LT.0.0E0) GO TO 70 115 IF (I.EQ.1) IFLAG = 3 116 70 CONTINUE 117C----------------------------------------------------------------------- 118C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 119C EXPONENT EXTREMES 120C----------------------------------------------------------------------- 121 CALL CAIRY(ARG, 0, 2, AI, NAI, IDUM) 122 CALL CAIRY(ARG, 1, 2, DAI, NDAI, IDUM) 123 S2 = PHI*(AI*ASUM+DAI*BSUM) 124 C2R = REAL(S1) 125 C2I = AIMAG(S1) 126 C2M = EXP(C2R)*REAL(CSS(IFLAG)) 127 S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I)) 128 S2 = S2*S1 129 IF (IFLAG.NE.1) GO TO 80 130 CALL CUCHK(S2, NW, BRY(1), TOL) 131 IF (NW.NE.0) GO TO 120 132 80 CONTINUE 133 IF (YY.LE.0.0E0) S2 = CONJG(S2) 134 J = ND - I + 1 135 S2 = S2*C2 136 CY(I) = S2 137 Y(J) = S2*CSR(IFLAG) 138 C2 = C2*CID 139 90 CONTINUE 140 IF (ND.LE.2) GO TO 110 141 RZ = CMPLX(2.0E0,0.0E0)/Z 142 BRY(2) = 1.0E0/BRY(1) 143 BRY(3) = R1MACH(2) 144 S1 = CY(1) 145 S2 = CY(2) 146 C1 = CSR(IFLAG) 147 ASCLE = BRY(IFLAG) 148 K = ND - 2 149 FN = FLOAT(K) 150 DO 100 I=3,ND 151 C2 = S2 152 S2 = S1 + CMPLX(FNU+FN,0.0E0)*RZ*S2 153 S1 = C2 154 C2 = S2*C1 155 Y(K) = C2 156 K = K - 1 157 FN = FN - 1.0E0 158 IF (IFLAG.GE.3) GO TO 100 159 C2R = REAL(C2) 160 C2I = AIMAG(C2) 161 C2R = ABS(C2R) 162 C2I = ABS(C2I) 163 C2M = AMAX1(C2R,C2I) 164 IF (C2M.LE.ASCLE) GO TO 100 165 IFLAG = IFLAG + 1 166 ASCLE = BRY(IFLAG) 167 S1 = S1*C1 168 S2 = C2 169 S1 = S1*CSS(IFLAG) 170 S2 = S2*CSS(IFLAG) 171 C1 = CSR(IFLAG) 172 100 CONTINUE 173 110 CONTINUE 174 RETURN 175 120 CONTINUE 176 IF (RS1.GT.0.0E0) GO TO 140 177C----------------------------------------------------------------------- 178C SET UNDERFLOW AND UPDATE PARAMETERS 179C----------------------------------------------------------------------- 180 Y(ND) = CZERO 181 NZ = NZ + 1 182 ND = ND - 1 183 IF (ND.EQ.0) GO TO 110 184 CALL CUOIK(Z, FNU, KODE, 1, ND, Y, NUF, TOL, ELIM, ALIM) 185 IF (NUF.LT.0) GO TO 140 186 ND = ND - NUF 187 NZ = NZ + NUF 188 IF (ND.EQ.0) GO TO 110 189 FN = FNU + FLOAT(ND-1) 190 IF (FN.LT.FNUL) GO TO 130 191C FN = AIMAG(CID) 192C J = NUF + 1 193C K = MOD(J,4) + 1 194C S1 = CIP(K) 195C IF (FN.LT.0.0E0) S1 = CONJG(S1) 196C C2 = C2*S1 197 IN = INU + ND - 1 198 IN = MOD(IN,4) + 1 199 C2 = ZAR*CIP(IN) 200 IF (YY.LE.0.0E0)C2=CONJG(C2) 201 GO TO 40 202 130 CONTINUE 203 NLAST = ND 204 RETURN 205 140 CONTINUE 206 NZ = -1 207 RETURN 208 150 CONTINUE 209 IF (RS1.GT.0.0E0) GO TO 140 210 NZ = N 211 DO 160 I=1,N 212 Y(I) = CZERO 213 160 CONTINUE 214 RETURN 215 END 216