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

Go to the documentation of this file.
00001       SUBROUTINE CBUNI(Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL,
00002      * ELIM, ALIM)
00003 C***BEGIN PROLOGUE  CBUNI
00004 C***REFER TO  CBESI,CBESK
00005 C
00006 C     CBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
00007 C     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
00008 C     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
00009 C     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
00010 C     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
00011 C
00012 C***ROUTINES CALLED  CUNI1,CUNI2,R1MACH
00013 C***END PROLOGUE  CBUNI
00014       COMPLEX CSCL, CSCR, CY, RZ, ST, S1, S2, Y, Z
00015       REAL ALIM, AX, AY, DFNU, ELIM, FNU, FNUI, FNUL, GNU, TOL, XX, YY,
00016      * ASCLE, BRY, STR, STI, STM, R1MACH
00017       INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
00018       DIMENSION Y(N), CY(2), BRY(3)
00019       NZ = 0
00020       XX = REAL(Z)
00021       YY = AIMAG(Z)
00022       AX = ABS(XX)*1.7321E0
00023       AY = ABS(YY)
00024       IFORM = 1
00025       IF (AY.GT.AX) IFORM = 2
00026       IF (NUI.EQ.0) GO TO 60
00027       FNUI = FLOAT(NUI)
00028       DFNU = FNU + FLOAT(N-1)
00029       GNU = DFNU + FNUI
00030       IF (IFORM.EQ.2) GO TO 10
00031 C-----------------------------------------------------------------------
00032 C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
00033 C     -PI/3.LE.ARG(Z).LE.PI/3
00034 C-----------------------------------------------------------------------
00035       CALL CUNI1(Z, GNU, KODE, 2, CY, NW, NLAST, FNUL, TOL, ELIM, ALIM)
00036       GO TO 20
00037    10 CONTINUE
00038 C-----------------------------------------------------------------------
00039 C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
00040 C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
00041 C     AND HPI=PI/2
00042 C-----------------------------------------------------------------------
00043       CALL CUNI2(Z, GNU, KODE, 2, CY, NW, NLAST, FNUL, TOL, ELIM, ALIM)
00044    20 CONTINUE
00045       IF (NW.LT.0) GO TO 50
00046       IF (NW.NE.0) GO TO 90
00047       AY = CABS(CY(1))
00048 C----------------------------------------------------------------------
00049 C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
00050 C----------------------------------------------------------------------
00051       BRY(1) = 1.0E+3*R1MACH(1)/TOL
00052       BRY(2) = 1.0E0/BRY(1)
00053       BRY(3) = BRY(2)
00054       IFLAG = 2
00055       ASCLE = BRY(2)
00056       AX = 1.0E0
00057       CSCL = CMPLX(AX,0.0E0)
00058       IF (AY.GT.BRY(1)) GO TO 21
00059       IFLAG = 1
00060       ASCLE = BRY(1)
00061       AX = 1.0E0/TOL
00062       CSCL = CMPLX(AX,0.0E0)
00063       GO TO 25
00064    21 CONTINUE
00065       IF (AY.LT.BRY(2)) GO TO 25
00066       IFLAG = 3
00067       ASCLE = BRY(3)
00068       AX = TOL
00069       CSCL = CMPLX(AX,0.0E0)
00070    25 CONTINUE
00071       AY = 1.0E0/AX
00072       CSCR = CMPLX(AY,0.0E0)
00073       S1 = CY(2)*CSCL
00074       S2 = CY(1)*CSCL
00075       RZ = CMPLX(2.0E0,0.0E0)/Z
00076       DO 30 I=1,NUI
00077         ST = S2
00078         S2 = CMPLX(DFNU+FNUI,0.0E0)*RZ*S2 + S1
00079         S1 = ST
00080         FNUI = FNUI - 1.0E0
00081         IF (IFLAG.GE.3) GO TO 30
00082         ST = S2*CSCR
00083         STR = REAL(ST)
00084         STI = AIMAG(ST)
00085         STR = ABS(STR)
00086         STI = ABS(STI)
00087         STM = AMAX1(STR,STI)
00088         IF (STM.LE.ASCLE) GO TO 30
00089         IFLAG = IFLAG+1
00090         ASCLE = BRY(IFLAG)
00091         S1 = S1*CSCR
00092         S2 = ST
00093         AX = AX*TOL
00094         AY = 1.0E0/AX
00095         CSCL = CMPLX(AX,0.0E0)
00096         CSCR = CMPLX(AY,0.0E0)
00097         S1 = S1*CSCL
00098         S2 = S2*CSCL
00099    30 CONTINUE
00100       Y(N) = S2*CSCR
00101       IF (N.EQ.1) RETURN
00102       NL = N - 1
00103       FNUI = FLOAT(NL)
00104       K = NL
00105       DO 40 I=1,NL
00106         ST = S2
00107         S2 = CMPLX(FNU+FNUI,0.0E0)*RZ*S2 + S1
00108         S1 = ST
00109         ST = S2*CSCR
00110         Y(K) = ST
00111         FNUI = FNUI - 1.0E0
00112         K = K - 1
00113         IF (IFLAG.GE.3) GO TO 40
00114         STR = REAL(ST)
00115         STI = AIMAG(ST)
00116         STR = ABS(STR)
00117         STI = ABS(STI)
00118         STM = AMAX1(STR,STI)
00119         IF (STM.LE.ASCLE) GO TO 40
00120         IFLAG = IFLAG+1
00121         ASCLE = BRY(IFLAG)
00122         S1 = S1*CSCR
00123         S2 = ST
00124         AX = AX*TOL
00125         AY = 1.0E0/AX
00126         CSCL = CMPLX(AX,0.0E0)
00127         CSCR = CMPLX(AY,0.0E0)
00128         S1 = S1*CSCL
00129         S2 = S2*CSCL
00130    40 CONTINUE
00131       RETURN
00132    50 CONTINUE
00133       NZ = -1
00134       IF(NW.EQ.(-2)) NZ=-2
00135       RETURN
00136    60 CONTINUE
00137       IF (IFORM.EQ.2) GO TO 70
00138 C-----------------------------------------------------------------------
00139 C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
00140 C     -PI/3.LE.ARG(Z).LE.PI/3
00141 C-----------------------------------------------------------------------
00142       CALL CUNI1(Z, FNU, KODE, N, Y, NW, NLAST, FNUL, TOL, ELIM, ALIM)
00143       GO TO 80
00144    70 CONTINUE
00145 C-----------------------------------------------------------------------
00146 C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
00147 C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
00148 C     AND HPI=PI/2
00149 C-----------------------------------------------------------------------
00150       CALL CUNI2(Z, FNU, KODE, N, Y, NW, NLAST, FNUL, TOL, ELIM, ALIM)
00151    80 CONTINUE
00152       IF (NW.LT.0) GO TO 50
00153       NZ = NW
00154       RETURN
00155    90 CONTINUE
00156       NLAST = N
00157       RETURN
00158       END