Remarks on the numerical implementation

The solution of Eq. (2.7) requires the evaluation of Bessel and Hankel functions of complex order, where typically values of $\nu = \gamma R$ in intervals $\Re(\nu) \in [10^1, 10^4]$ and $\Im(\nu) \in [-10^{-1}, -10^{-12}]$ and arguments $nkr$ up to $10^{5}$ are encountered. The relevant Bessel functions are well behaved for the above values of order and argument; in particular, no terms with numerically harmful growth appear. We found that subroutines for complex order Bessel functions are not included in standard numerical libraries. MATLAB includes routines for these functions with complex argument with real order. MAPLE has built in procedures for complex argument and complex order, but they turn out to be prohibitively slow for large orders. Therefore we had to resort to own implementations.

The implemented procedures are based on ``uniform asymptotic expansions'' of Bessel functions and their derivatives in terms of Airy functions [73,74]. More specifically, Eqs. (9.3.35), (9.3.36), (9.3.43), (9.3.44) from Ref. [73] were encoded, restricted to the first two terms of the summations, which we observed to be sufficient for the present examples. Hankel functions are computed as linear combination of Bessel functions. These 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 ($\gamma R$) with negative imaginary parts ($-\alpha$), and real arguments ($nkr$). While cases with approximate equality of order and argument could in principle occur for bend modes with extremely low loss (see e.g. bent waveguide configurations in Section 4.4.3.), by assigning a negligible small value to imaginary part of the argument, we circumvent any numerical problems associated with this transition region. Otherwise, supplementary routines (see e.g. [74]) would have to be incorporated that cover that region of parameters.

Routines for Airy functions with complex arguments according to Ref. [92] were adopted. The procedures concerning the bent modes, as well as all parts of resonators model discussed in Chapters 3, 4, were encoded with C++ programming language. Since the Airy function routines were available in FORTRAN 77, this choice lead to using FORTRAN 77 code within the C++ code. By comparison with the numerical values obtained with MAPLE, and also by validating the relations that held for the Wronskians of these functions (see Ref. [73], Eqs. (9.1.15), (9.1.16), (9.1.17)), the performance of the implemented routines for the desired range of orders and arguments is assessed.

Two techniques were implemented to find the roots $\nu$ of Eq. (2.7) in the complex plane. As a heuristic search procedure, a suitably selected region of the complex plane is divided into a number of rectangles, the lower left and upper right corner points of which are then supplied as initial guesses to a root finding routine based on the secant method. Repeated roots are rejected, the remaining unique roots are sorted in descending order of their imaginary part $-\alpha$. By refining the subdivision into rectangles, it can be ascertained with reasonable robustness that all the roots in the given region are captured.

Alternatively, a root tracking procedure can be implemented to solve the dispersion equation. Starting with the propagation constants of straight waveguides with a refractive index profile cross section equal to that of the given bent waveguide, a series of bends with decreasing radius are considered, with the roots found for each configuration used as initial guesses for the subsequent one. In this way, the bend propagation constants are followed in the complex plane.