amos/zuni1.f

   1       SUBROUTINE ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
   2      * TOL, ELIM, ALIM)
   3 C***BEGIN PROLOGUE  ZUNI1
   4 C***REFER TO  ZBESI,ZBESK
   5 C
   6 C     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
   7 C     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
   8 C
   9 C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
  10 C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
  11 C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
  12 C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
  13 C     Y(I)=CZERO FOR I=NLAST+1,N
  14 C
  15 C***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,AZABS
  16 C***END PROLOGUE  ZUNI1
  17 C     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
  18 C    *S2,Y,Z,ZETA1,ZETA2
  19       DOUBLE PRECISION ALIM, APHI, ASCLE, BRY, CONER, CRSC,
  20      * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN,
  21      * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI,
  22      * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I,
  23      * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, AZABS
  24       INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ
  25       DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3),
  26      * CSRR(3), CYR(2), CYI(2)
  27       DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
  28 C
  29       NZ = 0
  30       ND = N
  31       NLAST = 0
  32 C-----------------------------------------------------------------------
  33 C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
  34 C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
  35 C     EXP(ALIM)=EXP(ELIM)*TOL
  36 C-----------------------------------------------------------------------
  37       CSCL = 1.0D0/TOL
  38       CRSC = TOL
  39       CSSR(1) = CSCL
  40       CSSR(2) = CONER
  41       CSSR(3) = CRSC
  42       CSRR(1) = CRSC
  43       CSRR(2) = CONER
  44       CSRR(3) = CSCL
  45       BRY(1) = 1.0D+3*D1MACH(1)/TOL
  46 C-----------------------------------------------------------------------
  47 C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
  48 C-----------------------------------------------------------------------
  49       FN = DMAX1(FNU,1.0D0)
  50       INIT = 0
  51       CALL ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R,
  52      * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
  53       IF (KODE.EQ.1) GO TO 10
  54       STR = ZR + ZETA2R
  55       STI = ZI + ZETA2I
  56       RAST = FN/AZABS(STR,STI)
  57       STR = STR*RAST*RAST
  58       STI = -STI*RAST*RAST
  59       S1R = -ZETA1R + STR
  60       S1I = -ZETA1I + STI
  61       GO TO 20
  62    10 CONTINUE
  63       S1R = -ZETA1R + ZETA2R
  64       S1I = -ZETA1I + ZETA2I
  65    20 CONTINUE
  66       RS1 = S1R
  67       IF (DABS(RS1).GT.ELIM) GO TO 130
  68    30 CONTINUE
  69       NN = MIN0(2,ND)
  70       DO 80 I=1,NN
  71         FN = FNU + DBLE(FLOAT(ND-I))
  72         INIT = 0
  73         CALL ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R,
  74      *   ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
  75         IF (KODE.EQ.1) GO TO 40
  76         STR = ZR + ZETA2R
  77         STI = ZI + ZETA2I
  78         RAST = FN/AZABS(STR,STI)
  79         STR = STR*RAST*RAST
  80         STI = -STI*RAST*RAST
  81         S1R = -ZETA1R + STR
  82         S1I = -ZETA1I + STI + ZI
  83         GO TO 50
  84    40   CONTINUE
  85         S1R = -ZETA1R + ZETA2R
  86         S1I = -ZETA1I + ZETA2I
  87    50   CONTINUE
  88 C-----------------------------------------------------------------------
  89 C     TEST FOR UNDERFLOW AND OVERFLOW
  90 C-----------------------------------------------------------------------
  91         RS1 = S1R
  92         IF (DABS(RS1).GT.ELIM) GO TO 110
  93         IF (I.EQ.1) IFLAG = 2
  94         IF (DABS(RS1).LT.ALIM) GO TO 60
  95 C-----------------------------------------------------------------------
  96 C     REFINE  TEST AND SCALE
  97 C-----------------------------------------------------------------------
  98         APHI = AZABS(PHIR,PHII)
  99         RS1 = RS1 + DLOG(APHI)
 100         IF (DABS(RS1).GT.ELIM) GO TO 110
 101         IF (I.EQ.1) IFLAG = 1
 102         IF (RS1.LT.0.0D0) GO TO 60
 103         IF (I.EQ.1) IFLAG = 3
 104    60   CONTINUE
 105 C-----------------------------------------------------------------------
 106 C     SCALE S1 IF CABS(S1).LT.ASCLE
 107 C-----------------------------------------------------------------------
 108         S2R = PHIR*SUMR - PHII*SUMI
 109         S2I = PHIR*SUMI + PHII*SUMR
 110         STR = DEXP(S1R)*CSSR(IFLAG)
 111         S1R = STR*DCOS(S1I)
 112         S1I = STR*DSIN(S1I)
 113         STR = S2R*S1R - S2I*S1I
 114         S2I = S2R*S1I + S2I*S1R
 115         S2R = STR
 116         IF (IFLAG.NE.1) GO TO 70
 117         CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 118         IF (NW.NE.0) GO TO 110
 119    70   CONTINUE
 120         CYR(I) = S2R
 121         CYI(I) = S2I
 122         M = ND - I + 1
 123         YR(M) = S2R*CSRR(IFLAG)
 124         YI(M) = S2I*CSRR(IFLAG)
 125    80 CONTINUE
 126       IF (ND.LE.2) GO TO 100
 127       RAST = 1.0D0/AZABS(ZR,ZI)
 128       STR = ZR*RAST
 129       STI = -ZI*RAST
 130       RZR = (STR+STR)*RAST
 131       RZI = (STI+STI)*RAST
 132       BRY(2) = 1.0D0/BRY(1)
 133       BRY(3) = D1MACH(2)
 134       S1R = CYR(1)
 135       S1I = CYI(1)
 136       S2R = CYR(2)
 137       S2I = CYI(2)
 138       C1R = CSRR(IFLAG)
 139       ASCLE = BRY(IFLAG)
 140       K = ND - 2
 141       FN = DBLE(FLOAT(K))
 142       DO 90 I=3,ND
 143         C2R = S2R
 144         C2I = S2I
 145         S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
 146         S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
 147         S1R = C2R
 148         S1I = C2I
 149         C2R = S2R*C1R
 150         C2I = S2I*C1R
 151         YR(K) = C2R
 152         YI(K) = C2I
 153         K = K - 1
 154         FN = FN - 1.0D0
 155         IF (IFLAG.GE.3) GO TO 90
 156         STR = DABS(C2R)
 157         STI = DABS(C2I)
 158         C2M = DMAX1(STR,STI)
 159         IF (C2M.LE.ASCLE) GO TO 90
 160         IFLAG = IFLAG + 1
 161         ASCLE = BRY(IFLAG)
 162         S1R = S1R*C1R
 163         S1I = S1I*C1R
 164         S2R = C2R
 165         S2I = C2I
 166         S1R = S1R*CSSR(IFLAG)
 167         S1I = S1I*CSSR(IFLAG)
 168         S2R = S2R*CSSR(IFLAG)
 169         S2I = S2I*CSSR(IFLAG)
 170         C1R = CSRR(IFLAG)
 171    90 CONTINUE
 172   100 CONTINUE
 173       RETURN
 174 C-----------------------------------------------------------------------
 175 C     SET UNDERFLOW AND UPDATE PARAMETERS
 176 C-----------------------------------------------------------------------
 177   110 CONTINUE
 178       IF (RS1.GT.0.0D0) GO TO 120
 179       YR(ND) = ZEROR
 180       YI(ND) = ZEROI
 181       NZ = NZ + 1
 182       ND = ND - 1
 183       IF (ND.EQ.0) GO TO 100
 184       CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
 185       IF (NUF.LT.0) GO TO 120
 186       ND = ND - NUF
 187       NZ = NZ + NUF
 188       IF (ND.EQ.0) GO TO 100
 189       FN = FNU + DBLE(FLOAT(ND-1))
 190       IF (FN.GE.FNUL) GO TO 30
 191       NLAST = ND
 192       RETURN
 193   120 CONTINUE
 194       NZ = -1
 195       RETURN
 196   130 CONTINUE
 197       IF (RS1.GT.0.0D0) GO TO 120
 198       NZ = N
 199       DO 140 I=1,N
 200         YR(I) = ZEROR
 201         YI(I) = ZEROI
 202   140 CONTINUE
 203       RETURN
 204       END