GNU Octave: libcruft/amos/cseri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSERI(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
00002 C***BEGIN PROLOGUE  CSERI
00003 C***REFER TO  CBESI,CBESK
00004 C
00005 C     CSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
00006 C     MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
00007 C     REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
00008 C     NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
00009 C     DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
00010 C     CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
00011 C     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
00012 C
00013 C***ROUTINES CALLED  CUCHK,GAMLN,R1MACH
00014 C***END PROLOGUE  CSERI
00015       COMPLEX AK1, CK, COEF, CONE, CRSC, CZ, CZERO, HZ, RZ, S1, S2, W,
00016      * Y, Z
00017       REAL AA, ACZ, AK, ALIM, ARM, ASCLE, ATOL, AZ, DFNU, ELIM, FNU,
00018      * FNUP, RAK1, RS, RTR1, S, SS, TOL, X, GAMLN, R1MACH
00019       INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NW, NZ
00020       DIMENSION Y(N), W(2)
00021       DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
00022 C
00023       NZ = 0
00024       AZ = CABS(Z)
00025       IF (AZ.EQ.0.0E0) GO TO 150
00026       X = REAL(Z)
00027       ARM = 1.0E+3*R1MACH(1)
00028       RTR1 = SQRT(ARM)
00029       CRSC = CMPLX(1.0E0,0.0E0)
00030       IFLAG = 0
00031       IF (AZ.LT.ARM) GO TO 140
00032       HZ = Z*CMPLX(0.5E0,0.0E0)
00033       CZ = CZERO
00034       IF (AZ.GT.RTR1) CZ = HZ*HZ
00035       ACZ = CABS(CZ)
00036       NN = N
00037       CK = CLOG(HZ)
00038    10 CONTINUE
00039       DFNU = FNU + FLOAT(NN-1)
00040       FNUP = DFNU + 1.0E0
00041 C-----------------------------------------------------------------------
00042 C     UNDERFLOW TEST
00043 C-----------------------------------------------------------------------
00044       AK1 = CK*CMPLX(DFNU,0.0E0)
00045       AK = GAMLN(FNUP,IDUM)
00046       AK1 = AK1 - CMPLX(AK,0.0E0)
00047       IF (KODE.EQ.2) AK1 = AK1 - CMPLX(X,0.0E0)
00048       RAK1 = REAL(AK1)
00049       IF (RAK1.GT.(-ELIM)) GO TO 30
00050    20 CONTINUE
00051       NZ = NZ + 1
00052       Y(NN) = CZERO
00053       IF (ACZ.GT.DFNU) GO TO 170
00054       NN = NN - 1
00055       IF (NN.EQ.0) RETURN
00056       GO TO 10
00057    30 CONTINUE
00058       IF (RAK1.GT.(-ALIM)) GO TO 40
00059       IFLAG = 1
00060       SS = 1.0E0/TOL
00061       CRSC = CMPLX(TOL,0.0E0)
00062       ASCLE = ARM*SS
00063    40 CONTINUE
00064       AK = AIMAG(AK1)
00065       AA = EXP(RAK1)
00066       IF (IFLAG.EQ.1) AA = AA*SS
00067       COEF = CMPLX(AA,0.0E0)*CMPLX(COS(AK),SIN(AK))
00068       ATOL = TOL*ACZ/FNUP
00069       IL = MIN0(2,NN)
00070       DO 80 I=1,IL
00071         DFNU = FNU + FLOAT(NN-I)
00072         FNUP = DFNU + 1.0E0
00073         S1 = CONE
00074         IF (ACZ.LT.TOL*FNUP) GO TO 60
00075         AK1 = CONE
00076         AK = FNUP + 2.0E0
00077         S = FNUP
00078         AA = 2.0E0
00079    50   CONTINUE
00080         RS = 1.0E0/S
00081         AK1 = AK1*CZ*CMPLX(RS,0.0E0)
00082         S1 = S1 + AK1
00083         S = S + AK
00084         AK = AK + 2.0E0
00085         AA = AA*ACZ*RS
00086         IF (AA.GT.ATOL) GO TO 50
00087    60   CONTINUE
00088         M = NN - I + 1
00089         S2 = S1*COEF
00090         W(I) = S2
00091         IF (IFLAG.EQ.0) GO TO 70
00092         CALL CUCHK(S2, NW, ASCLE, TOL)
00093         IF (NW.NE.0) GO TO 20
00094    70   CONTINUE
00095         Y(M) = S2*CRSC
00096         IF (I.NE.IL) COEF = COEF*CMPLX(DFNU,0.0E0)/HZ
00097    80 CONTINUE
00098       IF (NN.LE.2) RETURN
00099       K = NN - 2
00100       AK = FLOAT(K)
00101       RZ = (CONE+CONE)/Z
00102       IF (IFLAG.EQ.1) GO TO 110
00103       IB = 3
00104    90 CONTINUE
00105       DO 100 I=IB,NN
00106         Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2)
00107         AK = AK - 1.0E0
00108         K = K - 1
00109   100 CONTINUE
00110       RETURN
00111 C-----------------------------------------------------------------------
00112 C     RECUR BACKWARD WITH SCALED VALUES
00113 C-----------------------------------------------------------------------
00114   110 CONTINUE
00115 C-----------------------------------------------------------------------
00116 C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
00117 C     UNDERFLOW LIMIT = ASCLE = R1MACH(1)*CSCL*1.0E+3
00118 C-----------------------------------------------------------------------
00119       S1 = W(1)
00120       S2 = W(2)
00121       DO 120 L=3,NN
00122         CK = S2
00123         S2 = S1 + CMPLX(AK+FNU,0.0E0)*RZ*S2
00124         S1 = CK
00125         CK = S2*CRSC
00126         Y(K) = CK
00127         AK = AK - 1.0E0
00128         K = K - 1
00129         IF (CABS(CK).GT.ASCLE) GO TO 130
00130   120 CONTINUE
00131       RETURN
00132   130 CONTINUE
00133       IB = L + 1
00134       IF (IB.GT.NN) RETURN
00135       GO TO 90
00136   140 CONTINUE
00137       NZ = N
00138       IF (FNU.EQ.0.0E0) NZ = NZ - 1
00139   150 CONTINUE
00140       Y(1) = CZERO
00141       IF (FNU.EQ.0.0E0) Y(1) = CONE
00142       IF (N.EQ.1) RETURN
00143       DO 160 I=2,N
00144         Y(I) = CZERO
00145   160 CONTINUE
00146       RETURN
00147 C-----------------------------------------------------------------------
00148 C     RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
00149 C     THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
00150 C-----------------------------------------------------------------------
00151   170 CONTINUE
00152       NZ = -NZ
00153       RETURN
00154       END