SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, * TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZUNI2 C***REFER TO ZBESI,ZBESK C C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. C C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. C Y(I)=CZERO FOR I=NLAST+1,N C C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,AZABS C***END PROLOGUE ZUNI2 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI, * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR, * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII, * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, * CYI, D1MACH, AZABS, CAR, SAR INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, * NN, NUF, NW, NZ, IDUM DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), * CSRR(3), CYR(2), CYI(2) DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/ DATA HPI, AIC / 1 1.57079632679489662D+00, 1.265512123484645396D+00/ C NZ = 0 ND = N NLAST = 0 C----------------------------------------------------------------------- C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, C EXP(ALIM)=EXP(ELIM)*TOL C----------------------------------------------------------------------- CSCL = 1.0D0/TOL CRSC = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CRSC CSRR(1) = CRSC CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = 1.0D+3*D1MACH(1)/TOL C----------------------------------------------------------------------- C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI C----------------------------------------------------------------------- ZNR = ZI ZNI = -ZR ZBR = ZR ZBI = ZI CIDI = -CONER INU = INT(SNGL(FNU)) ANG = HPI*(FNU-DBLE(FLOAT(INU))) C2R = DCOS(ANG) C2I = DSIN(ANG) CAR = C2R SAR = C2I IN = INU + N - 1 IN = MOD(IN,4) + 1 STR = C2R*CIPR(IN) - C2I*CIPI(IN) C2I = C2R*CIPI(IN) + C2I*CIPR(IN) C2R = STR IF (ZI.GT.0.0D0) GO TO 10 ZNR = -ZNR ZBI = -ZBI CIDI = -CIDI C2I = -C2I 10 CONTINUE C----------------------------------------------------------------------- C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER C----------------------------------------------------------------------- FN = DMAX1(FNU,1.0D0) CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) IF (KODE.EQ.1) GO TO 20 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI GO TO 30 20 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 30 CONTINUE RS1 = S1R IF (DABS(RS1).GT.ELIM) GO TO 150 40 CONTINUE NN = MIN0(2,ND) DO 90 I=1,NN FN = FNU + DBLE(FLOAT(ND-I)) CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) IF (KODE.EQ.1) GO TO 50 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN/AZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI + DABS(ZI) GO TO 60 50 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 60 CONTINUE C----------------------------------------------------------------------- C TEST FOR UNDERFLOW AND OVERFLOW C----------------------------------------------------------------------- RS1 = S1R IF (DABS(RS1).GT.ELIM) GO TO 120 IF (I.EQ.1) IFLAG = 2 IF (DABS(RS1).LT.ALIM) GO TO 70 C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- C----------------------------------------------------------------------- APHI = AZABS(PHIR,PHII) AARG = AZABS(ARGR,ARGI) RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC IF (DABS(RS1).GT.ELIM) GO TO 120 IF (I.EQ.1) IFLAG = 1 IF (RS1.LT.0.0D0) GO TO 70 IF (I.EQ.1) IFLAG = 3 70 CONTINUE C----------------------------------------------------------------------- C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR C EXPONENT EXTREMES C----------------------------------------------------------------------- CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM) CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM) STR = DAIR*BSUMR - DAII*BSUMI STI = DAIR*BSUMI + DAII*BSUMR STR = STR + (AIR*ASUMR-AII*ASUMI) STI = STI + (AIR*ASUMI+AII*ASUMR) S2R = PHIR*STR - PHII*STI S2I = PHIR*STI + PHII*STR STR = DEXP(S1R)*CSSR(IFLAG) S1R = STR*DCOS(S1I) S1I = STR*DSIN(S1I) STR = S2R*S1R - S2I*S1I S2I = S2R*S1I + S2I*S1R S2R = STR IF (IFLAG.NE.1) GO TO 80 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) IF (NW.NE.0) GO TO 120 80 CONTINUE IF (ZI.LE.0.0D0) S2I = -S2I STR = S2R*C2R - S2I*C2I S2I = S2R*C2I + S2I*C2R S2R = STR CYR(I) = S2R CYI(I) = S2I J = ND - I + 1 YR(J) = S2R*CSRR(IFLAG) YI(J) = S2I*CSRR(IFLAG) STR = -C2I*CIDI C2I = C2R*CIDI C2R = STR 90 CONTINUE IF (ND.LE.2) GO TO 110 RAZ = 1.0D0/AZABS(ZR,ZI) STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ BRY(2) = 1.0D0/BRY(1) BRY(3) = D1MACH(2) S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) C1R = CSRR(IFLAG) ASCLE = BRY(IFLAG) K = ND - 2 FN = DBLE(FLOAT(K)) DO 100 I=3,ND C2R = S2R C2I = S2I S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) S1R = C2R S1I = C2I C2R = S2R*C1R C2I = S2I*C1R YR(K) = C2R YI(K) = C2I K = K - 1 FN = FN - 1.0D0 IF (IFLAG.GE.3) GO TO 100 STR = DABS(C2R) STI = DABS(C2I) C2M = DMAX1(STR,STI) IF (C2M.LE.ASCLE) GO TO 100 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R*C1R S1I = S1I*C1R S2R = C2R S2I = C2I S1R = S1R*CSSR(IFLAG) S1I = S1I*CSSR(IFLAG) S2R = S2R*CSSR(IFLAG) S2I = S2I*CSSR(IFLAG) C1R = CSRR(IFLAG) 100 CONTINUE 110 CONTINUE RETURN 120 CONTINUE IF (RS1.GT.0.0D0) GO TO 140 C----------------------------------------------------------------------- C SET UNDERFLOW AND UPDATE PARAMETERS C----------------------------------------------------------------------- YR(ND) = ZEROR YI(ND) = ZEROI NZ = NZ + 1 ND = ND - 1 IF (ND.EQ.0) GO TO 110 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) IF (NUF.LT.0) GO TO 140 ND = ND - NUF NZ = NZ + NUF IF (ND.EQ.0) GO TO 110 FN = FNU + DBLE(FLOAT(ND-1)) IF (FN.LT.FNUL) GO TO 130 C FN = CIDI C J = NUF + 1 C K = MOD(J,4) + 1 C S1R = CIPR(K) C S1I = CIPI(K) C IF (FN.LT.0.0D0) S1I = -S1I C STR = C2R*S1R - C2I*S1I C C2I = C2R*S1I + C2I*S1R C C2R = STR IN = INU + ND - 1 IN = MOD(IN,4) + 1 C2R = CAR*CIPR(IN) - SAR*CIPI(IN) C2I = CAR*CIPI(IN) + SAR*CIPR(IN) IF (ZI.LE.0.0D0) C2I = -C2I GO TO 40 130 CONTINUE NLAST = ND RETURN 140 CONTINUE NZ = -1 RETURN 150 CONTINUE IF (RS1.GT.0.0D0) GO TO 140 NZ = N DO 160 I=1,N YR(I) = ZEROR YI(I) = ZEROI 160 CONTINUE RETURN END