amos/zseri.f

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