1 SUBROUTINE CUNK1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM) 2C***BEGIN PROLOGUE CUNK1 3C***REFER TO CBESK 4C 5C CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE 6C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE 7C UNIFORM ASYMPTOTIC EXPANSION. 8C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. 9C NZ=-1 MEANS AN OVERFLOW WILL OCCUR 10C 11C***ROUTINES CALLED CS1S2,CUCHK,CUNIK,R1MACH 12C***END PROLOGUE CUNK1 13 COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS, 14 * CWRK, CY, CZERO, C1, C2, PHI, RZ, SUM, S1, S2, Y, Z, 15 * ZETA1, ZETA2, ZR, PHID, ZETA1D, ZETA2D, SUMD 16 REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM, 17 * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH 18 INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG, 19 * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC 20 DIMENSION BRY(3), INIT(2), Y(N), SUM(2), PHI(2), ZETA1(2), 21 * ZETA2(2), CY(2), CWRK(16,3), CSS(3), CSR(3) 22 DATA CZERO, CONE / (0.0E0,0.0E0) , (1.0E0,0.0E0) / 23 DATA PI / 3.14159265358979324E0 / 24C 25 KDFLG = 1 26 NZ = 0 27C----------------------------------------------------------------------- 28C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN 29C THE UNDERFLOW LIMIT 30C----------------------------------------------------------------------- 31 CSCL = CMPLX(1.0E0/TOL,0.0E0) 32 CRSC = CMPLX(TOL,0.0E0) 33 CSS(1) = CSCL 34 CSS(2) = CONE 35 CSS(3) = CRSC 36 CSR(1) = CRSC 37 CSR(2) = CONE 38 CSR(3) = CSCL 39 BRY(1) = 1.0E+3*R1MACH(1)/TOL 40 BRY(2) = 1.0E0/BRY(1) 41 BRY(3) = R1MACH(2) 42 X = REAL(Z) 43 ZR = Z 44 IF (X.LT.0.0E0) ZR = -Z 45 J=2 46 DO 70 I=1,N 47C----------------------------------------------------------------------- 48C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J 49C----------------------------------------------------------------------- 50 J = 3 - J 51 FN = FNU + FLOAT(I-1) 52 INIT(J) = 0 53 CALL CUNIK(ZR, FN, 2, 0, TOL, INIT(J), PHI(J), ZETA1(J), 54 * ZETA2(J), SUM(J), CWRK(1,J)) 55 IF (KODE.EQ.1) GO TO 20 56 CFN = CMPLX(FN,0.0E0) 57 S1 = ZETA1(J) - CFN*(CFN/(ZR+ZETA2(J))) 58 GO TO 30 59 20 CONTINUE 60 S1 = ZETA1(J) - ZETA2(J) 61 30 CONTINUE 62C----------------------------------------------------------------------- 63C TEST FOR UNDERFLOW AND OVERFLOW 64C----------------------------------------------------------------------- 65 RS1 = REAL(S1) 66 IF (ABS(RS1).GT.ELIM) GO TO 60 67 IF (KDFLG.EQ.1) KFLAG = 2 68 IF (ABS(RS1).LT.ALIM) GO TO 40 69C----------------------------------------------------------------------- 70C REFINE TEST AND SCALE 71C----------------------------------------------------------------------- 72 APHI = CABS(PHI(J)) 73 RS1 = RS1 + ALOG(APHI) 74 IF (ABS(RS1).GT.ELIM) GO TO 60 75 IF (KDFLG.EQ.1) KFLAG = 1 76 IF (RS1.LT.0.0E0) GO TO 40 77 IF (KDFLG.EQ.1) KFLAG = 3 78 40 CONTINUE 79C----------------------------------------------------------------------- 80C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 81C EXPONENT EXTREMES 82C----------------------------------------------------------------------- 83 S2 = PHI(J)*SUM(J) 84 C2R = REAL(S1) 85 C2I = AIMAG(S1) 86 C2M = EXP(C2R)*REAL(CSS(KFLAG)) 87 S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I)) 88 S2 = S2*S1 89 IF (KFLAG.NE.1) GO TO 50 90 CALL CUCHK(S2, NW, BRY(1), TOL) 91 IF (NW.NE.0) GO TO 60 92 50 CONTINUE 93 CY(KDFLG) = S2 94 Y(I) = S2*CSR(KFLAG) 95 IF (KDFLG.EQ.2) GO TO 75 96 KDFLG = 2 97 GO TO 70 98 60 CONTINUE 99 IF (RS1.GT.0.0E0) GO TO 290 100C----------------------------------------------------------------------- 101C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 102C----------------------------------------------------------------------- 103 IF (X.LT.0.0E0) GO TO 290 104 KDFLG = 1 105 Y(I) = CZERO 106 NZ=NZ+1 107 IF (I.EQ.1) GO TO 70 108 IF (Y(I-1).EQ.CZERO) GO TO 70 109 Y(I-1) = CZERO 110 NZ=NZ+1 111 70 CONTINUE 112 I=N 113 75 CONTINUE 114 RZ = CMPLX(2.0E0,0.0E0)/ZR 115 CK = CMPLX(FN,0.0E0)*RZ 116 IB = I+1 117 IF (N.LT.IB) GO TO 160 118C----------------------------------------------------------------------- 119C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO 120C ON UNDERFLOW 121C----------------------------------------------------------------------- 122 FN = FNU+FLOAT(N-1) 123 IPARD = 1 124 IF (MR.NE.0) IPARD = 0 125 INITD = 0 126 CALL CUNIK(ZR,FN,2,IPARD,TOL,INITD,PHID,ZETA1D,ZETA2D,SUMD, 127 *CWRK(1,3)) 128 IF (KODE.EQ.1) GO TO 80 129 CFN=CMPLX(FN,0.0E0) 130 S1=ZETA1D-CFN*(CFN/(ZR+ZETA2D)) 131 GO TO 90 132 80 CONTINUE 133 S1=ZETA1D-ZETA2D 134 90 CONTINUE 135 RS1=REAL(S1) 136 IF (ABS(RS1).GT.ELIM) GO TO 95 137 IF (ABS(RS1).LT.ALIM) GO TO 100 138C----------------------------------------------------------------------- 139C REFINE ESTIMATE AND TEST 140C----------------------------------------------------------------------- 141 APHI=CABS(PHID) 142 RS1=RS1+ALOG(APHI) 143 IF (ABS(RS1).LT.ELIM) GO TO 100 144 95 CONTINUE 145 IF (RS1.GT.0.0E0) GO TO 290 146C----------------------------------------------------------------------- 147C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 148C----------------------------------------------------------------------- 149 IF (X.LT.0.0E0) GO TO 290 150 NZ=N 151 DO 96 I=1,N 152 Y(I) = CZERO 153 96 CONTINUE 154 RETURN 155 100 CONTINUE 156C----------------------------------------------------------------------- 157C RECUR FORWARD FOR REMAINDER OF THE SEQUENCE 158C----------------------------------------------------------------------- 159 S1 = CY(1) 160 S2 = CY(2) 161 C1 = CSR(KFLAG) 162 ASCLE = BRY(KFLAG) 163 DO 120 I=IB,N 164 C2 = S2 165 S2 = CK*S2 + S1 166 S1 = C2 167 CK = CK + RZ 168 C2 = S2*C1 169 Y(I) = C2 170 IF (KFLAG.GE.3) GO TO 120 171 C2R = REAL(C2) 172 C2I = AIMAG(C2) 173 C2R = ABS(C2R) 174 C2I = ABS(C2I) 175 C2M = AMAX1(C2R,C2I) 176 IF (C2M.LE.ASCLE) GO TO 120 177 KFLAG = KFLAG + 1 178 ASCLE = BRY(KFLAG) 179 S1 = S1*C1 180 S2 = C2 181 S1 = S1*CSS(KFLAG) 182 S2 = S2*CSS(KFLAG) 183 C1 = CSR(KFLAG) 184 120 CONTINUE 185 160 CONTINUE 186 IF (MR.EQ.0) RETURN 187C----------------------------------------------------------------------- 188C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0 189C----------------------------------------------------------------------- 190 NZ = 0 191 FMR = FLOAT(MR) 192 SGN = -SIGN(PI,FMR) 193C----------------------------------------------------------------------- 194C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP. 195C----------------------------------------------------------------------- 196 CSGN = CMPLX(0.0E0,SGN) 197 INU = INT(FNU) 198 FNF = FNU - FLOAT(INU) 199 IFN = INU + N - 1 200 ANG = FNF*SGN 201 CPN = COS(ANG) 202 SPN = SIN(ANG) 203 CSPN = CMPLX(CPN,SPN) 204 IF (MOD(IFN,2).EQ.1) CSPN = -CSPN 205 ASC = BRY(1) 206 KK = N 207 IUF = 0 208 KDFLG = 1 209 IB = IB-1 210 IC = IB-1 211 DO 260 K=1,N 212 FN = FNU + FLOAT(KK-1) 213C----------------------------------------------------------------------- 214C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K 215C FUNCTION ABOVE 216C----------------------------------------------------------------------- 217 M=3 218 IF (N.GT.2) GO TO 175 219 170 CONTINUE 220 INITD = INIT(J) 221 PHID = PHI(J) 222 ZETA1D = ZETA1(J) 223 ZETA2D = ZETA2(J) 224 SUMD = SUM(J) 225 M = J 226 J = 3 - J 227 GO TO 180 228 175 CONTINUE 229 IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180 230 IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 170 231 INITD = 0 232 180 CONTINUE 233 CALL CUNIK(ZR, FN, 1, 0, TOL, INITD, PHID, ZETA1D, 234 * ZETA2D, SUMD, CWRK(1,M)) 235 IF (KODE.EQ.1) GO TO 190 236 CFN = CMPLX(FN,0.0E0) 237 S1 = -ZETA1D + CFN*(CFN/(ZR+ZETA2D)) 238 GO TO 200 239 190 CONTINUE 240 S1 = -ZETA1D + ZETA2D 241 200 CONTINUE 242C----------------------------------------------------------------------- 243C TEST FOR UNDERFLOW AND OVERFLOW 244C----------------------------------------------------------------------- 245 RS1 = REAL(S1) 246 IF (ABS(RS1).GT.ELIM) GO TO 250 247 IF (KDFLG.EQ.1) IFLAG = 2 248 IF (ABS(RS1).LT.ALIM) GO TO 210 249C----------------------------------------------------------------------- 250C REFINE TEST AND SCALE 251C----------------------------------------------------------------------- 252 APHI = CABS(PHID) 253 RS1 = RS1 + ALOG(APHI) 254 IF (ABS(RS1).GT.ELIM) GO TO 250 255 IF (KDFLG.EQ.1) IFLAG = 1 256 IF (RS1.LT.0.0E0) GO TO 210 257 IF (KDFLG.EQ.1) IFLAG = 3 258 210 CONTINUE 259 S2 = CSGN*PHID*SUMD 260 C2R = REAL(S1) 261 C2I = AIMAG(S1) 262 C2M = EXP(C2R)*REAL(CSS(IFLAG)) 263 S1 = CMPLX(C2M,0.0E0)*CMPLX(COS(C2I),SIN(C2I)) 264 S2 = S2*S1 265 IF (IFLAG.NE.1) GO TO 220 266 CALL CUCHK(S2, NW, BRY(1), TOL) 267 IF (NW.NE.0) S2 = CMPLX(0.0E0,0.0E0) 268 220 CONTINUE 269 CY(KDFLG) = S2 270 C2 = S2 271 S2 = S2*CSR(IFLAG) 272C----------------------------------------------------------------------- 273C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N 274C----------------------------------------------------------------------- 275 S1 = Y(KK) 276 IF (KODE.EQ.1) GO TO 240 277 CALL CS1S2(ZR, S1, S2, NW, ASC, ALIM, IUF) 278 NZ = NZ + NW 279 240 CONTINUE 280 Y(KK) = S1*CSPN + S2 281 KK = KK - 1 282 CSPN = -CSPN 283 IF (C2.NE.CZERO) GO TO 245 284 KDFLG = 1 285 GO TO 260 286 245 CONTINUE 287 IF (KDFLG.EQ.2) GO TO 265 288 KDFLG = 2 289 GO TO 260 290 250 CONTINUE 291 IF (RS1.GT.0.0E0) GO TO 290 292 S2 = CZERO 293 GO TO 220 294 260 CONTINUE 295 K = N 296 265 CONTINUE 297 IL = N - K 298 IF (IL.EQ.0) RETURN 299C----------------------------------------------------------------------- 300C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE 301C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP 302C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. 303C----------------------------------------------------------------------- 304 S1 = CY(1) 305 S2 = CY(2) 306 CS = CSR(IFLAG) 307 ASCLE = BRY(IFLAG) 308 FN = FLOAT(INU+IL) 309 DO 280 I=1,IL 310 C2 = S2 311 S2 = S1 + CMPLX(FN+FNF,0.0E0)*RZ*S2 312 S1 = C2 313 FN = FN - 1.0E0 314 C2 = S2*CS 315 CK = C2 316 C1 = Y(KK) 317 IF (KODE.EQ.1) GO TO 270 318 CALL CS1S2(ZR, C1, C2, NW, ASC, ALIM, IUF) 319 NZ = NZ + NW 320 270 CONTINUE 321 Y(KK) = C1*CSPN + C2 322 KK = KK - 1 323 CSPN = -CSPN 324 IF (IFLAG.GE.3) GO TO 280 325 C2R = REAL(CK) 326 C2I = AIMAG(CK) 327 C2R = ABS(C2R) 328 C2I = ABS(C2I) 329 C2M = AMAX1(C2R,C2I) 330 IF (C2M.LE.ASCLE) GO TO 280 331 IFLAG = IFLAG + 1 332 ASCLE = BRY(IFLAG) 333 S1 = S1*CS 334 S2 = CK 335 S1 = S1*CSS(IFLAG) 336 S2 = S2*CSS(IFLAG) 337 CS = CSR(IFLAG) 338 280 CONTINUE 339 RETURN 340 290 CONTINUE 341 NZ = -1 342 RETURN 343 END 344