amos/zuoik.f

   1       SUBROUTINE ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
   2      * ELIM, ALIM)
   3 C***BEGIN PROLOGUE  ZUOIK
   4 C***REFER TO  ZBESI,ZBESK,ZBESH
   5 C
   6 C     ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
   7 C     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
   8 C     (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
   9 C     WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
  10 C     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
  11 C     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
  12 C     MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
  13 C     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
  14 C     EXP(-ELIM)/TOL
  15 C
  16 C     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
  17 C          =2 MEANS THE K SEQUENCE IS TESTED
  18 C     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
  19 C         =-1 MEANS AN OVERFLOW WOULD OCCUR
  20 C     IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
  21 C             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
  22 C     IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
  23 C     IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
  24 C             ANOTHER ROUTINE
  25 C
  26 C***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,AZABS,AZLOG
  27 C***END PROLOGUE  ZUOIK
  28 C     COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
  29 C    *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
  59 C-----------------------------------------------------------------------
  60 C     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
  61 C     REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
  62 C     THE SIGN OF THE IMAGINARY PART CORRECT.
  63 C-----------------------------------------------------------------------
  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
  93 C-----------------------------------------------------------------------
  94 C     OVERFLOW TEST
  95 C-----------------------------------------------------------------------
  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
 103 C-----------------------------------------------------------------------
 104 C     UNDERFLOW TEST
 105 C-----------------------------------------------------------------------
 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
 137 C-----------------------------------------------------------------------
 138 C     SET UNDERFLOWS ON I SEQUENCE
 139 C-----------------------------------------------------------------------
 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