Package gov.nih.mipav.model.algorithms

Class Bessel

java.lang.Object
gov.nih.mipav.model.algorithms.Bessel
public class Bessel extends 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 invalid input: '&' 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 final int
    AIRY_AI
    Airy function Ai.
    static final int
    AIRY_BI
    Airy function Bi.
    private double
    alim
    DOCUMENT ME!
    static final int
    BESSEL_H
    Bessel functions of the third kind, Hankel functions.
    static final int
    BESSEL_I
    modified Bessel function of the first kind.
    static final int
    BESSEL_J
    Bessel function of the first kind.
    static final int
    BESSEL_K
    modified Bessel function of the second kind.
    static final int
    BESSEL_Y
    Bessel function of the second kind.
    private final 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 final 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,...,sequenceNumber SCALED_FUNCTION returns cy[j-1] = besselFunction(initialOrder+j-1,z)*exp(-abs(Real(z))), j = 1,...,sequenceNumber.
    private int
    sequenceNumber
    number of members of the sequence.
    private double
    tiny
    2**-1022 = D1MACH(1).
    private double
    tol
    DOCUMENT ME!
    static final 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

    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 Details

    • BESSEL_H

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

    • BESSEL_I

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

    • BESSEL_J

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

    • BESSEL_K

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

    • BESSEL_Y

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

    • AIRY_AI

      public static final int AIRY_AI
      Airy function Ai.
      See Also:

    • AIRY_BI

      public static final int AIRY_BI
      Airy function Bi.
      See Also:

    • UNSCALED_FUNCTION

      public static final int UNSCALED_FUNCTION
      DOCUMENT ME!
      See Also:

    • SCALED_FUNCTION

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

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

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

    • finalize

      public void finalize() throws Throwable
      Cleanup memory.
      Overrides:
      finalize in class Object
      Throws:
      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[]