Computation of complex order Bessel functions

Introduction

• To solve the dispersion equation (see Ref.[1, Eq. 7]), one requires evaluation of Bessel and Hankel functions of complex order. We have the required routines, which are based on 'uniform asymptotic expansions' of Bessel functions and their derivatives in terms of Airy functions (Ref. [2],[3]).
  More specifically, Eqs. 9.3.35, 9.3.36, 9.3.43, 9.3.44 from Ref.[2] were encoded, restricted to the first two terms of the summations, which we observed to be sufficient for the present examples. Routines for Airy functions with complex arguments according to Ref. ,[4]) were adopted. Further details can be found in Ref. [1].

• These are double precision routines, and can be used for real/complex order and real/complex argument Bessel/Hankel functions.

• The expansions are not applicable in a regime where the order is close to the argument of the cylindrical functions. Fortunately we are interested here in configurations that involve complex orders γ R with negative imaginary parts -α, and real arguments n k r. While cases with approximate equality of order and argument could in principle occur for bend modes with extremely low loss, for the evaluation of propagation constants and the profile plots of the examples below this transition region appeared not to be relevant. Otherwise, supplementary routines (see e.g. Ref. [4]) for details) would have to be incorporated that cover that region of parameters.

Source files

Define

• #define Pi M_PI
• #define nan sqrt(-1.0)

Functions

Description of variables and functions

• #define Pi M_PI
  Constant Pi

• #define nan sqrt(-1.0)
  Not-a-Number

• void bessel_error(char *s)
  Exception handling function.

• extern complex<double> BesselJ( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> BesselY( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> HankelH1( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> HankelH2( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  Compute Bessel function of first kind (BesselJ), Bessel function of second kind (BesselY), Hankel function of first kind (HankelH1), and Hankel funciton of second kind (HankelH2) respectively.
  

  order	= order of a function
  argument	= argument of a function
  KODE	= scaling parameter
  	= 1 without scaling
  	= 2 with scaling
  IERR	= Error code
  	= 0 Success
  	= 1 Inout error
  	= 2 Overflow. No computation. Re(zeta) too large for KODE = 1 (for definition of zeta, see BesselJ() routine code).
  	= 3 |z| too large. Answer may be wrong, where z = argument/order.
  	= 4 |z| too large, where z = argument/order. No computation.
  	= 5 No computation, algorithm termination condition did not meet
  	= 11 Situation beyond scope of routines
  	= 12 Division by zero

• extern complex<double> DBesselJ( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> DBesselY( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> DHankelH1( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  extern complex<double> DHankelH2( complex<double> order, complex<double> argument, int KODE, int *NZ, int *IERR)
  Compute first order derivative with respect to argument of a Bessel function of first kind (BesselJ), Bessel function of second kind (BesselY), Hankel function of first kind (HankelH1), and Hankel funciton of second kind (HankelH2) respectively.
  Explaination of varibales is same as that of for corresponding Bessel/Hankel funcions.

• double lambda(int s)
  double mu(int s)
  complex<double> u(int k, complex<double> t)
  complex<double> v(int k, complex<double> t)
  complex<double> a(int k, complex<double> zeta)
  complex<double> b(int k, complex<double> zeta)
  complex<double> c(int k, complex<double> zeta, complex<double> z)
  complex<double> d(int k, complex<double> zeta, complex<double> z)
  Helper functions for computation of Bessel functions. See [2, Chapt. 9] for details.

• extern "C"
  {
  extern void zairy_( double *zr, double *zi, int *ID, int *KODE, double *air, double *aii, int *NZ, int *IERR);
  extern void zbiry_( double *zr, double *zi, int *ID, int *KODE, double *air, double *aii, int *IERR);
  extern double dgamma_(double *x);
  }
  FORTRAN routines to compute complex Airy functions AiryAi (zairy_), AiryBi (zbiry_) and Gamma function (dgamma_).

References

1. K. R. Hiremath, M. Hammer, S. Stoffer, L. Prkna, and J. Ctyroký.
   Analytic approach to dielectric optical bent slab waveguides.
   Optical and Quantum Electronics, 37(1-3):37-61, January 2005. PDF
2. M. Abramowitz and I. A. Stegun.
   Handbook of Mathematical Functions (Applied Mathematics Series 55).
   National Bureau of Standards, Washington, D.C., 1964.
3. N. M. Temme.
   Numerical algorithms for uniform Airy-type asymptotic expansions.
   Technical Report MAS-R9706, Centrum voor Wiskunde en Informatica, Amsterdam, The Netherlands, 1997.
4. D. E. Amos.
   A portable package for Bessel functions of a complex argument and nonnegative order, 1983.
   http://www.netlib.org/amos/.