src/numerical/slatec/fortran/zseri.f

   1 *DECK ZSERI
   2       SUBROUTINE ZSERI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
   3      +   ALIM)
   4 C***BEGIN PROLOGUE  ZSERI
   5 C***SUBSIDIARY
   6 C***PURPOSE  Subsidiary to ZBESI and ZBESK
   7 C***LIBRARY   SLATEC
   8 C***TYPE      ALL (CSERI-A, ZSERI-A)
   9 C***AUTHOR  Amos, D. E., (SNL)
  10 C***DESCRIPTION
  11 C
  12 C     ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
  13 C     MEANS OF THE POWER SERIES FOR LARGE ABS(Z) IN THE
  14 C     REGION ABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
  15 C     NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
  16 C     DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
  17 C     CONDITION ABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
  18 C     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
  19 C
  20 C***SEE ALSO  ZBESI, ZBESK
  21 C***ROUTINES CALLED  D1MACH, DGAMLN, ZABS, ZDIV, ZLOG, ZMLT, ZUCHK
  22 C***REVISION HISTORY  (YYMMDD)
  23 C   830501  DATE WRITTEN
  24 C   910415  Prologue converted to Version 4.0 format.  (BAB)
  25 C   930122  Added ZLOG to EXTERNAL statement.  (RWC)
  26 C***END PROLOGUE  ZSERI
  27 C     COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z
  28       DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL,
  29      * AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU,
  30      * ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI,
  31      * STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI,
  32      * ZR, DGAMLN, D1MACH, ZABS
  33       INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW
  34       DIMENSION YR(N), YI(N), WR(2), WI(2)
  35       EXTERNAL ZABS, ZLOG
  36       DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
  37 C***FIRST EXECUTABLE STATEMENT  ZSERI
  38       NZ = 0
  39       AZ = ZABS(ZR,ZI)
  40       IF (AZ.EQ.0.0D0) GO TO 160
  41       ARM = 1.0D+3*D1MACH(1)
  42       RTR1 = SQRT(ARM)
  43       CRSCR = 1.0D0
  44       IFLAG = 0
  45       IF (AZ.LT.ARM) GO TO 150
  46       HZR = 0.5D0*ZR
  47       HZI = 0.5D0*ZI
  48       CZR = ZEROR
  49       CZI = ZEROI
  50       IF (AZ.LE.RTR1) GO TO 10
  51       CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI)
  52    10 CONTINUE
  53       ACZ = ZABS(CZR,CZI)
  54       NN = N
  55       CALL ZLOG(HZR, HZI, CKR, CKI, IDUM)
  56    20 CONTINUE
  57       DFNU = FNU + (NN-1)
  58       FNUP = DFNU + 1.0D0
  59 C-----------------------------------------------------------------------
  60 C     UNDERFLOW TEST
  61 C-----------------------------------------------------------------------
  62       AK1R = CKR*DFNU
  63       AK1I = CKI*DFNU
  64       AK = DGAMLN(FNUP,IDUM)
  65       AK1R = AK1R - AK
  66       IF (KODE.EQ.2) AK1R = AK1R - ZR
  67       IF (AK1R.GT.(-ELIM)) GO TO 40
  68    30 CONTINUE
  69       NZ = NZ + 1
  70       YR(NN) = ZEROR
  71       YI(NN) = ZEROI
  72       IF (ACZ.GT.DFNU) GO TO 190
  73       NN = NN - 1
  74       IF (NN.EQ.0) RETURN
  75       GO TO 20
  76    40 CONTINUE
  77       IF (AK1R.GT.(-ALIM)) GO TO 50
  78       IFLAG = 1
  79       SS = 1.0D0/TOL
  80       CRSCR = TOL
  81       ASCLE = ARM*SS
  82    50 CONTINUE
  83       AA = EXP(AK1R)
  84       IF (IFLAG.EQ.1) AA = AA*SS
  85       COEFR = AA*COS(AK1I)
  86       COEFI = AA*SIN(AK1I)
  87       ATOL = TOL*ACZ/FNUP
  88       IL = MIN(2,NN)
  89       DO 90 I=1,IL
  90         DFNU = FNU + (NN-I)
  91         FNUP = DFNU + 1.0D0
  92         S1R = CONER
  93         S1I = CONEI
  94         IF (ACZ.LT.TOL*FNUP) GO TO 70
  95         AK1R = CONER
  96         AK1I = CONEI
  97         AK = FNUP + 2.0D0
  98         S = FNUP
  99         AA = 2.0D0
 100    60   CONTINUE
 101         RS = 1.0D0/S
 102         STR = AK1R*CZR - AK1I*CZI
 103         STI = AK1R*CZI + AK1I*CZR
 104         AK1R = STR*RS
 105         AK1I = STI*RS
 106         S1R = S1R + AK1R
 107         S1I = S1I + AK1I
 108         S = S + AK
 109         AK = AK + 2.0D0
 110         AA = AA*ACZ*RS
 111         IF (AA.GT.ATOL) GO TO 60
 112    70   CONTINUE
 113         S2R = S1R*COEFR - S1I*COEFI
 114         S2I = S1R*COEFI + S1I*COEFR
 115         WR(I) = S2R
 116         WI(I) = S2I
 117         IF (IFLAG.EQ.0) GO TO 80
 118         CALL ZUCHK(S2R, S2I, NW, ASCLE, TOL)
 119         IF (NW.NE.0) GO TO 30
 120    80   CONTINUE
 121         M = NN - I + 1
 122         YR(M) = S2R*CRSCR
 123         YI(M) = S2I*CRSCR
 124         IF (I.EQ.IL) GO TO 90
 125         CALL ZDIV(COEFR, COEFI, HZR, HZI, STR, STI)
 126         COEFR = STR*DFNU
 127         COEFI = STI*DFNU
 128    90 CONTINUE
 129       IF (NN.LE.2) RETURN
 130       K = NN - 2
 131       AK = K
 132       RAZ = 1.0D0/AZ
 133       STR = ZR*RAZ
 134       STI = -ZI*RAZ
 135       RZR = (STR+STR)*RAZ
 136       RZI = (STI+STI)*RAZ
 137       IF (IFLAG.EQ.1) GO TO 120
 138       IB = 3
 139   100 CONTINUE
 140       DO 110 I=IB,NN
 141         YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
 142         YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
 143         AK = AK - 1.0D0
 144         K = K - 1
 145   110 CONTINUE
 146       RETURN
 147 C-----------------------------------------------------------------------
 148 C     RECUR BACKWARD WITH SCALED VALUES
 149 C-----------------------------------------------------------------------
 150   120 CONTINUE
 151 C-----------------------------------------------------------------------
 152 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
 153 C     UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3
 154 C-----------------------------------------------------------------------
 155       S1R = WR(1)
 156       S1I = WI(1)
 157       S2R = WR(2)
 158       S2I = WI(2)
 159       DO 130 L=3,NN
 160         CKR = S2R
 161         CKI = S2I
 162         S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI)
 163         S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR)
 164         S1R = CKR
 165         S1I = CKI
 166         CKR = S2R*CRSCR
 167         CKI = S2I*CRSCR
 168         YR(K) = CKR
 169         YI(K) = CKI
 170         AK = AK - 1.0D0
 171         K = K - 1
 172         IF (ZABS(CKR,CKI).GT.ASCLE) GO TO 140
 173   130 CONTINUE
 174       RETURN
 175   140 CONTINUE
 176       IB = L + 1
 177       IF (IB.GT.NN) RETURN
 178       GO TO 100
 179   150 CONTINUE
 180       NZ = N
 181       IF (FNU.EQ.0.0D0) NZ = NZ - 1
 182   160 CONTINUE
 183       YR(1) = ZEROR
 184       YI(1) = ZEROI
 185       IF (FNU.NE.0.0D0) GO TO 170
 186       YR(1) = CONER
 187       YI(1) = CONEI
 188   170 CONTINUE
 189       IF (N.EQ.1) RETURN
 190       DO 180 I=2,N
 191         YR(I) = ZEROR
 192         YI(I) = ZEROI
 193   180 CONTINUE
 194       RETURN
 195 C-----------------------------------------------------------------------
 196 C     RETURN WITH NZ.LT.0 IF ABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
 197 C     THE CALCULATION IN CBINU WITH N=N-ABS(NZ)
 198 C-----------------------------------------------------------------------
 199   190 CONTINUE
 200       NZ = -NZ
 201       RETURN
 202       END