Package gov.nih.mipav.model.algorithms

Class Bessel

  • public class Bessel
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 doublechecking 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 Summary

      Fields 
      Modifier and Type Field Description
      static int AIRY_AI
      Airy function Ai.
      static int AIRY_BI
      Airy function Bi.
      private double alim
      DOCUMENT ME!
      static int BESSEL_H
      Bessel functions of the third kind, Hankel functions.
      static int BESSEL_I
      modified Bessel function of the first kind.
      static int BESSEL_J
      Bessel function of the first kind.
      static int BESSEL_K
      modified Bessel function of the second kind.
      static int BESSEL_Y
      Bessel function of the second kind.
      private int BesselType
      DOCUMENT ME!
      private double[] cyi
      DOCUMENT ME!
      private double[] cyr
      DOCUMENT ME!
      private int derivativeOrder
      order of derivative, 0 or 1.
      private boolean doTest
      DOCUMENT ME!
      private int doubleDigits
      I1MACH(14) = 53.
      private double elim
      DOCUMENT ME!
      private int emax
      I1MACH(16) = 1024 The largest exponent for double precision.
      private int emin
      I1MACH(15) = -1021; The smallest exponent for double precision.
      private double epsilon
      D1MACH(4).
      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.
      private double fnul
      DOCUMENT ME!
      private int HankelKind
      real and imaginary outputs Length of arrays passed in must be sequenceNumber Kind of Kankel function = 1 or 2.
      private double initialOrder
      order of initial function.
      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.
      private boolean overflowTest
      DOCUMENT ME!
      private double r1m5
      D1MACH(5) = log10(2).
      private double rl
      DOCUMENT ME!
      static int SCALED_FUNCTION
      multiplied by exp(-abs(Real(z))).
      private int scalingOption
      UNSCALED_FUNCTION returns cy[j-1] = besselFunction(initialOrder+j-1,z), j = 1,...
      private int sequenceNumber
      number of members of the sequence.
      private double tiny
      2**-1022 = D1MACH(1).
      private double tol
      DOCUMENT ME!
      static int UNSCALED_FUNCTION
      DOCUMENT ME!
      private double zi
      DOCUMENT ME!
      private double zr
      zr and zi are the real and imaginary parts of the complex argument.

    • Constructor Summary

      Constructors 
      Constructor Description
      Bessel​(int BesselType, boolean overflowTest)
      Creates a new Bessel object.
      Bessel​(int BesselType, double zr, double zi, double initialOrder, int scalingOption, int sequenceNumber, double[] cyr, double[] cyi, int[] nz, int[] errorFlag)
      This construtor used for I, J, K, and Y Bessel functions.
      Bessel​(int BesselType, double zr, double zi, double initialOrder, int scalingOption, int HankelKind, int sequenceNumber, double[] cyr, double[] cyi, int[] nz, int[] errorFlag)
      This constructor used for H functions of kind 1 and 2.
      Bessel​(int BesselType, double zr, double zi, int derivativeOrder, int scalingOption, double[] cyr, double[] cyi, int[] nz, int[] errorFlag)
      This constructor is used for Ai and Bi Airy functions.

    • Method Summary

      All Methods Instance Methods Concrete Methods 
      Modifier and Type Method Description
      private double dgamln​(double z, int[] ierr)
      DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR Z.GT.0.
      void finalize()
      Cleanup memory.
      void run()
      DOCUMENT ME!
      private void runTest()
      DOCUMENT ME!
      private double zabs​(double zr, double zi)
      zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.
      private void zacai​(double zr, double zi, double fnu, int kode, int mr, int n, double[] yr, double[] yi, int[] nz)
      ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA.
      private void zacon​(double zr, double zi, double fnu, int kode, int mr, int n, double[] yr, double[] yi, int[] nz)
      ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA.
      private void zairy​(double zr, double zi, int id, int kode, double[] air, double[] 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.
      private void zasyi​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] 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).
      private void zbesh​(double zr, double zi, double fnu, int kode, int m, int n, double[] cyr, double[] cyi, int[] nz, int[] ierr)
      DOCUMENT ME!
      private void zbesi​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz, int[] ierr)
      DOCUMENT ME!
      private void zbesj​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz, int[] ierr)
      DOCUMENT ME!
      private void zbesk​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz, int[] ierr)
      DOCUMENT ME!
      private void zbesy​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz, double[] cwrkr, double[] cwrki, int[] ierr)
      DOCUMENT ME!
      private void zbesyh​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz, double[] cwrkr, double[] cwrki, int[] ierr)
      An older version of zbesy used by the ZQCBY quick check routine.
      private void zbinu​(double zr, double zi, double fnu, int kode, int n, double[] cyr, double[] cyi, int[] nz)
      zbinu calculates the I function in the right half z plane.
      private void zbiry​(double zr, double zi, int id, int kode, double[] bir, double[] bii, int[] ierr)
      DOCUMENT ME!
      private void zbknu​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz)
      zbknu computes the k Bessel function in the right half z plane.
      private void zbuni​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz, int nui, int[] nlast, double fnul)
      ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
      private void zbunk​(double zr, double zi, double fnu, int kode, int mr, int n, double[] yr, double[] yi, int[] nz)
      ZBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU.GT.FNUL.
      private void zdiv​(double ar, double ai, double br, double bi, double[] cr, double[] ci)
      complex divide c = a/b.
      private void zexp​(double ar, double ai, double[] br, double[] bi)
      complex exponential function b = exp(a).
      private void zkscl​(double zrr, double zri, double fnu, int n, double[] yr, double[] yi, int[] nz, double rzr, double rzi, double 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.
      private void zlog​(double ar, double ai, double[] br, double[] bi, int[] ierr)
      complex logarithm b = clog(a).
      private void zmlri​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz)
      zmlri computes the I Bessel function for the Real(z) >= 0.0 by the Miller Algorithm normalized by a Neumann series.
      private void zmlt​(double ar, double ai, double br, double bi, double[] cr, double[] ci)
      complex multiply c = a * b.
      private void zqcai​(int mqc)
      DOCUMENT ME!
      private void zqcbh​(int mqc)
      DOCUMENT ME!
      private void zqcbi​(int mqc)
      DOCUMENT ME!
      private void zqcbj​(int mqc)
      DOCUMENT ME!
      private void zqcbk​(int mqc)
      DOCUMENT ME!
      private void zqcby​(int mqc)
      DOCUMENT ME!
      private void zrati​(double zr, double zi, double fnu, int n, double[] cyr, double[] cyi)
      ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD RECURRENCE.
      private void zs1s2​(double zrr, double zri, double[] s1r, double[] s1i, double[] s2r, double[] s2i, int[] nz, double 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.
      private void zseri​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] 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).
      private void zshch​(double zr, double zi, double[] cshr, double[] cshi, double[] cchr, double[] cchi)
      zshch computes the complex hyperbolic functions csh = sinh(x+i*y) and cch = cosh(x+i*y).
      private void zsqrt​(double ar, double ai, double[] br, double[] bi)
      complex square root b = csqrt(a).
      private void zuchk​(double yr, double yi, int[] nz, double ascle)
      Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN EXP(-ALIM)=ASCLE=1.0E+3*D1MACH(1)/TOL.
      private void zunhj​(double zr, double zi, double fnu, int ipmtr, double[] phir, double[] phii, double[] argr, double[] argi, double[] zeta1r, double[] zeta1i, double[] zeta2r, double[] zeta2i, double[] asumr, double[] asumi, double[] bsumr, double[] bsumi)
      REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M.
      private void zuni1​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz, int[] nlast, double 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.
      private void zuni2​(double zr, double zi, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz, int[] nlast, double 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.
      private void zunik​(double zrr, double zri, double fnu, int ikflg, int ipmtr, int[] init, double[] phir, double[] phii, double[] zeta1r, double[] zeta1i, double[] zeta2r, double[] zeta2i, double[] sumr, double[] sumi, double[] cwrkr, double[] cwrki)
      ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 RESPECTIVELY BY.
      private void zunk1​(double zr, double zi, double fnu, int kode, int mr, int n, double[] yr, double[] 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.
      private void zunk2​(double zr, double zi, double fnu, int kode, int mr, int n, double[] yr, double[] 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.
      private void zuoik​(double zr, double zi, double fnu, int kode, int ikflg, int n, double[] yr, double[] 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.
      private void zwrsk​(double zrr, double zri, double fnu, int kode, int n, double[] yr, double[] yi, int[] nz, double[] cwr, double[] cwi)
      zwrsk computes the I Bessel function for Real(z) >- 0.0 by normalizing the I function ratios from zrati by the Wronskian.

      • Methods inherited from class java.lang.Object

        clone, equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

    • 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

      • alim

        private double alim
        DOCUMENT ME!

      • BesselType

        private final int BesselType
        DOCUMENT ME!

      • cyi

        private double[] cyi
        DOCUMENT ME!

      • cyr

        private double[] cyr
        DOCUMENT ME!

      • derivativeOrder

        private int derivativeOrder
        order of derivative, 0 or 1.

      • doTest

        private boolean doTest
        DOCUMENT ME!

      • doubleDigits

        private int doubleDigits
        I1MACH(14) = 53.

      • elim

        private double elim
        DOCUMENT ME!

      • emax

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

      • emin

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

      • epsilon

        private double epsilon
        D1MACH(4).

      • 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.

      • fnul

        private double fnul
        DOCUMENT ME!

      • 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 double r1m5
        D1MACH(5) = log10(2).

      • rl

        private double rl
        DOCUMENT ME!

      • 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.

      • tiny

        private double tiny
        2**-1022 = D1MACH(1).

      • tol

        private double tol
        DOCUMENT ME!

      • 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

      • Bessel

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

      • Bessel

        public Bessel​(int BesselType,
              double zr,
              double zi,
              int derivativeOrder,
              int scalingOption,
              double[] cyr,
              double[] 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 - double[]
        cyi - double[]
        nz - int[]
        errorFlag - int[]

      • Bessel

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

      • Bessel

        public Bessel​(int BesselType,
              double zr,
              double zi,
              double initialOrder,
              int scalingOption,
              int HankelKind,
              int sequenceNumber,
              double[] cyr,
              double[] 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 - double[]
        cyi - double[]
        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 double dgamln​(double 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 - double, z > 0.0
        ierr - int[] error flag IERR=0, NORMAL RETURN, COMPUTATION COMPLETED IERR=1, Z.LE.0.0D0, NO COMPUTATION
        Returns:
        double NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0

      • runTest

        private void runTest()
        DOCUMENT ME!

      • zabs

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

      • zacai

        private void zacai​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int mr,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        mr - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zacon

        private void zacon​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int mr,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        mr - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zairy

        private void zairy​(double zr,
                   double zi,
                   int id,
                   int kode,
                   double[] air,
                   double[] 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 - double
        zi - double
        id - int
        kode - int
        air - double[]
        aii - double[]
        nz - int[]
        ierr - int[]

      • zasyi

        private void zasyi​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zbesh

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

      • zbesi

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

      • zbesj

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

      • zbesk

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

      • zbesy

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

      • zbesyh

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

      • zbinu

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

      • zbiry

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

      • zbknu

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

      • zbuni

        private void zbuni​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] yi,
                   int[] nz,
                   int nui,
                   int[] nlast,
                   double 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]
        nui - int
        nlast - int[]
        fnul - double

      • zbunk

        private void zbunk​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int mr,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        mr - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zdiv

        private void zdiv​(double ar,
                  double ai,
                  double br,
                  double bi,
                  double[] cr,
                  double[] ci)
        complex divide c = a/b.
        Parameters:
        ar - double
        ai - double
        br - double
        bi - double
        cr - double[]
        ci - double[]

      • zexp

        private void zexp​(double ar,
                  double ai,
                  double[] br,
                  double[] bi)
        complex exponential function b = exp(a).
        Parameters:
        ar - double
        ai - double
        br - double[]
        bi - double[]

      • zkscl

        private void zkscl​(double zrr,
                   double zri,
                   double fnu,
                   int n,
                   double[] yr,
                   double[] yi,
                   int[] nz,
                   double rzr,
                   double rzi,
                   double 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 - double
        zri - double
        fnu - double
        n - int
        yr - double[]
        yi - double[]
        nz - int[]
        rzr - double
        rzi - double
        ascle - double

      • zlog

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

      • zmlri

        private void zmlri​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zmlt

        private void zmlt​(double ar,
                  double ai,
                  double br,
                  double bi,
                  double[] cr,
                  double[] ci)
        complex multiply c = a * b.
        Parameters:
        ar - double
        ai - double
        br - double
        bi - double
        cr - double[]
        ci - double[]

      • 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​(double zr,
                   double zi,
                   double fnu,
                   int n,
                   double[] cyr,
                   double[] 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 - double
        zi - double
        fnu - double
        n - int
        cyr - double[]
        cyi - double[]

      • zs1s2

        private void zs1s2​(double zrr,
                   double zri,
                   double[] s1r,
                   double[] s1i,
                   double[] s2r,
                   double[] s2i,
                   int[] nz,
                   double 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 - double
        zri - double
        s1r - double[]
        s1i - double[]
        s2r - double[]
        s2i - double[]
        nz - int[]
        ascle - double
        iuf - int[]

      • zseri

        private void zseri​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zshch

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

      • zsqrt

        private void zsqrt​(double ar,
                   double ai,
                   double[] br,
                   double[] bi)
        complex square root b = csqrt(a).
        Parameters:
        ar - double
        ai - double
        br - double[]
        bi - double[]

      • zuchk

        private void zuchk​(double yr,
                   double yi,
                   int[] nz,
                   double 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 - double
        yi - double
        nz - int[]
        ascle - double

      • zunhj

        private void zunhj​(double zr,
                   double zi,
                   double fnu,
                   int ipmtr,
                   double[] phir,
                   double[] phii,
                   double[] argr,
                   double[] argi,
                   double[] zeta1r,
                   double[] zeta1i,
                   double[] zeta2r,
                   double[] zeta2i,
                   double[] asumr,
                   double[] asumi,
                   double[] bsumr,
                   double[] 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 - double
        zi - double
        fnu - double
        ipmtr - int
        phir - double[]
        phii - double[]
        argr - double[]
        argi - double[]
        zeta1r - double[]
        zeta1i - double[]
        zeta2r - double[]
        zeta2i - double[]
        asumr - double[]
        asumi - double[]
        bsumr - double[]
        bsumi - double[]

      • zuni1

        private void zuni1​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] yi,
                   int[] nz,
                   int[] nlast,
                   double 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]
        nlast - int[]
        fnul - double

      • zuni2

        private void zuni2​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] yi,
                   int[] nz,
                   int[] nlast,
                   double 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 - double
        zi - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]
        nlast - int[]
        fnul - double

      • zunik

        private void zunik​(double zrr,
                   double zri,
                   double fnu,
                   int ikflg,
                   int ipmtr,
                   int[] init,
                   double[] phir,
                   double[] phii,
                   double[] zeta1r,
                   double[] zeta1i,
                   double[] zeta2r,
                   double[] zeta2i,
                   double[] sumr,
                   double[] sumi,
                   double[] cwrkr,
                   double[] 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 - double
        zri - double
        fnu - double
        ikflg - int
        ipmtr - int
        init - int
        phir - double[]
        phii - double[]
        zeta1r - double[]
        zeta1i - double[]
        zeta2r - double[]
        zeta2i - double[]
        sumr - double[]
        sumi - double[]
        cwrkr - double[]
        cwrki - double[]

      • zunk1

        private void zunk1​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int mr,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        mr - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zunk2

        private void zunk2​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int mr,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        mr - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]

      • zuoik

        private void zuoik​(double zr,
                   double zi,
                   double fnu,
                   int kode,
                   int ikflg,
                   int n,
                   double[] yr,
                   double[] 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 - double
        zi - double
        fnu - double
        kode - int
        ikflg - int
        n - int
        yr - double[]
        yi - double[]
        nuf - int[]

      • zwrsk

        private void zwrsk​(double zrr,
                   double zri,
                   double fnu,
                   int kode,
                   int n,
                   double[] yr,
                   double[] yi,
                   int[] nz,
                   double[] cwr,
                   double[] 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 - double
        zri - double
        fnu - double
        kode - int
        n - int
        yr - double[]
        yi - double[]
        nz - int[]
        cwr - double[]
        cwi - double[]