zbesi.f

Go to the documentation of this file.
00001       SUBROUTINE ZBESI(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
00002 C***BEGIN PROLOGUE  ZBESI
00003 C***DATE WRITTEN   830501   (YYMMDD)
00004 C***REVISION DATE  890801   (YYMMDD)
00005 C***CATEGORY NO.  B5K
00006 C***KEYWORDS  I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
00007 C             MODIFIED BESSEL FUNCTION OF THE FIRST KIND
00008 C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
00009 C***PURPOSE  TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
00010 C***DESCRIPTION
00011 C
00012 C                    ***A DOUBLE PRECISION ROUTINE***
00013 C         ON KODE=1, ZBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
00014 C         BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE
00015 C         ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE
00016 C         -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESI RETURNS THE SCALED
00017 C         FUNCTIONS
00018 C
00019 C         CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z)   J = 1,...,N , X=REAL(Z)
00020 C
00021 C         WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND
00022 C         RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
00023 C         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
00024 C         (REF. 1).
00025 C
00026 C         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
00027 C           ZR,ZI  - Z=CMPLX(ZR,ZI),  -PI.LT.ARG(Z).LE.PI
00028 C           FNU    - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0
00029 C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
00030 C                    KODE= 1  RETURNS
00031 C                             CY(J)=I(FNU+J-1,Z), J=1,...,N
00032 C                        = 2  RETURNS
00033 C                             CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N
00034 C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
00035 C
00036 C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
00037 C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
00038 C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
00039 C                    CY(J)=I(FNU+J-1,Z)  OR
00040 C                    CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X))  J=1,...,N
00041 C                    DEPENDING ON KODE, X=REAL(Z)
00042 C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
00043 C                    NZ= 0   , NORMAL RETURN
00044 C                    NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
00045 C                              TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0)
00046 C                              J = N-NZ+1,...,N
00047 C           IERR   - ERROR FLAG
00048 C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
00049 C                    IERR=1, INPUT ERROR   - NO COMPUTATION
00050 C                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(Z) TOO
00051 C                            LARGE ON KODE=1
00052 C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
00053 C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
00054 C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
00055 C                            ACCURACY
00056 C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
00057 C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
00058 C                            CANCE BY ARGUMENT REDUCTION
00059 C                    IERR=5, ERROR              - NO COMPUTATION,
00060 C                            ALGORITHM TERMINATION CONDITION NOT MET
00061 C
00062 C***LONG DESCRIPTION
00063 C
00064 C         THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR
00065 C         SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z),
00066 C         THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A
00067 C         NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE
00068 C         UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z)
00069 C         FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE
00070 C         SEQUENCES OR REDUCE ORDERS WHEN NECESSARY.
00071 C
00072 C         THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND
00073 C         CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA
00074 C
00075 C         I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z)  REAL(Z).GT.0.0
00076 C                       M = +I OR -I,  I**2=-1
00077 C
00078 C         FOR NEGATIVE ORDERS,THE FORMULA
00079 C
00080 C              I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
00081 C
00082 C         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
00083 C         THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
00084 C         INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE
00085 C         NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
00086 C         K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
00087 C         TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
00088 C         UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
00089 C         OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
00090 C         LARGE MEANS FNU.GT.CABS(Z).
00091 C
00092 C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
00093 C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
00094 C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
00095 C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
00096 C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
00097 C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
00098 C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
00099 C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
00100 C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
00101 C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
00102 C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
00103 C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
00104 C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
00105 C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
00106 C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
00107 C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
00108 C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
00109 C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
00110 C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
00111 C
00112 C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
00113 C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
00114 C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
00115 C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
00116 C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
00117 C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
00118 C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
00119 C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
00120 C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
00121 C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
00122 C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
00123 C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
00124 C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
00125 C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
00126 C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
00127 C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
00128 C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
00129 C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
00130 C         OR -PI/2+P.
00131 C
00132 C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
00133 C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
00134 C                 COMMERCE, 1955.
00135 C
00136 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
00137 C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
00138 C
00139 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
00140 C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
00141 C
00142 C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
00143 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
00144 C                 1018, MAY, 1985
00145 C
00146 C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
00147 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
00148 C                 MATH. SOFTWARE, 1986
00149 C
00150 C***ROUTINES CALLED  ZBINU,I1MACH,D1MACH
00151 C***END PROLOGUE  ZBESI
00152 C     COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN
00153       DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI,
00154      * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR,
00155      * ZR, D1MACH, AZ, BB, FN, XZABS, ASCLE, RTOL, ATOL, STI
00156       INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH
00157       DIMENSION CYR(N), CYI(N)
00158       DATA PI /3.14159265358979324D0/
00159       DATA CONER, CONEI /1.0D0,0.0D0/
00160 C
00161 C***FIRST EXECUTABLE STATEMENT  ZBESI
00162       IERR = 0
00163       NZ=0
00164       IF (FNU.LT.0.0D0) IERR=1
00165       IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
00166       IF (N.LT.1) IERR=1
00167       IF (IERR.NE.0) RETURN
00168 C-----------------------------------------------------------------------
00169 C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
00170 C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
00171 C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
00172 C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
00173 C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
00174 C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
00175 C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
00176 C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
00177 C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
00178 C-----------------------------------------------------------------------
00179       TOL = DMAX1(D1MACH(4),1.0D-18)
00180       K1 = I1MACH(15)
00181       K2 = I1MACH(16)
00182       R1M5 = D1MACH(5)
00183       K = MIN0(IABS(K1),IABS(K2))
00184       ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
00185       K1 = I1MACH(14) - 1
00186       AA = R1M5*DBLE(FLOAT(K1))
00187       DIG = DMIN1(AA,18.0D0)
00188       AA = AA*2.303D0
00189       ALIM = ELIM + DMAX1(-AA,-41.45D0)
00190       RL = 1.2D0*DIG + 3.0D0
00191       FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
00192 C-----------------------------------------------------------------------------
00193 C     TEST FOR PROPER RANGE
00194 C-----------------------------------------------------------------------
00195       AZ = XZABS(ZR,ZI)
00196       FN = FNU+DBLE(FLOAT(N-1))
00197       AA = 0.5D0/TOL
00198       BB=DBLE(FLOAT(I1MACH(9)))*0.5D0
00199       AA = DMIN1(AA,BB)
00200       IF (AZ.GT.AA) GO TO 260
00201       IF (FN.GT.AA) GO TO 260
00202       AA = DSQRT(AA)
00203       IF (AZ.GT.AA) IERR=3
00204       IF (FN.GT.AA) IERR=3
00205       ZNR = ZR
00206       ZNI = ZI
00207       CSGNR = CONER
00208       CSGNI = CONEI
00209       IF (ZR.GE.0.0D0) GO TO 40
00210       ZNR = -ZR
00211       ZNI = -ZI
00212 C-----------------------------------------------------------------------
00213 C     CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
00214 C     WHEN FNU IS LARGE
00215 C-----------------------------------------------------------------------
00216       INU = INT(SNGL(FNU))
00217       ARG = (FNU-DBLE(FLOAT(INU)))*PI
00218       IF (ZI.LT.0.0D0) ARG = -ARG
00219       CSGNR = DCOS(ARG)
00220       CSGNI = DSIN(ARG)
00221       IF (MOD(INU,2).EQ.0) GO TO 40
00222       CSGNR = -CSGNR
00223       CSGNI = -CSGNI
00224    40 CONTINUE
00225       CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL,
00226      * ELIM, ALIM)
00227       IF (NZ.LT.0) GO TO 120
00228       IF (ZR.GE.0.0D0) RETURN
00229 C-----------------------------------------------------------------------
00230 C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
00231 C-----------------------------------------------------------------------
00232       NN = N - NZ
00233       IF (NN.EQ.0) RETURN
00234       RTOL = 1.0D0/TOL
00235       ASCLE = D1MACH(1)*RTOL*1.0D+3
00236       DO 50 I=1,NN
00237 C       STR = CYR(I)*CSGNR - CYI(I)*CSGNI
00238 C       CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR
00239 C       CYR(I) = STR
00240         AA = CYR(I)
00241         BB = CYI(I)
00242         ATOL = 1.0D0
00243         IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55
00244           AA = AA*RTOL
00245           BB = BB*RTOL
00246           ATOL = TOL
00247    55   CONTINUE
00248         STR = AA*CSGNR - BB*CSGNI
00249         STI = AA*CSGNI + BB*CSGNR
00250         CYR(I) = STR*ATOL
00251         CYI(I) = STI*ATOL
00252         CSGNR = -CSGNR
00253         CSGNI = -CSGNI
00254    50 CONTINUE
00255       RETURN
00256   120 CONTINUE
00257       IF(NZ.EQ.(-2)) GO TO 130
00258       NZ = 0
00259       IERR=2
00260       RETURN
00261   130 CONTINUE
00262       NZ=0
00263       IERR=5
00264       RETURN
00265   260 CONTINUE
00266       NZ=0
00267       IERR=4
00268       RETURN
00269       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines