Package gov.nih.mipav.model.algorithms

Class BesselEP

  • public class BesselEP
extends java.lang.Object

    This module computes Bessel functions of complex arguments and a nonnegative order. This work is the port of FORTRAN code written by Donald E. Amos at the Sandia National Laboratories. This is an installation of his Algorithm 644, Collected Algorithms form ACM. The original FORTRAN code was downloaded from http://www.netlib.org/toms/644. The decision to port this code was prompted by the fact that the MATLAB links to this Amos FORTRAN package. For DoubleDoublechecking this code 2 possibilities exist: 1.) The 6 self tests included in the package - zqcbh, zqcbi, zqcbj, zqcbk, zqcby, and zqcai. The code has passed all 6 self tests. 2.) Bessel function tables for real arguments as given in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables Edited by Milton Abramowitz and Irene A. Stegun. Some values for Bessel functions of complex arguments can be found in Computation of Special Functions by Shanjie Zhang and Jianming Jin, copyright 1966 by John Wiley & Sons, Inc. The I and K Bessel functions obtained in this code were found to match the values in these 2 sources. Note that the Zhang and Jin book also contains code for Bessel functions of complex arguments.

    REFERENCES:

    1. Abramowitz, M. and Stegun, I. A., Handbook of Mathematical Functions, NBS Applied Math Series 55, U.S. Dept. of Commerce, Washington, D.C., 1955
    2. Amos, D. E., Algorithm 644, A Portable Package For Bessel Functions of a Complex Argument and Nonnegative Order, ACM Transactions on Mathematical Software, Vol. 12, No. 3, September 1986, Pages 265-273
    3. Amos, D. E., Remark on Algorithm 644, ACM Transactions on Mathematical Software, Vol. 16, No. 4, December 1990, Page 404
    4. 4. Amos, D. E., Remark on Algorithm 644 (Improvements in Algorithm 644), ACM Transactions on Mathematical Software, Vol. 21, No. 4, December 1995, Pages 388-393.
    5. 5. Cody, W. J., Algorithm 665, MACHAR: A Subroutine to Dynamically Determine Machine Parameters, ACM Transactions on Mathematical Software, Vol. 14, No. 4, December 1988, Pages 303-311

    From the ACM website: 

     ACM Software License Agreement
 
  All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
  Noncommercial Use
 
  The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
 
  1. User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
  2. ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
  3. User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
  4. User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
  5. User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
  6. This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
  7. This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
 
  Commercial Use
 
  Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
 
  Revised 6/98

    • Field Detail

      • BESSEL_H

        public static final int BESSEL_H
        Bessel functions of the third kind, Hankel functions.
        See Also:
        Constant Field Values

      • BESSEL_I

        public static final int BESSEL_I
        modified Bessel function of the first kind.
        See Also:
        Constant Field Values

      • BESSEL_J

        public static final int BESSEL_J
        Bessel function of the first kind.
        See Also:
        Constant Field Values

      • BESSEL_K

        public static final int BESSEL_K
        modified Bessel function of the second kind.
        See Also:
        Constant Field Values

      • BESSEL_Y

        public static final int BESSEL_Y
        Bessel function of the second kind.
        See Also:
        Constant Field Values

      • UNSCALED_FUNCTION

        public static final int UNSCALED_FUNCTION
        DOCUMENT ME!
        See Also:
        Constant Field Values

      • SCALED_FUNCTION

        public static final int SCALED_FUNCTION
        multiplied by exp(-abs(Real(z))).
        See Also:
        Constant Field Values

      • BesselType

        private final int BesselType
        DOCUMENT ME!

      • derivativeOrder

        private int derivativeOrder
        order of derivative, 0 or 1.

      • doTest

        private boolean doTest
        DOCUMENT ME!

      • DoubleDoubleDigits

        private int DoubleDoubleDigits
        I1MACH(14) = 106.

      • emax

        private int emax
        I1MACH(16) = 1024 The largest exponent for DoubleDouble precision.

      • emin

        private int emin
        I1MACH(15) = -1021; The smallest exponent for DoubleDouble precision.

      • epsilon

        private DoubleDouble epsilon
        The smallest representable relative difference between two {link @ DoubleDouble} values

      • errorFlag

        private int[] errorFlag
        errorFlag = 0, normal return - computation completed errorFlag = 1, input error - no computation errorFlag = 2, overflow - no computation, Real(z) large on UNSCALED_FUNCTION errorFlag = 3, ABS(Z) or initialOrder + sequenceNumber - 1 large, computation but losses of significance by arguemnt reduction produce less than half of machine accuracy errorFlag = 4, ABS(Z) or initialOrder + sequenceNumber - 1 too large, No computation because of complete losses of significance by argument reduction errorFlag = 5, Error, no computation, argument termination condition not met Length of array passed in must be 1.

      • HankelKind

        private int HankelKind
        real and imaginary outputs Length of arrays passed in must be sequenceNumber Kind of Kankel function = 1 or 2.

      • initialOrder

        private double initialOrder
        order of initial function.

      • nz

        private int[] nz
        Number of components set to zero due to underflow nz[0] = 0, normal return nz[0] > 0, Last nz components of cy set to zero due to underflow, cy[j] = 0.0 + i*0.0 for j = sequenceNumber - nz, ..., sequenceNumber-1 Length of array passed in must be 1.

      • overflowTest

        private boolean overflowTest
        DOCUMENT ME!

      • r1m5

        private DoubleDouble r1m5
        D1MACH(5) = log10(2). = log10(e) * loge(2.0)

      • scalingOption

        private int scalingOption
        UNSCALED_FUNCTION returns cy[j-1] = besselFunction(initialOrder+j-1,z), j = 1,...,sequenceNumber SCALED_FUNCTION returns cy[j-1] = besselFunction(initialOrder+j-1,z)*exp(-abs(Real(z))), j = 1,...,sequenceNumber.

      • sequenceNumber

        private int sequenceNumber
        number of members of the sequence.

      • zi

        private double zi
        DOCUMENT ME!

      • zr

        private double zr
        zr and zi are the real and imaginary parts of the complex argument.

    • Constructor Detail

      • BesselEP

        public BesselEP​(int BesselType,
                boolean overflowTest)
        Creates a new BesselEP object.
        Parameters:
        BesselType - DOCUMENT ME!
        overflowTest - DOCUMENT ME!

      • BesselEP

        public BesselEP​(int BesselType,
                double zr,
                double zi,
                int derivativeOrder,
                int scalingOption,
                DoubleDouble[] cyr,
                DoubleDouble[] cyi,
                int[] nz,
                int[] errorFlag)
        This constructor is used for Ai and Bi Airy functions.
        Parameters:
        BesselType - int
        zr - double
        zi - double
        derivativeOrder - int
        scalingOption - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        errorFlag - int[]

      • BesselEP

        public BesselEP​(int BesselType,
                double zr,
                double zi,
                double initialOrder,
                int scalingOption,
                int sequenceNumber,
                DoubleDouble[] cyr,
                DoubleDouble[] cyi,
                int[] nz,
                int[] errorFlag)
        This constructor used for I, J, K, and Y Bessel functions.
        Parameters:
        BesselType - int
        zr - double
        zi - double
        initialOrder - double
        scalingOption - int
        sequenceNumber - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        errorFlag - int[]

      • BesselEP

        public BesselEP​(int BesselType,
                double zr,
                double zi,
                double initialOrder,
                int scalingOption,
                int HankelKind,
                int sequenceNumber,
                DoubleDouble[] cyr,
                DoubleDouble[] cyi,
                int[] nz,
                int[] errorFlag)
        This constructor used for H functions of kind 1 and 2.
        Parameters:
        BesselType - int
        zr - double
        zi - double
        initialOrder - double
        scalingOption - int
        HankelKind - int
        sequenceNumber - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        errorFlag - int[]

    • Method Detail

      • finalize

        public void finalize()
              throws java.lang.Throwable
        Cleanup memory.
        Overrides:
        finalize in class java.lang.Object
        Throws:
        java.lang.Throwable - DOCUMENT ME!

      • run

        public void run()
        DOCUMENT ME!

      • dgamln

        private DoubleDouble dgamln​(DoubleDouble z,
                            int[] ierr)
        DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN. THE FUNCTION WAS MADE AS PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE 10 DIGITS IN A WORD, RLN=AMAX1(-ALOG10(R1MACH(4)),0.5E-18) LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.

        SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100 VALUES IS USED FOR SPEED OF EXECUTION.

        Parameters:
        z - DoubleDouble, z > 0.0
        ierr - int[] error flag IERR=0, NORMAL RETURN, COMPUTATION COMPLETED IERR=1, Z.LE.0.0D0, NO COMPUTATION
        Returns:
        DoubleDouble NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0

      • runTest

        private void runTest()
        DOCUMENT ME!

      • zabs

        private DoubleDouble zabs​(DoubleDouble zr,
                          DoubleDouble zi)
        zabs computes the absolute value or magnitude of a DoubleDouble precision complex variable zr + j*zi.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        Returns:
        DoubleDouble

      • zacai

        private void zacai​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int mr,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA.

        K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) MP=PI*MR*CMPLX(0.0,1.0)

        TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1. ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON IS CALLED FROM ZAIRY.

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        mr - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zacon

        private void zacon​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int mr,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA.

        K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) MP=PI*MR*CMPLX(0.0,1.0)

        TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT HALF Z PLANE

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        mr - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zairy

        private void zairy​(DoubleDouble zr,
                   DoubleDouble zi,
                   int id,
                   int kode,
                   DoubleDouble[] air,
                   DoubleDouble[] aii,
                   int[] nz,
                   int[] ierr)
        ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)* DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z).

        WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS. DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF C MATHEMATICAL FUNCTIONS (REF. 1). C C INPUT ZR,ZI ARE DOUBLE PRECISION C ZR,ZI - Z=CMPLX(ZR,ZI) C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C AI=AI(Z) ON ID=0 OR C AI=DAI(Z)/DZ ON ID=1 C = 2 RETURNS C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE C ZTA=(2/3)*Z*CSQRT(Z) C C OUTPUT AIR,AII ARE DOUBLE PRECISION C AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND C KODE C NZ - UNDERFLOW INDICATOR C NZ= 0 , NORMAL RETURN C NZ= 1 , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN C -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1 C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA) C TOO LARGE ON KODE=1 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION C PRODUCE LESS THAN HALF OF MACHINE ACCURACY C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION C COMPLETE LOSS OF ACCURACY BY ARGUMENT C REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL C FUNCTIONS BY C C AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA) C C=1.0/(PI*SQRT(3.0)) C ZTA=(2/3)*Z**(3/2) C C WITH THE POWER SERIES FOR CABS(Z).LE.1.0. C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR C FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER C MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, ACM C TRANS. MATH. SOFTWARE, VOL. 12, NO. 3, SEPTEMBER 1986, C PP 265-273.

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        id - int
        kode - int
        air - DoubleDouble[]
        aii - DoubleDouble[]
        nz - int[]
        ierr - int[]

      • zasyi

        private void zasyi​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE REGION CABS(Z).GT.MAX(RL,fnu*fnu/2). NZ=0 IS A NORMAL RETURN. NZ.LT.0 INDICATES AN OVERFLOW ON kode=UNSCALED_FUNCTION.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zbesh

        private void zbesh​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int m,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        m - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        ierr - int[]

      • zbesi

        private void zbesi​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        ierr - int[]

      • zbesj

        private void zbesj​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        ierr - int[]

      • zbesk

        private void zbesk​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        ierr - int[]

      • zbesy

        private void zbesy​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz,
                   DoubleDouble[] cwrkr,
                   DoubleDouble[] cwrki,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        cwrkr - DoubleDouble[]
        cwrki - DoubleDouble[]
        ierr - int[]

      • zbesyh

        private void zbesyh​(DoubleDouble zr,
                    DoubleDouble zi,
                    DoubleDouble fnu,
                    int kode,
                    int n,
                    DoubleDouble[] cyr,
                    DoubleDouble[] cyi,
                    int[] nz,
                    DoubleDouble[] cwrkr,
                    DoubleDouble[] cwrki,
                    int[] ierr)
        An older version of zbesy used by the ZQCBY quick check routine.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]
        cwrkr - DoubleDouble[]
        cwrki - DoubleDouble[]
        ierr - int[]

      • zbinu

        private void zbinu​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi,
                   int[] nz)
        zbinu calculates the I function in the right half z plane.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]
        nz - int[]

      • zbiry

        private void zbiry​(DoubleDouble zr,
                   DoubleDouble zi,
                   int id,
                   int kode,
                   DoubleDouble[] bir,
                   DoubleDouble[] bii,
                   int[] ierr)
        DOCUMENT ME!
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        id - int
        kode - int
        bir - DoubleDouble[]
        bii - DoubleDouble[]
        ierr - int[]

      • zbknu

        private void zbknu​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        zbknu computes the k Bessel function in the right half z plane.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zbuni

        private void zbuni​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz,
                   int nui,
                   int[] nlast,
                   DoubleDouble fnul)
        ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT. FNUL AND fnu+N-1.LT.FNUL. THE ORDER IS INCREASED FROM fnu+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(fnu,Z) ON IFORM=1 AND THE EXPANSION FOR J(fnu,Z) ON IFORM=2
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]
        nui - int
        nlast - int[]
        fnul - DoubleDouble

      • zbunk

        private void zbunk​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int mr,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU.GT.FNUL. ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR K(FNU,Z) IN ZUNK1 AND THE EXPANSION FOR H(2,FNU,Z) IN ZUNK2
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        mr - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zkscl

        private void zkscl​(DoubleDouble zrr,
                   DoubleDouble zri,
                   DoubleDouble fnu,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz,
                   DoubleDouble rzr,
                   DoubleDouble rzi,
                   DoubleDouble ascle)
        Set K functions to zero on underflow, continue recurrence on scaled functions until two members come on scale, then return with min(zn[0]+2,n) values scaled by 1/tol.
        Parameters:
        zrr - DoubleDouble
        zri - DoubleDouble
        fnu - DoubleDouble
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]
        rzr - DoubleDouble
        rzi - DoubleDouble
        ascle - DoubleDouble

      • zlog

        private void zlog​(DoubleDouble ar,
                  DoubleDouble ai,
                  DoubleDouble[] br,
                  DoubleDouble[] bi,
                  int[] ierr)
        complex logarithm b = clog(a).
        Parameters:
        ar - DoubleDouble
        ai - DoubleDouble
        br - DoubleDouble[]
        bi - DoubleDouble[]
        ierr - int[] ierr = 0, normal return ierr = 1, z = cmplx(0.0, 0.0)

      • zmlri

        private void zmlri​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        zmlri computes the I Bessel function for the Real(z) >= 0.0 by the Miller Algorithm normalized by a Neumann series.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zqcai

        private void zqcai​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zqcbh

        private void zqcbh​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zqcbi

        private void zqcbi​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zqcbj

        private void zqcbj​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zqcbk

        private void zqcbk​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zqcby

        private void zqcby​(int mqc)
        DOCUMENT ME!
        Parameters:
        mqc - int Default is mqc = 1 does not check underflow limits and overflow limits. mqc != 1 checks underflow and overflow limits

      • zrati

        private void zrati​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int n,
                   DoubleDouble[] cyr,
                   DoubleDouble[] cyi)
        ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, BY D. J. SOOKNE.
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        n - int
        cyr - DoubleDouble[]
        cyi - DoubleDouble[]

      • zs1s2

        private void zs1s2​(DoubleDouble zrr,
                   DoubleDouble zri,
                   DoubleDouble[] s1r,
                   DoubleDouble[] s1i,
                   DoubleDouble[] s2r,
                   DoubleDouble[] s2i,
                   int[] nz,
                   DoubleDouble ascle,
                   int[] iuf)
        ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON- TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION. ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE PRECISION ABOVE THE UNDERFLOW LIMIT.
        Parameters:
        zrr - DoubleDouble
        zri - DoubleDouble
        s1r - DoubleDouble[]
        s1i - DoubleDouble[]
        s2r - DoubleDouble[]
        s2i - DoubleDouble[]
        nz - int[]
        ascle - DoubleDouble
        iuf - int[]

      • zseri

        private void zseri​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE REGION CABS(Z).LE.2*SQRT(fnu+1). NZ=0 IS A NORMAL RETURN. NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE CONDITION CABS(Z).LE.2*SQRT(fnu+1) WAS VIOLATED AND THE COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zshch

        private void zshch​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble[] cshr,
                   DoubleDouble[] cshi,
                   DoubleDouble[] cchr,
                   DoubleDouble[] cchi)
        zshch computes the complex hyperbolic functions csh = sinh(x+i*y) and cch = cosh(x+i*y).
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        cshr - DoubleDouble[]
        cshi - DoubleDouble[]
        cchr - DoubleDouble[]
        cchi - DoubleDouble[]

      • zuchk

        private void zuchk​(DoubleDouble yr,
                   DoubleDouble yi,
                   int[] nz,
                   DoubleDouble ascle)
        Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN EXP(-ALIM)=ASCLE=1.0E+3*D1MACH(1)/TOL. THE TEST IS MADE TO SEE IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED.
        Parameters:
        yr - DoubleDouble
        yi - DoubleDouble
        nz - int[]
        ascle - DoubleDouble

      • zunhj

        private void zunhj​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int ipmtr,
                   DoubleDouble[] phir,
                   DoubleDouble[] phii,
                   DoubleDouble[] argr,
                   DoubleDouble[] argi,
                   DoubleDouble[] zeta1r,
                   DoubleDouble[] zeta1i,
                   DoubleDouble[] zeta2r,
                   DoubleDouble[] zeta2i,
                   DoubleDouble[] asumr,
                   DoubleDouble[] asumi,
                   DoubleDouble[] bsumr,
                   DoubleDouble[] bsumi)
        REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9.

        ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC PRESS, N.Y., 1974, PAGE 420

        ABSTRACT ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION

        C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) )

        FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE.

        (2/3)*FNU*ZETA**1.5 = ZETA1-ZETA2,

        ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY.

        MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR=1 COMPUTES ALL EXCEPT ASUM AND BSUM.

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        ipmtr - int
        phir - DoubleDouble[]
        phii - DoubleDouble[]
        argr - DoubleDouble[]
        argi - DoubleDouble[]
        zeta1r - DoubleDouble[]
        zeta1i - DoubleDouble[]
        zeta2r - DoubleDouble[]
        zeta2i - DoubleDouble[]
        asumr - DoubleDouble[]
        asumi - DoubleDouble[]
        bsumr - DoubleDouble[]
        bsumi - DoubleDouble[]

      • zuni1

        private void zuni1​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz,
                   int[] nlast,
                   DoubleDouble fnul)
        ZUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.

        NUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. Y(I)=CZERO FOR I=NLAST+1,N

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]
        nlast - int[]
        fnul - DoubleDouble

      • zuni2

        private void zuni2​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz,
                   int[] nlast,
                   DoubleDouble fnul)
        ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.

        FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. Y(I)=CZERO FOR I=NLAST+1,N

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]
        nlast - int[]
        fnul - DoubleDouble

      • zunik

        private void zunik​(DoubleDouble zrr,
                   DoubleDouble zri,
                   DoubleDouble fnu,
                   int ikflg,
                   int ipmtr,
                   int[] init,
                   DoubleDouble[] phir,
                   DoubleDouble[] phii,
                   DoubleDouble[] zeta1r,
                   DoubleDouble[] zeta1i,
                   DoubleDouble[] zeta2r,
                   DoubleDouble[] zeta2i,
                   DoubleDouble[] sumr,
                   DoubleDouble[] sumi,
                   DoubleDouble[] cwrkr,
                   DoubleDouble[] cwrki)
        ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 RESPECTIVELY BY.

        W(fnu,ZR) = PHI*EXP(ZETA)*SUM

        WHERE ZETA=-ZETA1 + ZETA2 OR ZETA1 - ZETA2

        THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, ZETA1,ZETA2.

        Parameters:
        zrr - DoubleDouble
        zri - DoubleDouble
        fnu - DoubleDouble
        ikflg - int
        ipmtr - int
        init - int
        phir - DoubleDouble[]
        phii - DoubleDouble[]
        zeta1r - DoubleDouble[]
        zeta1i - DoubleDouble[]
        zeta2r - DoubleDouble[]
        zeta2i - DoubleDouble[]
        sumr - DoubleDouble[]
        sumi - DoubleDouble[]
        cwrkr - DoubleDouble[]
        cwrki - DoubleDouble[]

      • zunk1

        private void zunk1​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int mr,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION. MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. NZ=-1 MEANS AN OVERFLOW WILL OCCUR
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        mr - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zunk2

        private void zunk2​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int mr,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz)
        ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN) WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC- ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. NZ=-1 MEANS AN OVERFLOW WILL OCCUR
        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        mr - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]

      • zuoik

        private void zuoik​(DoubleDouble zr,
                   DoubleDouble zi,
                   DoubleDouble fnu,
                   int kode,
                   int ikflg,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nuf)
        ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)= EXP(-ELIM)/TOL

        IKFLG=1 MEANS THE I SEQUENCE IS TESTED =2 MEANS THE K SEQUENCE IS TESTED NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE =-1 MEANS AN OVERFLOW WOULD OCCUR IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY ANOTHER ROUTINE

        Parameters:
        zr - DoubleDouble
        zi - DoubleDouble
        fnu - DoubleDouble
        kode - int
        ikflg - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nuf - int[]

      • zwrsk

        private void zwrsk​(DoubleDouble zrr,
                   DoubleDouble zri,
                   DoubleDouble fnu,
                   int kode,
                   int n,
                   DoubleDouble[] yr,
                   DoubleDouble[] yi,
                   int[] nz,
                   DoubleDouble[] cwr,
                   DoubleDouble[] cwi)
        zwrsk computes the I Bessel function for Real(z) >- 0.0 by normalizing the I function ratios from zrati by the Wronskian.
        Parameters:
        zrr - DoubleDouble
        zri - DoubleDouble
        fnu - DoubleDouble
        kode - int
        n - int
        yr - DoubleDouble[]
        yi - DoubleDouble[]
        nz - int[]
        cwr - DoubleDouble[]
        cwi - DoubleDouble[]