src/numerical/slatec/fortran/zuoik.f

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