amos.h

/*
 *
 * This file is a C++ translation of the Fortran code written by
 * D.E. Amos with the following original description:
 *
 *
 * A Portable Package for Bessel Functions of a Complex Argument
 * and Nonnegative Order
 *
 * This algorithm is a package of subroutines for computing Bessel
 * functions and Airy functions.  The routines are updated
 * versions of those routines found in TOMS algorithm 644.
 *
 * Disclaimer:
 *
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *                 ISSUED BY SANDIA LABORATORIES,
 *                   A PRIME CONTRACTOR TO THE
 *               UNITED STATES DEPARTMENT OF ENERGY
 * * * * * * * * * * * * * *  NOTICE   * * * * * * * * * * * * * * *
 * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
 * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
 * UNITED STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR
 * EMPLOYEES, NOR ANY OF THEIR CONTRACTORS, SUBCONTRACTORS, OR THEIR
 * EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
 * LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS
 * OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT OR PROCESS
 * DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE
 * PRIVATELY OWNED RIGHTS.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * THIS CODE HAS BEEN APPROVED FOR UNLIMITED RELEASE.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *
 *
 *
 * The original Fortran code can be found at https://www.netlib.org/amos/
 *
 * 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, DOI:10.1145/7921.214331
 *
 * [3]: Amos, D. E., Remark on Algorithm 644, ACM Transactions on
 *      Mathematical Software, Vol. 16, No. 4, December 1990, Page
 *      404, DOI:10.1145/98267.98299
 *
 * [4]: Amos, D. E., A remark on Algorithm 644: "A portable package
 *      for Bessel functions of a complex argument and nonnegative order",
 *      ACM Transactions on Mathematical Software, Vol. 21, No. 4,
 *      December 1995, Pages 388-393, DOI:10.1145/212066.212078
 *
 * [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, DOI:10.1145/50063.51907
 *
 */

/*
 * Copyright (C) 2024 SciPy developers
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * a. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 * b. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * c. Names of the SciPy Developers may not be used to endorse or promote
 *    products derived from this software without specific prior written
 *    permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */
#pragma once

#include <cmath>
#include <complex>
#include <cstdlib>
#include <memory> // unique_ptr

namespace xsf {
namespace amos {

    int acai(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double, double);
    int
    acon(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double, double, double);
    int asyi(std::complex<double>, double, int, int, std::complex<double> *, double, double, double, double);
    int
    binu(std::complex<double>, double fnu, int, int, std::complex<double> *, double, double, double, double, double);
    int bknu(std::complex<double>, double, int, int, std::complex<double> *, double, double, double);
    int
    buni(std::complex<double>, double, int, int, std::complex<double> *, int, int *, double, double, double, double);
    int bunk(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double);
    double gamln(double);
    int kscl(std::complex<double>, double, int, std::complex<double> *, std::complex<double>, double *, double, double);
    int mlri(std::complex<double>, double, int, int, std::complex<double> *, double);
    void rati(std::complex<double>, double, int, std::complex<double> *, double);
    int seri(std::complex<double>, double, int, int, std::complex<double> *, double, double, double);
    int s1s2(std::complex<double>, std::complex<double> *, std::complex<double> *, double, double, int *);
    int uchk(std::complex<double>, double, double);
    void unhj(
        std::complex<double>, double, int, double, std::complex<double> *, std::complex<double> *,
        std::complex<double> *, std::complex<double> *, std::complex<double> *, std::complex<double> *
    );
    void
    uni1(std::complex<double>, double, int, int, std::complex<double> *, int *, int *, double, double, double, double);
    void
    uni2(std::complex<double>, double, int, int, std::complex<double> *, int *, int *, double, double, double, double);
    void unik(
        std::complex<double>, double, int, int, double, int *, std::complex<double> *, std::complex<double> *,
        std::complex<double> *, std::complex<double> *, std::complex<double> *
    );
    int unk1(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double);
    int unk2(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double);
    int uoik(std::complex<double>, double, int, int, int, std::complex<double> *, double, double, double);
    int wrsk(
        std::complex<double>, double, int, int, std::complex<double> *, std::complex<double> *, double, double, double
    );

    constexpr double d1mach[5] = {
        2.2250738585072014e-308, /* np.finfo(np.float64).tiny      */
        1.7976931348623157e+308, /* np.finfo(np.float64).max       */
        1.1102230246251565e-16,  /* 0.5 * np.finfo(np.float64).eps */
        2.220446049250313e-16,   /* np.finfo(np.float64).eps       */
        0.3010299956639812       /* np.log10(2)                    */
    };

    constexpr double i1mach[16] = {
        5,          /* standard input         */
        6,          /* standard output        */
        7,          /* standard punch         */
        0,          /* standard error         */
        32,         /* bits per integer       */
        4,          /* sizeof(int);           */
        2,          /* base for integers      */
        31,         /* digits of integer base */
        2147483647, /* LONG MAX 2**31 - 1     */
        2,          /* FLT_RADIX;             */
        24,         /* FLT_MANT_DIG;          */
        -126,       /* FLT_MIN_EXP;           */
        128,        /* FLT_MAX_EXP;           */
        53,         /* DBL_MANT_DIG;          */
        -1021,      /* DBL_MIN_EXP;           */
        1024        /* DBL_MAX_EXP;           */
    };

    constexpr double zunhj_ar[14] = {
        1.00000000000000000e+00, 1.04166666666666667e-01, 8.35503472222222222e-02, 1.28226574556327160e-01, //  0
        2.91849026464140464e-01, 8.81627267443757652e-01, 3.32140828186276754e+00, 1.49957629868625547e+01, //  4
        7.89230130115865181e+01, 4.74451538868264323e+02, 3.20749009089066193e+03, 2.40865496408740049e+04, //  8
        1.98923119169509794e+05, 1.79190200777534383e+06                                                    // 12
    };

    constexpr double zunhj_br[14] = {
        1.00000000000000000e+00,  -1.45833333333333333e-01, -9.87413194444444444e-02, -1.43312053915895062e-01, //  0
        -3.17227202678413548e-01, -9.42429147957120249e-01, -3.51120304082635426e+00, -1.57272636203680451e+01, //  4
        -8.22814390971859444e+01, -4.92355370523670524e+02, -3.31621856854797251e+03, -2.48276742452085896e+04, //  8
        -2.04526587315129788e+05, -1.83844491706820990e+06                                                      // 12
    };

    constexpr double zunhj_c[105] = {
        1.00000000000000000e+00,  -2.08333333333333333e-01, 1.25000000000000000e-01,  3.34201388888888889e-01,  //   0
        -4.01041666666666667e-01, 7.03125000000000000e-02,  -1.02581259645061728e+00, 1.84646267361111111e+00,  //   4
        -8.91210937500000000e-01, 7.32421875000000000e-02,  4.66958442342624743e+00,  -1.12070026162229938e+01, //   8
        8.78912353515625000e+00,  -2.36408691406250000e+00, 1.12152099609375000e-01,  -2.82120725582002449e+01, //  12
        8.46362176746007346e+01,  -9.18182415432400174e+01, 4.25349987453884549e+01,  -7.36879435947963170e+00, //  16
        2.27108001708984375e-01,  2.12570130039217123e+02,  -7.65252468141181642e+02, 1.05999045252799988e+03,  //  20
        -6.99579627376132541e+02, 2.18190511744211590e+02,  -2.64914304869515555e+01, 5.72501420974731445e-01,  //  24
        -1.91945766231840700e+03, 8.06172218173730938e+03,  -1.35865500064341374e+04, 1.16553933368645332e+04,  //  28
        -5.30564697861340311e+03, 1.20090291321635246e+03,  -1.08090919788394656e+02, 1.72772750258445740e+00,  //  32
        2.02042913309661486e+04,  -9.69805983886375135e+04, 1.92547001232531532e+05,  -2.03400177280415534e+05, //  36
        1.22200464983017460e+05,  -4.11926549688975513e+04, 7.10951430248936372e+03,  -4.93915304773088012e+02, //  40
        6.07404200127348304e+00,  -2.42919187900551333e+05, 1.31176361466297720e+06,  -2.99801591853810675e+06, //  44
        3.76327129765640400e+06,  -2.81356322658653411e+06, 1.26836527332162478e+06,  -3.31645172484563578e+05, //  48
        4.52187689813627263e+04,  -2.49983048181120962e+03, 2.43805296995560639e+01,  3.28446985307203782e+06,  //  52
        -1.97068191184322269e+07, 5.09526024926646422e+07,  -7.41051482115326577e+07, 6.63445122747290267e+07,  //  56
        -3.75671766607633513e+07, 1.32887671664218183e+07,  -2.78561812808645469e+06, 3.08186404612662398e+05,  //  60
        -1.38860897537170405e+04, 1.10017140269246738e+02,  -4.93292536645099620e+07, 3.25573074185765749e+08,  //  64
        -9.39462359681578403e+08, 1.55359689957058006e+09,  -1.62108055210833708e+09, 1.10684281682301447e+09,  //  68
        -4.95889784275030309e+08, 1.42062907797533095e+08,  -2.44740627257387285e+07, 2.24376817792244943e+06,  //  72
        -8.40054336030240853e+04, 5.51335896122020586e+02,  8.14789096118312115e+08,  -5.86648149205184723e+09, //  76
        1.86882075092958249e+10,  -3.46320433881587779e+10, 4.12801855797539740e+10,  -3.30265997498007231e+10, //  80
        1.79542137311556001e+10,  -6.56329379261928433e+09, 1.55927986487925751e+09,  -2.25105661889415278e+08, //  84
        1.73951075539781645e+07,  -5.49842327572288687e+05, 3.03809051092238427e+03,  -1.46792612476956167e+10, //  88
        1.14498237732025810e+11,  -3.99096175224466498e+11, 8.19218669548577329e+11,  -1.09837515608122331e+12, //  92
        1.00815810686538209e+12,  -6.45364869245376503e+11, 2.87900649906150589e+11,  -8.78670721780232657e+10, //  96
        1.76347306068349694e+10,  -2.16716498322379509e+09, 1.43157876718888981e+08,  -3.87183344257261262e+06, // 100
        1.82577554742931747e+04                                                                                 // 104
    };

    constexpr double zunhj_alfa[180] = {
        -4.44444444444444444e-03, -9.22077922077922078e-04, -8.84892884892884893e-05, 1.65927687832449737e-04,  //   0
        2.46691372741792910e-04,  2.65995589346254780e-04,  2.61824297061500945e-04,  2.48730437344655609e-04,  //   4
        2.32721040083232098e-04,  2.16362485712365082e-04,  2.00738858762752355e-04,  1.86267636637545172e-04,  //   8
        1.73060775917876493e-04,  1.61091705929015752e-04,  1.50274774160908134e-04,  1.40503497391269794e-04,  //  12
        1.31668816545922806e-04,  1.23667445598253261e-04,  1.16405271474737902e-04,  1.09798298372713369e-04,  //  16
        1.03772410422992823e-04,  9.82626078369363448e-05,  9.32120517249503256e-05,  8.85710852478711718e-05,  //  20
        8.42963105715700223e-05,  8.03497548407791151e-05,  7.66981345359207388e-05,  7.33122157481777809e-05,  //  24
        7.01662625163141333e-05,  6.72375633790160292e-05,  6.93735541354588974e-04,  2.32241745182921654e-04,  //  28
        -1.41986273556691197e-05, -1.16444931672048640e-04, -1.50803558053048762e-04, -1.55121924918096223e-04, //  32
        -1.46809756646465549e-04, -1.33815503867491367e-04, -1.19744975684254051e-04, -1.06184319207974020e-04, //  36
        -9.37699549891194492e-05, -8.26923045588193274e-05, -7.29374348155221211e-05, -6.44042357721016283e-05, //  40
        -5.69611566009369048e-05, -5.04731044303561628e-05, -4.48134868008882786e-05, -3.98688727717598864e-05, //  44
        -3.55400532972042498e-05, -3.17414256609022480e-05, -2.83996793904174811e-05, -2.54522720634870566e-05, //  48
        -2.28459297164724555e-05, -2.05352753106480604e-05, -1.84816217627666085e-05, -1.66519330021393806e-05, //  52
        -1.50179412980119482e-05, -1.35554031379040526e-05, -1.22434746473858131e-05, -1.10641884811308169e-05, //  56
        -3.54211971457743841e-04, -1.56161263945159416e-04, 3.04465503594936410e-05,  1.30198655773242693e-04,  //  60
        1.67471106699712269e-04,  1.70222587683592569e-04,  1.56501427608594704e-04,  1.36339170977445120e-04,  //  64
        1.14886692029825128e-04,  9.45869093034688111e-05,  7.64498419250898258e-05,  6.07570334965197354e-05,  //  68
        4.74394299290508799e-05,  3.62757512005344297e-05,  2.69939714979224901e-05,  1.93210938247939253e-05,  //  72
        1.30056674793963203e-05,  7.82620866744496661e-06,  3.59257485819351583e-06,  1.44040049814251817e-07,  //  76
        -2.65396769697939116e-06, -4.91346867098485910e-06, -6.72739296091248287e-06, -8.17269379678657923e-06, //  80
        -9.31304715093561232e-06, -1.02011418798016441e-05, -1.08805962510592880e-05, -1.13875481509603555e-05, //  84
        -1.17519675674556414e-05, -1.19987364870944141e-05, 3.78194199201772914e-04,  2.02471952761816167e-04,  //  88
        -6.37938506318862408e-05, -2.38598230603005903e-04, -3.10916256027361568e-04, -3.13680115247576316e-04, //  92
        -2.78950273791323387e-04, -2.28564082619141374e-04, -1.75245280340846749e-04, -1.25544063060690348e-04, //  96
        -8.22982872820208365e-05, -4.62860730588116458e-05, -1.72334302366962267e-05, 5.60690482304602267e-06,  // 100
        2.31395443148286800e-05,  3.62642745856793957e-05,  4.58006124490188752e-05,  5.24595294959114050e-05,  // 104
        5.68396208545815266e-05,  5.94349820393104052e-05,  6.06478527578421742e-05,  6.08023907788436497e-05,  // 108
        6.01577894539460388e-05,  5.89199657344698500e-05,  5.72515823777593053e-05,  5.52804375585852577e-05,  // 112
        5.31063773802880170e-05,  5.08069302012325706e-05,  4.84418647620094842e-05,  4.60568581607475370e-05,  // 116
        -6.91141397288294174e-04, -4.29976633058871912e-04, 1.83067735980039018e-04,  6.60088147542014144e-04,  // 120
        8.75964969951185931e-04,  8.77335235958235514e-04,  7.49369585378990637e-04,  5.63832329756980918e-04,  // 124
        3.68059319971443156e-04,  1.88464535514455599e-04,  3.70663057664904149e-05,  -8.28520220232137023e-05, // 128
        -1.72751952869172998e-04, -2.36314873605872983e-04, -2.77966150694906658e-04, -3.02079514155456919e-04, // 132
        -3.12594712643820127e-04, -3.12872558758067163e-04, -3.05678038466324377e-04, -2.93226470614557331e-04, // 136
        -2.77255655582934777e-04, -2.59103928467031709e-04, -2.39784014396480342e-04, -2.20048260045422848e-04, // 140
        -2.00443911094971498e-04, -1.81358692210970687e-04, -1.63057674478657464e-04, -1.45712672175205844e-04, // 144
        -1.29425421983924587e-04, -1.14245691942445952e-04, 1.92821964248775885e-03,  1.35592576302022234e-03,  // 148
        -7.17858090421302995e-04, -2.58084802575270346e-03, -3.49271130826168475e-03, -3.46986299340960628e-03, // 152
        -2.82285233351310182e-03, -1.88103076404891354e-03, -8.89531718383947600e-04, 3.87912102631035228e-06,  // 156
        7.28688540119691412e-04,  1.26566373053457758e-03,  1.62518158372674427e-03,  1.83203153216373172e-03,  // 160
        1.91588388990527909e-03,  1.90588846755546138e-03,  1.82798982421825727e-03,  1.70389506421121530e-03,  // 164
        1.55097127171097686e-03,  1.38261421852276159e-03,  1.20881424230064774e-03,  1.03676532638344962e-03,  // 168
        8.71437918068619115e-04,  7.16080155297701002e-04,  5.72637002558129372e-04,  4.42089819465802277e-04,  // 172
        3.24724948503090564e-04,  2.20342042730246599e-04,  1.28412898401353882e-04,  4.82005924552095464e-05   // 176
    };

    constexpr double zunhj_beta[210] = {
        1.79988721413553309e-02,  5.59964911064388073e-03,  2.88501402231132779e-03,  1.80096606761053941e-03,  //   0
        1.24753110589199202e-03,  9.22878876572938311e-04,  7.14430421727287357e-04,  5.71787281789704872e-04,  //   4
        4.69431007606481533e-04,  3.93232835462916638e-04,  3.34818889318297664e-04,  2.88952148495751517e-04,  //   8
        2.52211615549573284e-04,  2.22280580798883327e-04,  1.97541838033062524e-04,  1.76836855019718004e-04,  //  12
        1.59316899661821081e-04,  1.44347930197333986e-04,  1.31448068119965379e-04,  1.20245444949302884e-04,  //  16
        1.10449144504599392e-04,  1.01828770740567258e-04,  9.41998224204237509e-05,  8.74130545753834437e-05,  //  20
        8.13466262162801467e-05,  7.59002269646219339e-05,  7.09906300634153481e-05,  6.65482874842468183e-05,  //  24
        6.25146958969275078e-05,  5.88403394426251749e-05,  -1.49282953213429172e-03, -8.78204709546389328e-04, //  28
        -5.02916549572034614e-04, -2.94822138512746025e-04, -1.75463996970782828e-04, -1.04008550460816434e-04, //  32
        -5.96141953046457895e-05, -3.12038929076098340e-05, -1.26089735980230047e-05, -2.42892608575730389e-07, //  36
        8.05996165414273571e-06,  1.36507009262147391e-05,  1.73964125472926261e-05,  1.98672978842133780e-05,  //  40
        2.14463263790822639e-05,  2.23954659232456514e-05,  2.28967783814712629e-05,  2.30785389811177817e-05,  //  44
        2.30321976080909144e-05,  2.28236073720348722e-05,  2.25005881105292418e-05,  2.20981015361991429e-05,  //  48
        2.16418427448103905e-05,  2.11507649256220843e-05,  2.06388749782170737e-05,  2.01165241997081666e-05,  //  52
        1.95913450141179244e-05,  1.90689367910436740e-05,  1.85533719641636667e-05,  1.80475722259674218e-05,  //  56
        5.52213076721292790e-04,  4.47932581552384646e-04,  2.79520653992020589e-04,  1.52468156198446602e-04,  //  60
        6.93271105657043598e-05,  1.76258683069991397e-05,  -1.35744996343269136e-05, -3.17972413350427135e-05, //  64
        -4.18861861696693365e-05, -4.69004889379141029e-05, -4.87665447413787352e-05, -4.87010031186735069e-05, //  68
        -4.74755620890086638e-05, -4.55813058138628452e-05, -4.33309644511266036e-05, -4.09230193157750364e-05, //  72
        -3.84822638603221274e-05, -3.60857167535410501e-05, -3.37793306123367417e-05, -3.15888560772109621e-05, //  76
        -2.95269561750807315e-05, -2.75978914828335759e-05, -2.58006174666883713e-05, -2.41308356761280200e-05, //  80
        -2.25823509518346033e-05, -2.11479656768912971e-05, -1.98200638885294927e-05, -1.85909870801065077e-05, //  84
        -1.74532699844210224e-05, -1.63997823854497997e-05, -4.74617796559959808e-04, -4.77864567147321487e-04, //  88
        -3.20390228067037603e-04, -1.61105016119962282e-04, -4.25778101285435204e-05, 3.44571294294967503e-05,  //  92
        7.97092684075674924e-05,  1.03138236708272200e-04,  1.12466775262204158e-04,  1.13103642108481389e-04,  //  96
        1.08651634848774268e-04,  1.01437951597661973e-04,  9.29298396593363896e-05,  8.40293133016089978e-05,  // 100
        7.52727991349134062e-05,  6.69632521975730872e-05,  5.92564547323194704e-05,  5.22169308826975567e-05,  // 104
        4.58539485165360646e-05,  4.01445513891486808e-05,  3.50481730031328081e-05,  3.05157995034346659e-05,  // 108
        2.64956119950516039e-05,  2.29363633690998152e-05,  1.97893056664021636e-05,  1.70091984636412623e-05,  // 112
        1.45547428261524004e-05,  1.23886640995878413e-05,  1.04775876076583236e-05,  8.79179954978479373e-06,  // 116
        7.36465810572578444e-04,  8.72790805146193976e-04,  6.22614862573135066e-04,  2.85998154194304147e-04,  // 120
        3.84737672879366102e-06,  -1.87906003636971558e-04, -2.97603646594554535e-04, -3.45998126832656348e-04, // 124
        -3.53382470916037712e-04, -3.35715635775048757e-04, -3.04321124789039809e-04, -2.66722723047612821e-04, // 128
        -2.27654214122819527e-04, -1.89922611854562356e-04, -1.55058918599093870e-04, -1.23778240761873630e-04, // 132
        -9.62926147717644187e-05, -7.25178327714425337e-05, -5.22070028895633801e-05, -3.50347750511900522e-05, // 136
        -2.06489761035551757e-05, -8.70106096849767054e-06, 1.13698686675100290e-06,  9.16426474122778849e-06,  // 140
        1.56477785428872620e-05,  2.08223629482466847e-05,  2.48923381004595156e-05,  2.80340509574146325e-05,  // 144
        3.03987774629861915e-05,  3.21156731406700616e-05,  -1.80182191963885708e-03, -2.43402962938042533e-03, // 148
        -1.83422663549856802e-03, -7.62204596354009765e-04, 2.39079475256927218e-04,  9.49266117176881141e-04,  // 152
        1.34467449701540359e-03,  1.48457495259449178e-03,  1.44732339830617591e-03,  1.30268261285657186e-03,  // 156
        1.10351597375642682e-03,  8.86047440419791759e-04,  6.73073208165665473e-04,  4.77603872856582378e-04,  // 160
        3.05991926358789362e-04,  1.60315694594721630e-04,  4.00749555270613286e-05,  -5.66607461635251611e-05, // 164
        -1.32506186772982638e-04, -1.90296187989614057e-04, -2.32811450376937408e-04, -2.62628811464668841e-04, // 168
        -2.82050469867598672e-04, -2.93081563192861167e-04, -2.97435962176316616e-04, -2.96557334239348078e-04, // 172
        -2.91647363312090861e-04, -2.83696203837734166e-04, -2.73512317095673346e-04, -2.61750155806768580e-04, // 176
        6.38585891212050914e-03,  9.62374215806377941e-03,  7.61878061207001043e-03,  2.83219055545628054e-03,  // 180
        -2.09841352012720090e-03, -5.73826764216626498e-03, -7.70804244495414620e-03, -8.21011692264844401e-03, // 184
        -7.65824520346905413e-03, -6.47209729391045177e-03, -4.99132412004966473e-03, -3.45612289713133280e-03, // 188
        -2.01785580014170775e-03, -7.59430686781961401e-04, 2.84173631523859138e-04,  1.10891667586337403e-03,  // 192
        1.72901493872728771e-03,  2.16812590802684701e-03,  2.45357710494539735e-03,  2.61281821058334862e-03,  // 196
        2.67141039656276912e-03,  2.65203073395980430e-03,  2.57411652877287315e-03,  2.45389126236094427e-03,  // 200
        2.30460058071795494e-03,  2.13684837686712662e-03,  1.95896528478870911e-03,  1.77737008679454412e-03,  // 204
        1.59690280765839059e-03,  1.42111975664438546e-03                                                       // 208
    };

    constexpr double zunhj_gama[30] = {
        6.29960524947436582e-01, 2.51984209978974633e-01, 1.54790300415655846e-01, 1.10713062416159013e-01, //  0
        8.57309395527394825e-02, 6.97161316958684292e-02, 5.86085671893713576e-02, 5.04698873536310685e-02, //  4
        4.42600580689154809e-02, 3.93720661543509966e-02, 3.54283195924455368e-02, 3.21818857502098231e-02, //  8
        2.94646240791157679e-02, 2.71581677112934479e-02, 2.51768272973861779e-02, 2.34570755306078891e-02, // 12
        2.19508390134907203e-02, 2.06210828235646240e-02, 1.94388240897880846e-02, 1.83810633800683158e-02, // 16
        1.74293213231963172e-02, 1.65685837786612353e-02, 1.57865285987918445e-02, 1.50729501494095594e-02, // 20
        1.44193250839954639e-02, 1.38184805735341786e-02, 1.32643378994276568e-02, 1.27517121970498651e-02, // 24
        1.22761545318762767e-02, 1.18338262398482403e-02                                                    // 28
    };

    constexpr double zunik_c[120] = {
        1.00000000000000000e+00,  -2.08333333333333333e-01, 1.25000000000000000e-01,  3.34201388888888889e-01,  //   0
        -4.01041666666666667e-01, 7.03125000000000000e-02,  -1.02581259645061728e+00, 1.84646267361111111e+00,  //   4
        -8.91210937500000000e-01, 7.32421875000000000e-02,  4.66958442342624743e+00,  -1.12070026162229938e+01, //   8
        8.78912353515625000e+00,  -2.36408691406250000e+00, 1.12152099609375000e-01,  -2.82120725582002449e+01, //  12
        8.46362176746007346e+01,  -9.18182415432400174e+01, 4.25349987453884549e+01,  -7.36879435947963170e+00, //  16
        2.27108001708984375e-01,  2.12570130039217123e+02,  -7.65252468141181642e+02, 1.05999045252799988e+03,  //  20
        -6.99579627376132541e+02, 2.18190511744211590e+02,  -2.64914304869515555e+01, 5.72501420974731445e-01,  //  24
        -1.91945766231840700e+03, 8.06172218173730938e+03,  -1.35865500064341374e+04, 1.16553933368645332e+04,  //  28
        -5.30564697861340311e+03, 1.20090291321635246e+03,  -1.08090919788394656e+02, 1.72772750258445740e+00,  //  32
        2.02042913309661486e+04,  -9.69805983886375135e+04, 1.92547001232531532e+05,  -2.03400177280415534e+05, //  36
        1.22200464983017460e+05,  -4.11926549688975513e+04, 7.10951430248936372e+03,  -4.93915304773088012e+02, //  40
        6.07404200127348304e+00,  -2.42919187900551333e+05, 1.31176361466297720e+06,  -2.99801591853810675e+06, //  44
        3.76327129765640400e+06,  -2.81356322658653411e+06, 1.26836527332162478e+06,  -3.31645172484563578e+05, //  48
        4.52187689813627263e+04,  -2.49983048181120962e+03, 2.43805296995560639e+01,  3.28446985307203782e+06,  //  52
        -1.97068191184322269e+07, 5.09526024926646422e+07,  -7.41051482115326577e+07, 6.63445122747290267e+07,  //  56
        -3.75671766607633513e+07, 1.32887671664218183e+07,  -2.78561812808645469e+06, 3.08186404612662398e+05,  //  60
        -1.38860897537170405e+04, 1.10017140269246738e+02,  -4.93292536645099620e+07, 3.25573074185765749e+08,  //  64
        -9.39462359681578403e+08, 1.55359689957058006e+09,  -1.62108055210833708e+09, 1.10684281682301447e+09,  //  68
        -4.95889784275030309e+08, 1.42062907797533095e+08,  -2.44740627257387285e+07, 2.24376817792244943e+06,  //  72
        -8.40054336030240853e+04, 5.51335896122020586e+02,  8.14789096118312115e+08,  -5.86648149205184723e+09, //  76
        1.86882075092958249e+10,  -3.46320433881587779e+10, 4.12801855797539740e+10,  -3.30265997498007231e+10, //  80
        1.79542137311556001e+10,  -6.56329379261928433e+09, 1.55927986487925751e+09,  -2.25105661889415278e+08, //  84
        1.73951075539781645e+07,  -5.49842327572288687e+05, 3.03809051092238427e+03,  -1.46792612476956167e+10, //  88
        1.14498237732025810e+11,  -3.99096175224466498e+11, 8.19218669548577329e+11,  -1.09837515608122331e+12, //  92
        1.00815810686538209e+12,  -6.45364869245376503e+11, 2.87900649906150589e+11,  -8.78670721780232657e+10, //  96
        1.76347306068349694e+10,  -2.16716498322379509e+09, 1.43157876718888981e+08,  -3.87183344257261262e+06, // 100
        1.82577554742931747e+04,  2.86464035717679043e+11,  -2.40629790002850396e+12, 9.10934118523989896e+12,  // 104
        -2.05168994109344374e+13, 3.05651255199353206e+13,  -3.16670885847851584e+13, 2.33483640445818409e+13,  // 108
        -1.23204913055982872e+13, 4.61272578084913197e+12,  -1.19655288019618160e+12, 2.05914503232410016e+11,  // 112
        -2.18229277575292237e+10, 1.24700929351271032e+09,  -2.91883881222208134e+07, 1.18838426256783253e+05   // 116
    };

    constexpr double dgamln_gln[100] = {
        0.00000000000000000e+00, 0.00000000000000000e+00, 6.93147180559945309e-01, 1.79175946922805500e+00, //   0
        3.17805383034794562e+00, 4.78749174278204599e+00, 6.57925121201010100e+00, 8.52516136106541430e+00, //   4
        1.06046029027452502e+01, 1.28018274800814696e+01, 1.51044125730755153e+01, 1.75023078458738858e+01, //   8
        1.99872144956618861e+01, 2.25521638531234229e+01, 2.51912211827386815e+01, 2.78992713838408916e+01, //  12
        3.06718601060806728e+01, 3.35050734501368889e+01, 3.63954452080330536e+01, 3.93398841871994940e+01, //  16
        4.23356164607534850e+01, 4.53801388984769080e+01, 4.84711813518352239e+01, 5.16066755677643736e+01, //  20
        5.47847293981123192e+01, 5.80036052229805199e+01, 6.12617017610020020e+01, 6.45575386270063311e+01, //  24
        6.78897431371815350e+01, 7.12570389671680090e+01, 7.46582363488301644e+01, 7.80922235533153106e+01, //  28
        8.15579594561150372e+01, 8.50544670175815174e+01, 8.85808275421976788e+01, 9.21361756036870925e+01, //  32
        9.57196945421432025e+01, 9.93306124547874269e+01, 1.02968198614513813e+02, 1.06631760260643459e+02, //  36
        1.10320639714757395e+02, 1.14034211781461703e+02, 1.17771881399745072e+02, 1.21533081515438634e+02, //  40
        1.25317271149356895e+02, 1.29123933639127215e+02, 1.32952575035616310e+02, 1.36802722637326368e+02, //  44
        1.40673923648234259e+02, 1.44565743946344886e+02, 1.48477766951773032e+02, 1.52409592584497358e+02, //  48
        1.56360836303078785e+02, 1.60331128216630907e+02, 1.64320112263195181e+02, 1.68327445448427652e+02, //  52
        1.72352797139162802e+02, 1.76395848406997352e+02, 1.80456291417543771e+02, 1.84533828861449491e+02, //  56
        1.88628173423671591e+02, 1.92739047287844902e+02, 1.96866181672889994e+02, 2.01009316399281527e+02, //  60
        2.05168199482641199e+02, 2.09342586752536836e+02, 2.13532241494563261e+02, 2.17736934113954227e+02, //  64
        2.21956441819130334e+02, 2.26190548323727593e+02, 2.30439043565776952e+02, 2.34701723442818268e+02, //  68
        2.38978389561834323e+02, 2.43268849002982714e+02, 2.47572914096186884e+02, 2.51890402209723194e+02, //  72
        2.56221135550009525e+02, 2.60564940971863209e+02, 2.64921649798552801e+02, 2.69291097651019823e+02, //  76
        2.73673124285693704e+02, 2.78067573440366143e+02, 2.82474292687630396e+02, 2.86893133295426994e+02, //  80
        2.91323950094270308e+02, 2.95766601350760624e+02, 3.00220948647014132e+02, 3.04686856765668715e+02, //  84
        3.09164193580146922e+02, 3.13652829949879062e+02, 3.18152639620209327e+02, 3.22663499126726177e+02, //  88
        3.27185287703775217e+02, 3.31717887196928473e+02, 3.36261181979198477e+02, 3.40815058870799018e+02, //  92
        3.45379407062266854e+02, 3.49954118040770237e+02, 3.54539085519440809e+02, 3.59134205369575399e+02  //  96
    };

    constexpr double dgamln_cf[22] = {
        8.33333333333333333e-02, -2.77777777777777778e-03, 7.93650793650793651e-04, -5.95238095238095238e-04, //  0
        8.41750841750841751e-04, -1.91752691752691753e-03, 6.41025641025641026e-03, -2.95506535947712418e-02, //  4
        1.79644372368830573e-01, -1.39243221690590112e+00, 1.34028640441683920e+01, -1.56848284626002017e+02, //  8
        2.19310333333333333e+03, -3.61087712537249894e+04, 6.91472268851313067e+05, -1.52382215394074162e+07, // 12
        3.82900751391414141e+08, -1.08822660357843911e+10, 3.47320283765002252e+11, -1.23696021422692745e+13, // 16
        4.88788064793079335e+14, -2.13203339609193739e+16                                                     // 20
    };

    inline int acai(
        std::complex<double> z, double fnu, int kode, int mr, int n, std::complex<double> *y, double rl, double tol,
        double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZACAI
        //***REFER TO  ZAIRY
        //
        //     ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
        //
        //         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
        //                 MP=PI*MR*std::complex<double>(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.
        //
        //***ROUTINES CALLED  ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,AZABS
        //***END PROLOGUE  ZACAI

        std::complex<double> csgn, cspn, c1, c2, zn, cy[2];
        double arg, ascle, az, cpn, dfnu, fmr, sgn, spn, yy;
        int inu, iuf, nn, nw;
        double pi = 3.14159265358979324;
        int nz = 0;
        zn = -z;
        az = std::abs(z);
        nn = n;
        dfnu = fnu + (n - 1);
        if ((az > 2.0) && (az * az * 0.25 > dfnu + 1.0)) {
            /* 20 */
            if (az >= rl) {
                //
                // ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
                //
                nw = asyi(zn, fnu, kode, nn, y, rl, tol, elim, alim);
            } else {
                //
                // MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
                //
                nw = mlri(zn, fnu, kode, nn, y, tol);
            }
            if (nw < 0) {
                nz = -1;
                if (nw == -2) {
                    nz = -2;
                }
                return nz;
            }
        } else {
            //
            // POWER SERIES FOR THE I FUNCTION
            //
            seri(zn, fnu, kode, nn, y, tol, elim, alim);
        }
        /* 40 */
        //
        // ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
        //
        nw = bknu(zn, fnu, kode, 1, &cy[0], tol, elim, alim);
        if (nw != 0) {
            nz = -1;
            if (nw == -2) {
                nz = -2;
            }
            return nz;
        }
        fmr = mr;
        sgn = (fmr < 0.0 ? pi : -pi);
        csgn = std::complex<double>(0.0, sgn);
        if (kode != 1) {
            yy = -std::imag(zn);
            cpn = std::cos(yy);
            spn = std::sin(yy);
            csgn *= std::complex<double>(cpn, spn);
        }
        //
        // CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
        // WHEN FNU IS LARGE
        //
        inu = (int)fnu;
        arg = (fnu - inu) * sgn;
        cpn = std::cos(arg);
        spn = std::sin(arg);
        cspn = std::complex<double>(cpn, spn);
        if (inu % 2 == 1) {
            cspn = -cspn;
        }
        c1 = cy[0];
        c2 = y[0];
        if (kode != 1) {
            iuf = 0;
            ascle = 1e3 * d1mach[0] / tol;
            nw = s1s2(zn, &c1, &c2, ascle, alim, &iuf);
            nz += nw;
        }
        y[0] = cspn * c1 + csgn * c2;
        return nz;
    }

    inline int acon(
        std::complex<double> z, double fnu, int kode, int mr, int n, std::complex<double> *y, double rl, double fnul,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZACON
        //***REFER TO  ZBESK,ZBESH
        //
        //     ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
        //
        //         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
        //                 MP=PI*MR*std::complex<double>(0.0,1.0)
        //
        //     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
        //     HALF Z PLANE
        //
        //***ROUTINES CALLED  ZBINU,ZBKNU,ZS1S2,D1MACH,AZABS,ZMLT
        //***END PROLOGUE  ZACON

        std::complex<double> ck, cs, cscl, cscr, csgn, cspn, c1, c2, rz, sc1, sc2 = 0.0, st, s1, s2, zn;
        double arg, ascle, as2, bscle, c1i, c1m, c1r, fmr, sgn, yy;
        int i, inu, iuf, kflag, nn, nw, nz;
        double pi = 3.14159265358979324;
        std::complex<double> cy[2] = {0.0};
        std::complex<double> css[3] = {0.0};
        std::complex<double> csr[3] = {0.0};
        double bry[3] = {0.0};

        nz = 0;
        zn = -z;
        nn = n;
        nw = binu(zn, fnu, kode, nn, y, rl, fnul, tol, elim, alim);
        if (nw >= 0) {
            //
            // ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
            //
            nn = (n > 2 ? 2 : n);
            nw = bknu(zn, fnu, kode, nn, cy, tol, elim, alim);
            if (nw == 0) {
                s1 = cy[0];
                fmr = mr;
                sgn = (fmr < 0 ? pi : -pi);
                csgn = std::complex<double>(0.0, sgn);
                if (kode != 1) {
                    yy = -std::imag(zn);
                    csgn *= std::complex<double>(std::cos(yy), std::sin(yy));
                }
                inu = (int)fnu;
                arg = (fnu - inu) * sgn;
                cspn = std::complex<double>(std::cos(arg), std::sin(arg));
                if (inu % 2 == 1) {
                    cspn = -cspn;
                }
                iuf = 0;
                c1 = s1;
                c2 = y[0];
                ascle = 1e3 * d1mach[0] / tol;
                if (kode != 1) {
                    nw = s1s2(zn, &c1, &c2, ascle, alim, &iuf);
                    nz += nw;
                    sc1 = c1;
                }
                y[0] = cspn * c1 + csgn * c2;
                if (n == 1) {
                    return nz;
                }
                cspn = -cspn;
                s2 = cy[1];
                c1 = s2;
                c2 = y[1];
                if (kode != 1) {
                    nw = s1s2(zn, &c1, &c2, ascle, alim, &iuf);
                    nz += nw;
                    sc2 = c1;
                }
                y[1] = cspn * c1 + csgn * c2;
                if (n == 2) {
                    return nz;
                }
                cspn = -cspn;
                rz = 2.0 / zn;
                ck = (fnu + 1.0) * rz;
                //
                // SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
                //
                cscl = 1.0 / tol;
                cscr = tol;
                css[0] = cscl;
                css[1] = 1.0;
                css[2] = cscr;
                csr[0] = cscr;
                csr[1] = 1.0;
                csr[2] = cscl;
                bry[0] = ascle;
                bry[1] = 1.0 / ascle;
                bry[2] = d1mach[1];
                as2 = std::abs(s2);
                kflag = 2;
                if (as2 <= bry[0]) {
                    kflag = 1;
                } else {
                    if (as2 >= bry[1]) {
                        kflag = 3;
                    }
                }
                bscle = bry[kflag - 1];
                s1 *= css[kflag - 1];
                s2 *= css[kflag - 1];
                cs = csr[kflag - 1];
                for (i = 3; i < (n + 1); i++) {
                    st = s2;
                    s2 = ck * s2 + s1;
                    s1 = st;
                    c1 = s2 * cs;
                    st = c1;
                    c2 = y[i - 1];
                    if (kode != 1) {
                        if (iuf >= 0) {
                            nw = s1s2(zn, &c1, &c2, ascle, alim, &iuf);
                            nz += nw;
                            sc1 = sc2;
                            sc2 = c1;
                            if (iuf == 3) {
                                iuf = -4;
                                s1 = sc1 * css[kflag - 1];
                                s2 = sc2 * css[kflag - 1];
                                st = sc2;
                            }
                        }
                    }
                    y[i - 1] = cspn * c1 + csgn * c2;
                    ck += rz;
                    cspn = -cspn;
                    if (kflag < 3) {
                        c1r = std::fabs(std::real(c1));
                        c1i = std::fabs(std::imag(c1));
                        c1m = std::fmax(c1r, c1i);
                        if (c1m > bscle) {
                            kflag += 1;
                            bscle = bry[kflag - 1];
                            s1 *= cs;
                            s2 = st;
                            s1 *= css[kflag - 1];
                            s2 *= css[kflag - 1];
                            cs = csr[kflag - 1];
                        }
                    }
                }
                return nz;
            }
        }
        nz = -1;
        if (nw == -2) {
            nz = -2;
        }
        return nz;
    }

    inline std::complex<double> airy(std::complex<double> z, int id, int kode, int *nz, int *ierr) {

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

        std::complex<double> ai, csq, cy[1], s1, s2, trm1, trm2, zta, z3;
        double aa, ad, ak, alim, atrm, az, az3, bk, ck, dig, dk, d1, d2, elim, fid, fnu, rl, r1m5, sfac, tol, zi, zr,
            bb, alaz;
        int iflag, k, k1, k2, mr, nn;
        double tth = 2. / 3.;
        double c1 = 0.35502805388781723926;   /* 1/(Gamma(2/3) * 3**(2/3)) */
        double c2 = 0.25881940379280679841;   /* 1/(Gamma(1/3) * 3**(1/3)) */
        double coef = 0.18377629847393068317; /* 1 / (sqrt(3) * PI)        */

        *ierr = 0;
        *nz = 0;
        ai = 0.;
        if ((id < 0) || (id > 1)) {
            *ierr = 1;
        }
        if ((kode < 1) || (kode > 2)) {
            *ierr = 1;
        }
        if (*ierr != 0)
            return 0.;
        az = std::abs(z);
        tol = d1mach[3];
        fid = id;

        if (az <= 1.0) {
            //
            // POWER SERIES FOR ABS(Z) <= 1.
            //
            s1 = 1.0;
            s2 = 1.0;
            if (az < tol) {
                aa = 1e3 * d1mach[0];
                s1 = 0.;
                if (id != 1) {
                    if (az > aa) {
                        s1 = c2 * z;
                    }
                    ai = c1 - s1;
                    return ai;
                }
                ai = -c2;
                aa = std::sqrt(aa);
                if (az > aa) {
                    s1 = z * z * 0.5;
                }
                ai += s1 * c1;
                return ai;
            }
            aa = az * az;
            if (aa >= tol / az) {
                trm1 = 1.0;
                trm2 = 1.0;
                atrm = 1.0;
                z3 = z * z * z;
                az3 = az * aa;
                ak = 2.0 + fid;
                bk = 3.0 - fid - fid;
                ck = 4.0 - fid;
                dk = 3.0 + fid + fid;
                d1 = ak * dk;
                d2 = bk * ck;
                ad = (d1 > d2 ? d2 : d1);
                ak = 24.0 + 9.0 * fid;
                bk = 30.0 - 9.0 * fid;
                for (int k = 1; k < 26; k++) {
                    trm1 *= z3 / d1;
                    s1 += trm1;
                    trm2 *= z3 / d2;
                    s2 += trm2;
                    atrm *= az3 / ad;
                    d1 += ak;
                    d2 += bk;
                    ad = (d1 > d2 ? d2 : d1);
                    if (atrm < tol * ad) {
                        break;
                    }
                    ak += 18.0;
                    bk += 18.0;
                }
            }
            if (id != 1) {
                ai = s1 * c1 - z * s2 * c2;
                if (kode == 1) {
                    return ai;
                }
                zta = z * std::sqrt(z) * tth;
                ai *= std::exp(zta);
                return ai;
            }
            ai = -s2 * c2;
            if (az > tol) {
                ai += z * z * s1 * c1 / (1. + fid);
            }
            if (kode == 1) {
                return ai;
            }
            zta = z * std::sqrt(z) * tth;
            return ai * std::exp(zta);
        }
        //
        // CASE FOR ABS(Z) > 1.0
        //
        fnu = (1.0 + fid) / 3.0;
        //
        // SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        // TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        // ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        // EXP(-ELIM) < EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        // EXP(ELIM) > EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        // UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        // RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        // DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        //
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = (aa > 18.0 ? 18.0 : aa);
        aa *= 2.303;
        alim = elim + (-aa > -41.45 ? -aa : -41.45);
        rl = 1.2 * dig + 3.0;
        alaz = std::log(az);
        //
        // TEST FOR RANGE
        //
        aa = 0.5 / tol;
        bb = i1mach[8] * 0.5;
        aa = (aa > bb ? bb : aa);
        aa = std::pow(aa, tth);
        if (az > aa) {
            *ierr = 4;
            *nz = 0;
            return 0.;
        }
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        csq = std::sqrt(z);
        zta = z * csq * tth;
        //
        // RE(ZTA) <= 0 WHEN RE(Z) < 0, ESPECIALLY WHEN IM(Z) IS SMALL
        //
        iflag = 0;
        sfac = 1.0;
        zi = std::imag(z);
        zr = std::real(z);
        ak = std::imag(zta);
        if (zr < 0.0) {
            bk = std::real(zta);
            ck = -std::fabs(bk);
            zta = std::complex<double>(ck, ak);
        }
        if ((zi == 0.0) && (zr <= 0.0)) {
            zta = std::complex<double>(0.0, ak);
        }
        aa = std::real(zta);
        if ((aa < 0.0) || (zr <= 0.0)) {
            if (kode != 2) {
                //
                // OVERFLOW TEST
                //
                if (aa <= -alim) {
                    aa = -aa + 0.25 * alaz;
                    iflag = 1;
                    sfac = tol;
                    if (aa > elim) {
                        /* 270 */
                        *nz = 0;
                        *ierr = 2;
                        return ai;
                    }
                }
            }
            //
            // CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
            //
            mr = 1;
            if (zi < 0.0) {
                mr = -1;
            }
            nn = acai(zta, fnu, kode, mr, 1, &cy[0], rl, tol, elim, alim);
            if (nn < 0) {
                if (nn == -1) {
                    *nz = 1;
                    return 0.;
                } else {
                    *nz = 0;
                    *ierr = 5;
                    return 0.;
                }
            }
            *nz += nn;
        } else {
            if (kode != 2) {
                //
                // UNDERFLOW TEST
                //
                if (aa >= alim) {
                    aa = -aa - 0.25 * alaz;
                    iflag = 2;
                    sfac = 1.0 / tol;
                    if (aa < -elim) {
                        *nz = 1;
                        return 0.;
                    }
                }
            }
            *nz = bknu(zta, fnu, kode, 1, &cy[0], tol, elim, alim);
        }
        s1 = cy[0] * coef;

        if (iflag == 0) {
            if (id != 1) {
                return csq * s1;
            }
            return (-z * s1);
        }
        s1 *= sfac;
        if (id != 1) {
            s1 *= csq;
            return (s1 / sfac);
        }
        s1 *= -z;
        return (s1 / sfac);
    }

    inline int asyi(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, double rl, double tol,
        double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZASYI
        //***REFER TO  ZBESI,ZBESK
        //
        //     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=1.
        //
        //***ROUTINES CALLED  D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT
        //***END PROLOGUE  ZASYI

        std::complex<double> ak1, ck, cs1, cs2, cz, dk, ez, p1, rz, s2;
        double aa, acz, aez, ak, arg, arm, atol, az, bb, bk, dfnu;
        double dnu2, fdn, rtr1, s, sgn, sqk, x, yy;
        int ib, il, inu, j, jl, k, koded, m, nn;
        double pi = 3.14159265358979324;
        double rpi = 0.159154943091895336; /* (1 / pi) */
        int nz = 0;
        az = std::abs(z);
        x = std::real(z);
        arm = 1e3 * d1mach[0];
        rtr1 = std::sqrt(arm);
        il = (n > 2 ? 2 : n);
        dfnu = fnu + (n - il);
        // OVERFLOW TEST
        ak1 = std::sqrt(rpi / z);
        cz = z;
        if (kode == 2) {
            cz = std::complex<double>(0.0, std::imag(z));
        }
        acz = std::real(cz);
        if (std::fabs(acz) <= elim) {
            dnu2 = dfnu + dfnu;
            koded = 1;
            if (!((std::fabs(acz) > alim) && (n > 2))) {
                koded = 0;
                ak1 *= std::exp(cz);
            }
            fdn = 0.;
            if (dnu2 > rtr1) {
                fdn = dnu2 * dnu2;
            }
            ez = z * 8.;
            // WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE
            // RELATIVE TO THE FIRST RECIPROCAL POWER SINCE THIS
            // IS THE LEADING TERM OF THE EXPANSION FOR THE
            // IMAGINARY PART.
            aez = 8. * az;
            s = tol / aez;
            jl = (int)(rl + rl) + 2;
            yy = std::imag(z);
            p1 = 0.;
            if (yy != 0.) {
                inu = (int)fnu;
                arg = (fnu - inu) * pi;
                inu += n - il;
                ak = -std::sin(arg);
                bk = std::cos(arg);
                if (yy < 0.) {
                    bk = -bk;
                }
                p1 = std::complex<double>(ak, bk);
                if (inu % 2 == 1) {
                    p1 = -p1;
                }
            }
            for (int k = 1; k < (il + 1); k++) {
                sqk = fdn - 1.;
                atol = s * std::fabs(sqk);
                sgn = 1.;
                cs1 = 1.;
                cs2 = 1.;
                ck = 1.;
                ak = 0.;
                aa = 1.;
                bb = aez;
                dk = ez;
                j = 1;
                for (j = 1; j < (jl + 1); j++) {
                    ck *= sqk / dk;
                    cs2 += ck;
                    sgn = -sgn;
                    cs1 += ck * sgn;
                    dk += ez;
                    aa *= std::fabs(sqk) / bb;
                    bb += aez;
                    ak += 8.;
                    sqk -= ak;
                    if (aa <= atol) {
                        break;
                    }
                }
                if ((j == (jl + 1)) && (aa > atol)) {
                    return -2;
                }

                /* 50 */
                s2 = cs1;
                if (x + x < elim) {
                    s2 += p1 * cs2 * std::exp(-z - z);
                }
                fdn += 8. * dfnu + 4.;
                p1 = -p1;
                m = n - il + k;
                y[m - 1] = s2 * ak1;
            }
            if (n <= 2) {
                return nz;
            }
            nn = n;
            k = nn - 2;
            ak = k;
            rz = 2. / z;
            ib = 3;
            for (int i = ib; i < (nn + 1); i++) {
                y[k - 1] = (ak + fnu) * rz * y[k] + y[k + 1];
                ak -= 1.;
                k -= 1;
            }
            if (koded == 0) {
                return nz;
            }
            ck = std::exp(cz);
            for (int i = 0; i < nn; i++) {
                y[i] *= ck;
            }
            /* 90 */
            return nz;
        }
        /* 100 */
        return -1;
    }

    inline int besh(std::complex<double> z, double fnu, int kode, int m, int n, std::complex<double> *cy, int *ierr) {

        //***BEGIN PROLOGUE  ZBESH
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
        //             BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT
        //***DESCRIPTION
        //
        //                      ***A DOUBLE PRECISION ROUTINE***
        //         ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
        //         HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
        //         OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
        //         Z.NE.std::complex<double>(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
        //         ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS
        //
        //         CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z)       MM=3-2*M,   I**2=-1.
        //
        //         WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND
        //         LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE
        //         NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
        //
        //         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI), Z.NE.std::complex<double>(0.0D0,0.0D0),
        //                    -PT.LT.ARG(Z).LE.PI
        //           FNU    - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             CY(J)=H(M,FNU+J-1,Z),   J=1,...,N
        //                        = 2  RETURNS
        //                             CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
        //                                  J=1,...,N  ,  I**2=-1
        //           M      - KIND OF HANKEL FUNCTION, M=1 OR 2
        //           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
        //
        //         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
        //           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
        //                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
        //                    CY(J)=H(M,FNU+J-1,Z)  OR
        //                    CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))  J=1,...,N
        //                    DEPENDING ON KODE, I**2=-1.
        //           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
        //                    NZ= 0   , NORMAL RETURN
        //                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
        //                              TO UNDERFLOW, CY(J)=std::complex<double>(0.0D0,0.0D0)
        //                              J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
        //                              Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
        //                              HALF PLANES, NZ STATES ONLY THE NUMBER
        //                              OF UNDERFLOWS.
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU TOO
        //                            LARGE OR CABS(Z) TOO SMALL OR BOTH
        //                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
        //                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
        //                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
        //                            ACCURACY
        //                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
        //                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
        //                            CANCE BY ARGUMENT REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //
        //***LONG DESCRIPTION
        //
        //         THE COMPUTATION IS CARRIED OUT BY THE RELATION
        //
        //         H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP))
        //             MP=MM*HPI*I,  MM=3-2*M,  HPI=PI/2,  I**2=-1
        //
        //         FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE
        //         RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED
        //         TO THE LEFT HALF PLANE BY THE RELATION
        //
        //         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
        //         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
        //
        //         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
        //
        //         EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z
        //         PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2.  EXPONENTIAL
        //         GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES.  SCALING
        //         BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE
        //         WHOLE Z PLANE FOR Z TO INFINITY.
        //
        //         FOR NEGATIVE ORDERS,THE FORMULAE
        //
        //               H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I)
        //               H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I)
        //                         I**2=-1
        //
        //         CAN BE USED.
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
        //         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
        //         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
        //         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
        //         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
        //         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
        //         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
        //         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
        //         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
        //         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
        //         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
        //         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
        //         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
        //         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
        //         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
        //         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH
        //***END PROLOGUE  ZBESH

        std::complex<double> zn, zt, csgn;
        double aa, alim, aln, arg, az, cpn, dig, elim, fmm, fn, fnul, rhpi, rl, r1m5, sgn, spn, tol, ufl, xn, xx, yn,
            yy, bb, ascle, rtol, atol;
        int i, inu, inuh, ir, k, k1, k2, mm, mr, nn, nuf, nw, nz;

        double hpi = 1.57079632679489662; /* 0.5 PI */

        nz = 0;
        xx = std::real(z);
        yy = std::imag(z);
        *ierr = 0;

        if ((xx == 0.0) && (yy == 0.0)) {
            *ierr = 1;
        }
        if (fnu < 0.0) {
            *ierr = 1;
        }
        if ((m < 1) || (m > 2)) {
            *ierr = 1;
        }
        if ((kode < 1) || (kode > 2)) {
            *ierr = 1;
        }
        if (n < 1) {
            *ierr = 1;
        }
        if (*ierr != 0) {
            return nz;
        }
        nn = n;
        //
        //  SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        //  TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        //  ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        //  EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        //  EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        //  UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        //  RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        //  DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        //  FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
        //
        tol = std::fmax(d1mach[3], 1e-18);
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = std::fmin(aa, 18.0);
        aa *= 2.303;
        alim = elim + std::fmax(-aa, -41.45);
        fnul = 10.0 + 6.0 * (dig - 3.0);
        rl = 1.2 * dig + 3.0;
        fn = fnu + (nn - 1);
        mm = 3 - m - m;
        fmm = mm;
        zn = z * std::complex<double>(0.0, -fmm);
        xn = std::real(zn);
        yn = std::imag(zn);
        //
        // TEST FOR PROPER RANGE
        //
        az = std::abs(z);
        bb = d1mach[1] * 0.5;
        aa = std::fmin(0.5 / tol, bb);
        if ((az > aa) || (fn > aa)) {
            *ierr = 4;
            return 0;
        } /* GO TO 260 */
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        if (fn > aa) {
            *ierr = 3;
        }
        //
        // OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
        //
        ufl = d1mach[0] * 1.0e3;
        if (az < ufl) {
            *ierr = 2;
            return 0;
        } /* GO TO 230 */
        if (fnu <= fnul) {
            //
            // Untangling GOTOs with explicit conditions
            //
            if ((fn > 1.0) && (fn <= 2.0) && (az <= tol)) {
                /* Failed through all checks */
                arg = 0.5 * az;
                aln = -fn * std::log(arg);
                if (aln > elim) {
                    *ierr = 2;
                    return 0;
                } /* GO TO 230 */
                /* GO TO 70 */
            } else if ((fn > 1.0) && (fn <= 2.0) && (az > tol)) {
                /* Failed all but the az > tol hence do nothing and GO TO 70 */
            } else if ((fn > 1.0) && (fn > 2.0)) {
                /* GO TO 60 */
                nuf = uoik(zn, fnu, kode, 2, nn, cy, tol, elim, alim);
                if (nuf < 0) {
                    *ierr = 2;
                    return 0;
                } /* GO TO 230 */
                nz += nuf;
                nn -= nuf;
                //
                // HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
                // IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
                //
                if (nn == 0) {
                    /* GO TO 140 */
                    if (xn < 0.0) {
                        *ierr = 2;
                        return 0;
                    } /* GO TO 230 */
                    return nz;
                }
                /* GO TO 70 */
            } else {
                /* Passed the first hence GO TO 70 */
            }

            /* GO TO 70 */
            //
            // More GOTOs untangling
            //
            if ((xn < 0.0) || ((xn == 0.0) && (yn < 0.0) && (m == 2))) {
                /* GO TO 80 */
                mr = -mm;
                nw = acon(zn, fnu, kode, mr, nn, cy, rl, fnul, tol, elim, alim);
                if (nw < 0) {
                    /* GO TO 240 */
                    if (nw == -1) {
                        *ierr = 2;
                        return 0;
                    } /* GO TO 230 */
                    *ierr = 5;
                    return 0;
                }
                nz = nw;
                /* GO TO 110 */
            } else {
                //
                // RIGHT HALF PLANE COMPUTATION, XN >= 0.  .AND.  (XN.NE.0.  .OR.
                // YN >= 0.  .OR.  M=1)
                //
                nz = bknu(zn, fnu, kode, nn, cy, tol, elim, alim);
                /* GO TO 110 */
            }
        } else {
            /* GO TO 90 */
            //
            // UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU > FNUL
            //
            mr = 0;
            if (!((xn >= 0.0) && ((xn != 0.0) || (yn >= 0.0) || (m != 2)))) {
                mr = -mm;
                if ((xn == 0.0) && (yn < 0.0)) {
                    zn = -zn;
                }
            }
            /* GO TO 100 */
            nw = bunk(zn, fnu, kode, mr, nn, cy, tol, elim, alim);
            if (nw < 0) {
                /* GO TO 240 */
                if (nw == -1) {
                    *ierr = 2;
                    return 0;
                } /* GO TO 230 */
                *ierr = 5;
                return 0;
            }
            nz += nw;
        }
        /* 110 */
        //
        // H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
        // ZT=EXP(-FMM*HPI*I) = std::complex<double>(0.0,-FMM), FMM=3-2*M, M=1,2
        //
        sgn = (-fmm < 0 ? -hpi : hpi);
        //
        // CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
        // WHEN FNU IS LARGE
        //
        inu = (int)fnu;
        inuh = inu / 2;
        ir = inu - 2 * inuh;
        arg = (fnu - (inu - ir)) * sgn;
        rhpi = 1.0 / sgn;
        cpn = rhpi * std::cos(arg);
        spn = -rhpi * std::sin(arg);
        csgn = std::complex<double>(spn, cpn);
        if (inuh % 2 == 1) {
            csgn = -csgn;
        }
        zt = std::complex<double>(0.0, -fmm);
        rtol = 1.0 / tol;
        ascle = ufl * rtol;
        for (i = 1; i < (nn + 1); i++) {
            zn = cy[i - 1];
            atol = 1.0;
            if (std::fmax(std::fabs(std::real(zn)), std::fabs(std::imag(zn))) <= ascle) {
                zn *= rtol;
                atol = tol;
            }
            zn *= csgn;
            cy[i - 1] = zn * atol;
            csgn *= zt;
        }
        return nz;
    }

    inline int besi(std::complex<double> z, double fnu, int kode, int n, std::complex<double> *cy, int *ierr) {

        //***BEGIN PROLOGUE  ZBESI
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
        //             MODIFIED BESSEL FUNCTION OF THE FIRST KIND
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //***DESCRIPTION
        //
        //                    ***A DOUBLE PRECISION ROUTINE***
        //         ON KODE=1, ZBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
        //         BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE
        //         ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE
        //         -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESI RETURNS THE SCALED
        //         FUNCTIONS
        //
        //         CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z)   J = 1,...,N , X=REAL(Z)
        //
        //         WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND
        //         RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
        //         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
        //         (REF. 1).
        //
        //         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI),  -PI.LT.ARG(Z).LE.PI
        //           FNU    - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             CY(J)=I(FNU+J-1,Z), J=1,...,N
        //                        = 2  RETURNS
        //                             CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N
        //           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
        //
        //         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
        //           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
        //                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
        //                    CY(J)=I(FNU+J-1,Z)  OR
        //                    CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X))  J=1,...,N
        //                    DEPENDING ON KODE, X=REAL(Z)
        //           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
        //                    NZ= 0   , NORMAL RETURN
        //                    NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
        //                              TO UNDERFLOW, CY(J)=std::complex<double>(0.0D0,0.0D0)
        //                              J = N-NZ+1,...,N
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(Z) TOO
        //                            LARGE ON KODE=1
        //                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
        //                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
        //                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
        //                            ACCURACY
        //                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
        //                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
        //                            CANCE BY ARGUMENT REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //
        //***LONG DESCRIPTION
        //
        //         THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR
        //         SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z),
        //         THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A
        //         NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE
        //         UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z)
        //         FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE
        //         SEQUENCES OR REDUCE ORDERS WHEN NECESSARY.
        //
        //         THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND
        //         CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA
        //
        //         I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z)  REAL(Z).GT.0.0
        //                       M = +I OR -I,  I**2=-1
        //
        //         FOR NEGATIVE ORDERS,THE FORMULA
        //
        //              I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
        //
        //         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
        //         THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
        //         INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE
        //         NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
        //         K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
        //         TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
        //         UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
        //         OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
        //         LARGE MEANS FNU.GT.CABS(Z).
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
        //         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
        //         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
        //         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
        //         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
        //         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
        //         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
        //         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
        //         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
        //         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
        //         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
        //         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
        //         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
        //         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
        //         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
        //         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZBINU,I1MACH,D1MACH
        //***END PROLOGUE  ZBESI

        std::complex<double> csgn, zn;
        double aa, alim, arg, atol, ascle, az, bb, dig, elim, fn, fnul, rl, rtol, r1m5, tol, xx, yy;
        int i, inu, k, k1, k2, nn, nz;
        double pi = 3.14159265358979324;

        *ierr = 0;
        nz = 0;
        if (fnu < 0.0) {
            *ierr = 1;
        }
        if ((kode < 1) || (kode > 2)) {
            *ierr = 1;
        }
        if (n < 1) {
            *ierr = 1;
        }
        if (*ierr != 0) {
            return nz;
        }
        xx = std::real(z);
        yy = std::imag(z);
        //
        //  SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        //  TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        //  ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        //  EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        //  EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        //  UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        //  RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        //  DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        //  FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
        //
        tol = std::fmax(d1mach[3], 1e-18);
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = std::fmin(aa, 18.0);
        aa *= 2.303;
        alim = elim + std::fmax(-aa, -41.45);
        rl = 1.2 * dig + 3.0;
        fnul = 10.0 + 6.0 * (dig - 3.0);
        //
        // TEST FOR PROPER RANGE
        //
        az = std::abs(z);
        fn = fnu + (n - 1);
        aa = 0.5 / tol;
        bb = i1mach[8] * 0.5;
        aa = std::fmin(aa, bb);
        if ((az > aa) || (fn > aa)) {
            *ierr = 4;
            return 0;
        }
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        if (fn > aa) {
            *ierr = 3;
        }
        zn = z;
        csgn = 1.0;
        if (xx < 0.0) {
            zn = -z;
            //
            // CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
            // WHEN FNU IS LARGE
            //
            inu = (int)fnu;
            arg = (fnu - inu) * pi;
            if (yy < 0.0) {
                arg = -arg;
            }
            csgn = std::complex<double>(std::cos(arg), std::sin(arg));
            if (inu % 2 == 1) {
                csgn = -csgn;
            }
        }
        /* 40 */
        nz = binu(zn, fnu, kode, n, cy, rl, fnul, tol, elim, alim);
        if (nz < 0) {
            if (nz == -2) {
                *ierr = 5;
                return 0;
            }
            *ierr = 2;
            return 0;
        }
        if (xx >= 0.0) {
            return nz;
        }
        //
        // ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
        //
        nn = n - nz;
        if (nn == 0) {
            return nz;
        }
        rtol = 1.0 / tol;
        ascle = d1mach[0] * rtol * 1e3;
        for (i = 1; i < (nn + 1); i++) {
            zn = cy[i - 1];
            atol = 1.0;
            if (std::fmax(std::fabs(std::real(zn)), std::fabs(std::imag(zn))) <= ascle) {
                zn *= rtol;
                atol = tol;
            }
            cy[i - 1] = atol * (zn * csgn);
            csgn = -csgn;
        }
        *ierr = 0;
        return nz;
    }

    inline int besj(std::complex<double> z, double fnu, int kode, int n, std::complex<double> *cy, int *ierr) {

        //***BEGIN PROLOGUE  ZBESJ
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  J-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
        //             BESSEL FUNCTION OF FIRST KIND
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE THE J-BESSEL FUNCTION OF A COMPLEX ARGUMENT
        //***DESCRIPTION
        //
        //                      ***A DOUBLE PRECISION ROUTINE***
        //         ON KODE=1, CBESJ COMPUTES AN N MEMBER  SEQUENCE OF COMPLEX
        //         BESSEL FUNCTIONS CY(I)=J(FNU+I-1,Z) FOR REAL, NONNEGATIVE
        //         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
        //         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESJ RETURNS THE SCALED
        //         FUNCTIONS
        //
        //         CY(I)=EXP(-ABS(Y))*J(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
        //
        //         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
        //         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
        //         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
        //         (REF. 1).
        //
        //         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI),  -PI.LT.ARG(Z).LE.PI
        //           FNU    - ORDER OF INITIAL J FUNCTION, FNU.GE.0.0D0
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             CY(I)=J(FNU+I-1,Z), I=1,...,N
        //                        = 2  RETURNS
        //                             CY(I)=J(FNU+I-1,Z)EXP(-ABS(Y)), I=1,...,N
        //           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
        //
        //         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
        //           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
        //                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
        //                    CY(I)=J(FNU+I-1,Z)  OR
        //                    CY(I)=J(FNU+I-1,Z)EXP(-ABS(Y))  I=1,...,N
        //                    DEPENDING ON KODE, Y=AIMAG(Z).
        //           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
        //                    NZ= 0   , NORMAL RETURN
        //                    NZ.GT.0 , LAST NZ COMPONENTS OF CY SET  ZERO DUE
        //                              TO UNDERFLOW, CY(I)=std::complex<double>(0.0D0,0.0D0),
        //                              I = N-NZ+1,...,N
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, AIMAG(Z)
        //                            TOO LARGE ON KODE=1
        //                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
        //                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
        //                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
        //                            ACCURACY
        //                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
        //                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
        //                            CANCE BY ARGUMENT REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //
        //***LONG DESCRIPTION
        //
        //         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
        //
        //         J(FNU,Z)=EXP( FNU*PI*I/2)*I(FNU,-I*Z)    AIMAG(Z).GE.0.0
        //
        //         J(FNU,Z)=EXP(-FNU*PI*I/2)*I(FNU, I*Z)    AIMAG(Z).LT.0.0
        //
        //         WHERE I**2 = -1 AND I(FNU,Z) IS THE I BESSEL FUNCTION.
        //
        //         FOR NEGATIVE ORDERS,THE FORMULA
        //
        //              J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU)
        //
        //         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
        //         THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
        //         INTEGER,THE MAGNITUDE OF J(-FNU,Z)=J(FNU,Z)*COS(PI*FNU) IS A
        //         LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
        //         Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
        //         TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
        //         UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
        //         OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
        //         LARGE MEANS FNU.GT.CABS(Z).
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
        //         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
        //         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
        //         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
        //         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
        //         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
        //         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
        //         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
        //         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
        //         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
        //         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
        //         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
        //         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
        //         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
        //         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
        //         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZBINU,I1MACH,D1MACH
        //***END PROLOGUE  ZBESJ

        std::complex<double> ci, csgn, zn;
        double aa, alim, arg, dig, elim, fnul, rl, r1, r1m5, r2, tol, yy, az, fn, bb, ascle, rtol, atol;
        int i, inu, inuh, ir, k1, k2, nl, nz, k;
        double hpi = 1.570796326794896619;

        *ierr = 0;
        nz = 0;
        if (fnu < 0.0)
            *ierr = 1;
        if (kode < 1 || kode > 2)
            *ierr = 1;
        if (n < 1)
            *ierr = 1;
        if (*ierr != 0)
            return nz;
        //
        // SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        // TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        // ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        // EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        // EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        // UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        // RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        // DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        // FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
        //
        tol = std::fmax(d1mach[3], 1e-18);
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = std::fmin(aa, 18.0);
        aa *= 2.303;
        alim = elim + std::fmax(-aa, -41.45);
        fnul = 10.0 + 6.0 * (dig - 3.0);
        rl = 1.2 * dig + 3.0;
        //
        // TEST FOR PROPER RANGE
        //
        yy = std::imag(z);
        az = std::abs(z);
        fn = fnu + (n - 1);

        aa = 0.5 / tol;
        bb = d1mach[1] * 0.5;
        aa = std::fmin(aa, bb);
        if ((az > aa) || (fn > aa)) {
            *ierr = 4;
            return 0;
        }
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        if (fn > aa) {
            *ierr = 3;
        }
        //
        // CALCULATE CSGN = EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
        // WHEN FNU IS LARGE
        //
        ci.imag(1);
        inu = (int)fnu;
        inuh = inu / 2;
        ir = inu - 2 * inuh;
        arg = (fnu - (inu - ir)) * hpi;
        r1 = std::cos(arg);
        r2 = std::sin(arg);
        csgn = std::complex<double>(r1, r2);
        if (inuh % 2 == 1) {
            csgn = -csgn;
        }
        //
        // ZN IS IN THE RIGHT HALF PLANE
        //
        zn = -z * ci;
        if (yy < 0.0) {
            zn = -zn;
            csgn = std::conj(csgn);
            ci = std::conj(ci);
        }
        nz = binu(zn, fnu, kode, n, cy, rl, fnul, tol, elim, alim);
        if (nz < 0) {
            if (nz == -2) {
                *ierr = 5;
                return 0;
            }
            *ierr = 2;
            return 0;
        }
        nl = n - nz;
        if (nl == 0) {
            return nz;
        }
        rtol = 1.0 / tol;
        ascle = d1mach[0] * rtol * 1e3;
        for (i = 1; i < (nl + 1); i++) {
            zn = cy[i - 1];
            aa = std::real(zn);
            bb = std::imag(zn);
            atol = 1.0;
            if (std::fmax(std::fabs(aa), std::fabs(bb)) <= ascle) {
                zn *= rtol;
                atol = tol;
            }
            cy[i - 1] = atol * (zn * csgn);
            csgn = csgn * ci;
        }
        return nz;
    }

    inline int besk(std::complex<double> z, double fnu, int kode, int n, std::complex<double> *cy, int *ierr) {

        //***BEGIN PROLOGUE  ZBESK
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
        //             MODIFIED BESSEL FUNCTION OF THE SECOND KIND,
        //             BESSEL FUNCTION OF THE THIRD KIND
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //***DESCRIPTION
        //
        //                      ***A DOUBLE PRECISION ROUTINE***
        //
        //         ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
        //         BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE
        //         ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.std::complex<double>(0.0,0.0)
        //         IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK
        //         RETURNS THE SCALED K FUNCTIONS,
        //
        //         CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N,
        //
        //         WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND
        //         RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND
        //         NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
        //         FUNCTIONS (REF. 1).
        //
        //         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI), Z.NE.std::complex<double>(0.0D0,0.0D0),
        //                    -PI.LT.ARG(Z).LE.PI
        //           FNU    - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0
        //           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             CY(I)=K(FNU+I-1,Z), I=1,...,N
        //                        = 2  RETURNS
        //                             CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
        //
        //         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
        //           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
        //                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
        //                    CY(I)=K(FNU+I-1,Z), I=1,...,N OR
        //                    CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
        //                    DEPENDING ON KODE
        //           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW.
        //                    NZ= 0   , NORMAL RETURN
        //                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
        //                              TO UNDERFLOW, CY(I)=std::complex<double>(0.0D0,0.0D0),
        //                              I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0
        //                              NZ STATES ONLY THE NUMBER OF UNDERFLOWS
        //                              IN THE SEQUENCE.
        //
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
        //                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
        //                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
        //                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
        //                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
        //                            ACCURACY
        //                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
        //                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
        //                            CANCE BY ARGUMENT REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //
        //***LONG DESCRIPTION
        //
        //         EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS
        //         DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD
        //         RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT
        //         HALF PLANE BY THE RELATION
        //
        //         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
        //         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
        //
        //         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
        //
        //         FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED
        //         BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS.
        //
        //         FOR NEGATIVE ORDERS, THE FORMULA
        //
        //                       K(-FNU,Z) = K(FNU,Z)
        //
        //         CAN BE USED.
        //
        //         CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS
        //         AVAILABLE.
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
        //         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
        //         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
        //         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
        //         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
        //         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
        //         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
        //         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
        //         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
        //         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
        //         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
        //         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
        //         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
        //         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
        //         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
        //         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983.
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH
        //***END PROLOGUE  ZBESK

        double xx = std::real(z);
        double yy = std::imag(z);
        double aa, alim, aln, arg, az, dig, elim, fn, fnul, rl, r1m5, tol, ufl, bb;
        int k, k1, k2, mr, nn, nuf, nw, nz;

        *ierr = 0;
        nz = 0;

        if ((yy == 0.0) && (xx == 0.0)) {
            *ierr = 1;
        }
        if (fnu < 0.0) {
            *ierr = 1;
        }
        if (kode < 1 || kode > 2) {
            *ierr = 1;
        }
        if (n < 1) {
            *ierr = 1;
        }
        if (*ierr != 0) {
            return nz;
        }

        nn = n;
        //
        //  SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        //  TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        //  ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        //  EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        //  EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        //  UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        //  RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        //  DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        //  FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
        //
        tol = std::fmax(d1mach[3], 1e-18);
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = std::fmin(aa, 18.0);
        aa *= 2.303;
        alim = elim + std::fmax(-aa, -41.45);
        fnul = 10.0 + 6.0 * (dig - 3.0);
        rl = 1.2 * dig + 3.0;
        //
        // TEST FOR PROPER RANGE
        //
        az = std::abs(z);
        fn = fnu + (nn - 1);
        aa = 0.5 / tol;
        bb = i1mach[8] * 0.5;
        aa = std::fmin(aa, bb);
        if ((az > aa) || (fn > aa)) {
            *ierr = 4;
            return 0;
        }
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        if (fn > aa) {
            *ierr = 3;
        }
        //
        // OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
        //
        ufl = d1mach[0] * 1.0E+3;
        if (az < ufl) {
            *ierr = 2;
            return 0;
        }
        if (fnu <= fnul) {
            if (fn > 1.0) {
                if (fn <= 2.0) {
                    if (az <= tol) {
                        arg = 0.5 * az;
                        aln = -fn * std::log(arg);
                        if (aln > elim) {
                            *ierr = 2;
                            return 0;
                        }
                    }
                    /* GO TO 60 */
                } else {
                    nuf = uoik(z, fnu, kode, 2, nn, cy, tol, elim, alim);
                    if (nuf < 0) {
                        *ierr = 2;
                        return 0;
                    }
                    nz += nuf;
                    nn -= nuf;
                    //
                    // HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
                    // IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
                    //
                    if (nn == 0) {
                        if (xx < 0.0) {
                            *ierr = 2;
                            return 0;
                        }
                        return nz;
                    }
                }
            }

            /* 60 */
            if (xx >= 0.0) {
                //
                // RIGHT HALF PLANE COMPUTATION, REAL(Z) >= 0.
                //
                nw = bknu(z, fnu, kode, nn, cy, tol, elim, alim);
                if (nw < 0) {
                    if (nw == -1) {
                        *ierr = 2;
                    } else {
                        *ierr = 5;
                    }
                    return 0;
                }
                return nw;
            }
            /* 70 */
            //
            // LEFT HALF PLANE COMPUTATION
            // PI/2 < ARG(Z) <= PI AND -PI < ARG(Z) < -PI/2.
            //
            if (nz != 0) {
                *ierr = 2;
                return 0;
            }
            mr = 1;
            if (yy < 0.0) {
                mr = -1;
            }
            nw = acon(z, fnu, kode, mr, nn, cy, rl, fnul, tol, elim, alim);
            if (nw < 0) {
                if (nw == -1) {
                    *ierr = 2;
                } else {
                    *ierr = 5;
                }
                return 0;
            }
            return nw;
        }

        /* 80 */
        //
        // UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU > FNUL
        //
        mr = 0;
        if (xx < 0.0) {
            mr = 1;
            if (yy < 0.0) {
                mr = -1;
            }
        }
        nw = bunk(z, fnu, kode, mr, nn, cy, tol, elim, alim);
        if (nw < 0) {
            if (nw == -1) {
                *ierr = 2;
            } else {
                *ierr = 5;
            }
            return 0;
        }
        nz += nw;
        return nz;
    }

    inline int besy(std::complex<double> z, double fnu, int kode, int n, std::complex<double> *cy, int *ierr) {

        //***BEGIN PROLOGUE  ZBESY
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
        //             BESSEL FUNCTION OF SECOND KIND
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
        //***DESCRIPTION
        //
        //                      ***A DOUBLE PRECISION ROUTINE***
        //
        //         ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
        //         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
        //         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
        //         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
        //         FUNCTIONS
        //
        //         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
        //
        //         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
        //         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
        //         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
        //         (REF. 1).
        //
        //         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI), Z.NE.std::complex<double>(0.0D0,0.0D0),
        //                    -PI.LT.ARG(Z).LE.PI
        //           FNU    - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
        //                        = 2  RETURNS
        //                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
        //                             WHERE Y=AIMAG(Z)
        //           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
        //           CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
        //           CWRKI    AT LEAST N
        //
        //         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
        //           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
        //                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
        //                    CY(I)=Y(FNU+I-1,Z)  OR
        //                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
        //                    DEPENDING ON KODE.
        //           NZ     - NZ=0 , A NORMAL RETURN
        //                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
        //                    UNDERFLOW (GENERALLY ON KODE=2)
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
        //                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
        //                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
        //                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
        //                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
        //                            ACCURACY
        //                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
        //                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
        //                            CANCE BY ARGUMENT REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //                    IERR=6, Memory allocation failed.
        //
        //***LONG DESCRIPTION
        //
        //         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
        //
        //         Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
        //
        //         WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
        //         AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
        //
        //         FOR NEGATIVE ORDERS,THE FORMULA
        //
        //              Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
        //
        //         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
        //         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
        //         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
        //         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
        //         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
        //         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
        //         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
        //         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
        //         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
        //         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
        //         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
        //         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
        //         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
        //         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
        //         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
        //         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
        //         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
        //         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
        //         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
        //         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
        //         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
        //         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
        //         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
        //         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZBESH,I1MACH,D1MACH
        //***END PROLOGUE  ZBESY

        std::complex<double> c1, c2, hci, st;
        double elim, exr, exi, ey, tay, xx, yy, ascle, rtol, atol, tol, aa, bb, r1m5;
        int i, k, k1, k2, nz, nz1, nz2;

        xx = std::real(z);
        yy = std::imag(z);
        *ierr = 0;
        nz = 0;

        if ((xx == 0.0) && (yy == 0.0)) {
            *ierr = 1;
        }
        if (fnu < 0.0) {
            *ierr = 1;
        }
        if ((kode < 1) || (kode > 2)) {
            *ierr = 1;
        }
        if (n < 1) {
            *ierr = 1;
        }
        if (*ierr != 0) {
            return nz;
        }

        hci = std::complex<double>(0.0, 0.5);
        nz1 = besh(z, fnu, kode, 1, n, cy, ierr);
        if ((*ierr != 0) && (*ierr != 3)) {
            return 0;
        }

        auto cwrk = std::unique_ptr<std::complex<double>[]>{new (std::nothrow) std::complex<double>[n]};
        if (cwrk == nullptr) {
            *ierr = 6; // Memory allocation failed.
            return 0;
        }

        nz2 = besh(z, fnu, kode, 2, n, cwrk.get(), ierr);
        if ((*ierr != 0) && (*ierr != 3)) {
            return 0;
        }

        nz = (nz1 > nz2 ? nz2 : nz1);
        if (kode != 2) {
            for (i = 1; i < (n + 1); i++) {
                cy[i - 1] = hci * (cwrk[i - 1] - cy[i - 1]);
            }
            return nz;
        }

        tol = std::fmax(d1mach[3], 1e-18);
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        //
        // ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
        //
        elim = 2.303 * (k * r1m5 - 3.0);
        exr = std::cos(xx);
        exi = std::sin(xx);
        ey = 0.0;
        tay = std::fabs(yy + yy);
        if (tay < elim) {
            ey = std::exp(-tay);
        }
        if (yy < 0.0) {
            /* 90 */
            c1 = std::complex<double>(exr, exi);
            c2 = ey * std::complex<double>(exr, -exi);
        } else {
            c1 = ey * std::complex<double>(exr, exi);
            c2 = std::complex<double>(exr, -exi);
        }

        nz = 0;
        rtol = 1.0 / tol;
        ascle = 1e3 * d1mach[0] * rtol;
        for (i = 1; i < (n + 1); i++) {
            aa = std::real(cwrk[i - 1]);
            bb = std::imag(cwrk[i - 1]);
            atol = 1.0;
            if (std::fmax(std::fabs(aa), std::fabs(bb)) <= ascle) {
                aa *= rtol;
                bb *= rtol;
                atol = tol;
            }

            st = std::complex<double>(aa, bb) * c2 * atol;
            aa = std::real(cy[i - 1]);
            bb = std::imag(cy[i - 1]);
            atol = 1.0;
            if (std::fmax(std::fabs(aa), std::fabs(bb)) <= ascle) {
                aa *= rtol;
                bb *= rtol;
                atol = tol;
            }

            st -= std::complex<double>(aa, bb) * c1 * atol;
            cy[i - 1] = st * hci;
            if ((st == 0.0) && (ey == 0.0)) {
                nz += 1;
            }
        }

        return nz;
    }

    inline int binu(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *cy, double rl, double fnul,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZBINU
        //***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY
        //
        //     ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
        //
        //***ROUTINES CALLED  AZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
        //***END PROLOGUE  ZBINU

        std::complex<double> cw[2] = {0.};
        double az, dfnu;
        int inw, nlast, nn, nui, nw, nz;

        nz = 0;
        az = std::abs(z);
        nn = n;
        dfnu = fnu + n - 1;
        if ((az <= 2.) || (az * az * 0.25 <= (dfnu + 1.0))) {
            /* GOTO 10 */
            nw = seri(z, fnu, kode, nn, cy, tol, elim, alim);
            inw = std::abs(nw);
            nz += inw;
            nn -= inw;
            if (nn == 0) {
                return nz;
            }
            if (nw >= 0) {
                return nz;
            }
            dfnu = fnu + nn - 1;
        }
        /* GOTO 30 conditions*/
        //
        // ASYMPTOTIC EXPANSION FOR LARGE Z
        //
        if (az < rl) {
            /* 40 */
            if (dfnu <= 1.0) {
                /* 70 */
                //
                // MILLER ALGORITHM NORMALIZED BY THE SERIES
                //
                nw = mlri(z, fnu, kode, nn, cy, tol);
                if (nw < 0) {
                    nz = -1;
                    if (nw == -2) {
                        nz = -2;
                    }
                    return nz;
                }
                return nz;
            }
            /* GO TO 50 */
        } else {
            if ((dfnu <= 1.0) || (az + az >= dfnu * dfnu)) {
                /* 30 */
                //
                // ASYMPTOTIC EXPANSION FOR LARGE Z
                //
                nw = asyi(z, fnu, kode, nn, cy, rl, tol, elim, alim);
                if (nw < 0) {
                    nz = -1;
                    if (nw == -2) {
                        nz = -2;
                    }
                    return nz;
                }
                return nz;
            }
            /* GO TO 50 */
        }
        /* 50 */
        //
        // OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
        //
        nw = uoik(z, fnu, kode, 1, nn, cy, tol, elim, alim);
        if (nw < 0) {
            nz = -1;
            if (nw == -2) {
                nz = -2;
            }
            return nz;
        }
        nz += nw;
        nn -= nw;
        if (nn == 0) {
            return nz;
        }
        dfnu = fnu + (nn - 1);
        /* GOTO 110s handled here */
        if ((dfnu > fnul) || (az > fnul)) {
            //
            // INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
            //
            nui = (int)(fnul - dfnu) + 1;
            nui = (nui > 0 ? nui : 0);
            nw = buni(z, fnu, kode, nn, cy, nui, &nlast, fnul, tol, elim, alim);
            if (nw < 0) {
                nz = -1;
                if (nw == -2) {
                    nz = -2;
                }
                return nz;
            }
            nz += nw;
            if (nlast == 0) {
                return nz;
            }
            nn = nlast;
        }
        /* 60 */
        if (az <= rl) {
            /* 70 */
            nw = mlri(z, fnu, kode, nn, cy, tol);
            if (nw < 0) {
                nz = -1;
                if (nw == -2) {
                    nz = -2;
                }
                return nz;
            }
            return nz;
        }
        /* 80 */
        //
        // MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
        //
        //
        // OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
        //
        nw = uoik(z, fnu, kode, 2, 2, cw, tol, elim, alim);
        if (nw < 0) {
            nz = nn;
            /* 90 */
            for (int i = 0; i < nn; i++) {
                cy[i] = 0.0;
            }
            return nz;
        }
        /* 100 */
        if (nw > 0) {
            return -1;
        }
        nw = wrsk(z, fnu, kode, nn, cy, cw, tol, elim, alim);
        if (nw < 0) {
            nz = -1;
            if (nw == -2) {
                nz = -2;
            }
            return nz;
        }
        return nz;
    }

    inline std::complex<double> biry(std::complex<double> z, int id, int kode, int *ierr) {

        //***BEGIN PROLOGUE  ZBIRY
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  890801   (YYMMDD)
        //***CATEGORY NO.  B5K
        //***KEYWORDS  AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE AIRY FUNCTIONS BI(Z) AND DBI(Z) FOR COMPLEX Z
        //***DESCRIPTION
        //
        //                      ***A DOUBLE PRECISION ROUTINE***
        //         ON KODE=1, CBIRY COMPUTES THE COMPLEX AIRY FUNCTION BI(Z) OR
        //         ITS DERIVATIVE DBI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
        //         KODE=2, A SCALING OPTION CEXP(-AXZTA)*BI(Z) OR CEXP(-AXZTA)*
        //         DBI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN
        //         BOTH THE LEFT AND RIGHT HALF PLANES WHERE
        //         ZTA=(2/3)*Z*CSQRT(Z)=std::complex<double>(XZTA,YZTA) AND AXZTA=ABS(XZTA).
        //         DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
        //         MATHEMATICAL FUNCTIONS (REF. 1).
        //
        //         INPUT      ZR,ZI ARE DOUBLE PRECISION
        //           ZR,ZI  - Z=std::complex<double>(ZR,ZI)
        //           ID     - ORDER OF DERIVATIVE, ID=0 OR ID=1
        //           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
        //                    KODE= 1  RETURNS
        //                             BI=BI(Z)                 ON ID=0 OR
        //                             BI=DBI(Z)/DZ             ON ID=1
        //                        = 2  RETURNS
        //                             BI=CEXP(-AXZTA)*BI(Z)     ON ID=0 OR
        //                             BI=CEXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE
        //                             ZTA=(2/3)*Z*CSQRT(Z)=std::complex<double>(XZTA,YZTA)
        //                             AND AXZTA=ABS(XZTA)
        //
        //         OUTPUT     BIR,BII ARE DOUBLE PRECISION
        //           BIR,BII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
        //                    KODE
        //           IERR   - ERROR FLAG
        //                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
        //                    IERR=1, INPUT ERROR   - NO COMPUTATION
        //                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(Z)
        //                            TOO LARGE ON KODE=1
        //                    IERR=3, CABS(Z) LARGE      - COMPUTATION COMPLETED
        //                            LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
        //                            PRODUCE LESS THAN HALF OF MACHINE ACCURACY
        //                    IERR=4, CABS(Z) TOO LARGE  - NO COMPUTATION
        //                            COMPLETE LOSS OF ACCURACY BY ARGUMENT
        //                            REDUCTION
        //                    IERR=5, ERROR              - NO COMPUTATION,
        //                            ALGORITHM TERMINATION CONDITION NOT MET
        //
        //***LONG DESCRIPTION
        //
        //         BI AND DBI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE I BESSEL
        //         FUNCTIONS BY
        //
        //                BI(Z)=C*SQRT(Z)*( I(-1/3,ZTA) + I(1/3,ZTA) )
        //               DBI(Z)=C *  Z  * ( I(-2/3,ZTA) + I(2/3,ZTA) )
        //                               C=1.0/SQRT(3.0)
        //                             ZTA=(2/3)*Z**(3/2)
        //
        //         WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
        //
        //         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
        //         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
        //         OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
        //         THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
        //         THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
        //         FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
        //         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
        //         ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
        //         ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
        //         FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
        //         LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
        //         MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
        //         AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
        //         PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
        //         PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
        //         ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
        //         NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
        //         DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
        //         EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
        //         NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
        //         PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
        //         MACHINES.
        //
        //         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
        //         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
        //         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
        //         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
        //         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
        //         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
        //         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
        //         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
        //         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
        //         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
        //         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
        //         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
        //         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
        //         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
        //         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
        //         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
        //         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
        //         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
        //         OR -PI/2+P.
        //
        //***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
        //                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
        //                 COMMERCE, 1955.
        //
        //               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
        //
        //               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
        //                 1018, MAY, 1985
        //
        //               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
        //                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
        //                 MATH. SOFTWARE, 1986
        //
        //***ROUTINES CALLED  ZBINU,AZABS,ZDIV,AZSQRT,D1MACH,I1MACH
        //***END PROLOGUE  ZBIRY

        std::complex<double> bi, csq, s1, s2, trm1, trm2, zta, z3;
        double aa, ad, ak, alim, atrm, az, az3, bb, bk, ck, dig, dk, d1, d2, elim, fid, fmr, fnu, fnul, rl, r1m5, sfac,
            tol, zi, zr;
        int k, k1, k2, nz;
        std::complex<double> cy[2] = {0.0};
        double tth = 2. / 3.;
        double c1 = 0.614926627446000735150922369;   /* 1/( 3**(1/6) Gamma(2/3)) */
        double c2 = 0.448288357353826357914823710;   /* 3**(1/6) / Gamma(1/3)    */
        double coef = 0.577350269189625764509148780; /* sqrt( 1 / 3)             */
        double pi = 3.141592653589793238462643383;

        *ierr = 0;
        nz = 0;
        if ((id < 0) || (id > 1)) {
            *ierr = 1;
        }
        if ((kode < 1) || (kode > 2)) {
            *ierr = 1;
        }
        if (*ierr != 0) {
            return 0.0;
        }
        az = std::abs(z);
        tol = std::fmax(d1mach[3], 1e-18);
        fid = id;
        if (az <= 1.0) {
            //
            // POWER SERIES FOR ABS(Z) <= 1.
            //
            s1 = 1.0;
            s2 = 1.0;
            if (az < tol) {
                aa = c1 * (1.0 - fid) + fid * c2;
                return aa;
            }
            aa = az * az;
            if (aa >= tol / az) {
                trm1 = 1.0;
                trm2 = 1.0;
                atrm = 1.0;
                z3 = z * z * z;
                az3 = az * aa;
                ak = 2.0 + fid;
                bk = 3.0 - fid - fid;
                ck = 4.0 - fid;
                dk = 3.0 + fid + fid;
                d1 = ak * dk;
                d2 = bk * ck;
                ad = std::fmin(d1, d2);
                ak = 24.0 + 9.0 * fid;
                bk = 30.0 - 9.0 * fid;
                for (k = 1; k < 26; k++) {
                    trm1 *= z3 / d1;
                    s1 += trm1;
                    trm2 *= z3 / d2;
                    s2 += trm2;
                    atrm *= az3 / ad;
                    d1 += ak;
                    d2 += bk;
                    ad = std::fmin(d1, d2);
                    if (atrm < tol * ad) {
                        break;
                    }
                    ak += 18.0;
                    bk += 18.0;
                }
                /* 30 */
            }
            /* 40 */
            if (id != 1) {
                bi = s1 * c1 + z * s2 * c2;
                if (kode == 1) {
                    return bi;
                }
                zta = z * std::sqrt(z) * tth;
                aa = -std::fabs(std::real(zta));
                bi *= std::exp(aa);
                return bi;
            }
            /* 50 */
            bi = s2 * c2;
            if (az > tol) {
                bi += z * z * s1 * c1 / (1.0 + fid);
            }
            if (kode == 1) {
                return bi;
            }
            zta = z * std::sqrt(z) * tth;
            aa = -std::fabs(std::real(zta));
            bi *= std::exp(aa);
            return bi;
        }
        /* 70 */
        //
        // CASE FOR ABS(Z) > 1.0
        //
        fnu = (1.0 + fid) / 3.0;
        //
        // SET PARAMETERS RELATED TO MACHINE CONSTANTS.
        // TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
        // ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
        // EXP(-ELIM) < EXP(-ALIM)=EXP(-ELIM)/TOL    AND
        // EXP(ELIM) > EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
        // UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
        // RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
        // DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
        // FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
        //
        k1 = i1mach[14];
        k2 = i1mach[15];
        r1m5 = d1mach[4];
        k = (std::abs(k1) > std::abs(k2) ? std::abs(k2) : std::abs(k1));
        elim = 2.303 * (k * r1m5 - 3.0);
        k1 = i1mach[13] - 1;
        aa = r1m5 * k1;
        dig = std::fmin(aa, 18.0);
        aa *= 2.303;
        alim = elim + std::fmax(-aa, -41.45);
        rl = 1.2 * dig + 3.0;
        fnul = 10.0 + 6.0 * (dig - 3.0);
        //
        // TEST FOR RANGE
        //
        aa = 0.5 / tol;
        bb = i1mach[8] * 0.5;
        aa = std::fmin(aa, bb);
        aa = std::pow(aa, tth);
        if (az > aa) {
            *ierr = 4;
            return 0.0;
        }
        aa = std::sqrt(aa);
        if (az > aa) {
            *ierr = 3;
        }
        csq = std::sqrt(z);
        zta = z * csq * tth;
        //
        // RE(ZTA) <= 0 WHEN RE(Z) < 0, ESPECIALLY WHEN IM(Z) IS SMALL
        //
        sfac = 1.0;
        zi = std::imag(z);
        zr = std::real(z);
        ak = std::imag(zta);
        if (zr < 0.0) {
            bk = std::real(zta);
            ck = -std::fabs(bk);
            zta = std::complex<double>(ck, ak);
        }
        /* 80 */
        if ((zi == 0.0) && (zr <= 0.0)) {
            zta = std::complex<double>(0.0, ak);
        }
        /* 90 */
        aa = std::real(zta);
        if (kode != 2) {
            //
            // OVERFLOW TEST
            //
            bb = std::fabs(aa);
            if (bb >= alim) {
                bb += 0.25 * std::log(az);
                sfac = tol;
                if (bb > elim) {
                    *ierr = 2;
                    return 0.0;
                }
            }
        }
        /* 100 */
        fmr = 0.0;
        if ((aa < 0.0) || (zr <= 0.0)) {
            fmr = pi;
            if (zi < 0.0) {
                fmr = -pi;
            }
            zta = -zta;
        }
        /* 110 */
        //
        // AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA)
        // KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBINU
        //
        nz = binu(zta, fnu, kode, 1, cy, rl, fnul, tol, elim, alim);
        if (nz < 0) {
            if (nz == -1) {
                *ierr = 2;
            } else {
                *ierr = 5;
            }
            return 0.0;
        }
        aa = fmr * fnu;
        z3 = sfac;
        s1 = cy[0] * std::complex<double>(std::cos(aa), std::sin(aa)) * z3;
        fnu = (2.0 - fid) / 3.0;
        nz = binu(zta, fnu, kode, 2, cy, rl, fnul, tol, elim, alim);
        cy[0] *= z3;
        cy[1] *= z3;
        //
        // BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3
        //
        s2 = cy[0] * (fnu + fnu) / zta + cy[1];
        aa = fmr * (fnu - 1.0);
        s1 = (s1 + s2 * std::complex<double>(std::cos(aa), std::sin(aa))) * coef;
        if (id != 1) {
            s1 *= csq;
            bi = s1 / sfac;
            return bi;
        }
        /* 120 */
        s1 *= z;
        bi = s1 / sfac;
        return bi;
    }

    inline int bknu(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, double tol, double elim,
        double alim
    ) {

        //***BEGIN PROLOGUE  ZBKNU
        //***REFER TO  ZBESI,ZBESK,ZAIRY,ZBESH
        //
        //     ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
        //
        //***ROUTINES CALLED  DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,AZABS,ZDIV,
        //                    AZEXP,AZLOG,ZMLT,AZSQRT
        //***END PROLOGUE  ZBKNU

        std::complex<double> cch, ck, coef, crsc, cs, cscl, csh, cz, f, fmu, p, pt, p1, p2, q, rz, smu, st, s1, s2, zd;
        double aa, ak, ascle, a1, a2, bb, bk, caz, dnu, dnu2, etest, fc, fhs, fk, fks, g1, g2, p2i, p2m, p2r, rk, s, tm,
            t1, t2, xx, yy, elm, xd, yd, alas, as;
        int iflag, inu, k, kflag, kk, koded, j, ic, inub, i = 1;
        std::complex<double> cy[2];

        int kmax = 30;
        double r1 = 2.;
        double pi = 3.14159265358979324;
        double rthpi = 1.25331413731550025;
        double spi = 1.90985931710274403;
        double fpi = 1.89769999331517738;
        double tth = 2. / 3.;
        double cc[8] = {5.77215664901532861e-01, -4.20026350340952355e-02, -4.21977345555443367e-02,
                        7.21894324666309954e-03, -2.15241674114950973e-04, -2.01348547807882387e-05,
                        1.13302723198169588e-06, 6.11609510448141582e-09};
        xx = std::real(z);
        yy = std::imag(z);
        caz = std::abs(z);
        cscl = 1. / tol;
        crsc = tol;
        std::complex<double> css[3] = {cscl, 1., crsc};
        std::complex<double> csr[3] = {crsc, 1., cscl};
        double bry[3] = {1e3 * d1mach[0] / tol, tol / (1e3 * d1mach[0]), d1mach[1]};
        int nz = 0;
        iflag = 0;
        koded = kode;
        rz = 2. / z;
        inu = (int)(fnu + 0.5);
        dnu = fnu - inu;
        // Definitions for silencing initialization warnings.
        s1 = 0.0;
        s2 = 0.0;
        ck = 0.0;
        dnu2 = 0.0;
        if (std::fabs(dnu) != 0.5) {
            if (std::fabs(dnu) > tol) {
                dnu2 = dnu * dnu;
            }
            if (caz <= r1) {
                //
                //    SERIES FOR ABS(Z) <= R1
                //
                fc = 1.;
                smu = std::log(rz);
                fmu = smu * dnu;
                csh = std::sinh(fmu);
                cch = std::cosh(fmu);
                if (dnu != 0.0) {
                    fc = dnu * pi;
                    fc *= 1. / std::sin(fc);
                    smu = csh / dnu;
                }
                a2 = 1. + dnu;
                //
                // GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
                //
                t2 = std::exp(-gamln(a2));
                t1 = 1. / (t2 * fc);
                if (std::fabs(dnu) <= 0.1) {
                    //
                    // SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
                    //
                    ak = 1.;
                    s = cc[0];
                    for (int k = 2; k < 9; k++) {
                        ak *= dnu2;
                        tm = cc[k - 1] * ak;
                        s += tm;
                        if (std::fabs(tm) < tol) {
                            break;
                        }
                    }
                    g1 = -s;
                } else {
                    g1 = (t1 - t2) / (dnu + dnu);
                }
                g2 = 0.5 * (t1 + t2);
                f = fc * (g1 * cch + smu * g2);
                pt = std::exp(fmu);
                p = (0.5 / t2) * pt;
                q = (0.5 / t1) / pt;
                s1 = f;
                s2 = p;
                ak = 1.0;
                a1 = 1.0;
                ck = 1.0;
                bk = 1.0 - dnu2;
                if ((inu <= 0) && (n <= 1)) {
                    //
                    // GENERATE K(FNU,Z), 0.0D0  <=  FNU  <  0.5D0 AND N=1
                    //
                    if (caz >= tol) {
                        cz = z * z * 0.25;
                        t1 = 0.25 * caz * caz;
                        do {
                            f = (f * ak + p + q) / bk;
                            p = p / (ak - dnu);
                            q = q / (ak + dnu);
                            rk = 1.0 / ak;
                            ck *= cz * rk;
                            s1 += ck * f;
                            a1 *= t1 * rk;
                            bk += ak + ak + 1.0;
                            ak += 1.0;
                        } while (a1 > tol);
                    }
                    y[0] = s1;
                    if (koded == 1) {
                        return nz;
                    }
                    y[0] = s1 * std::exp(z);
                    return nz;
                }
                //
                // GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
                //
                if (caz >= tol) {
                    cz = z * z * 0.25;
                    t1 = 0.25 * caz * caz;
                    do {
                        f = (f * ak + p + q) / bk;
                        p *= 1.0 / (ak - dnu);
                        q *= 1.0 / (ak + dnu);
                        rk = 1. / ak;
                        ck *= cz * rk;
                        s1 += ck * f;
                        s2 += ck * (p - f * ak);
                        a1 *= t1 * rk;
                        bk += ak + ak + 1.0;
                        ak += 1.0;
                    } while (a1 > tol);
                }
                kflag = 2;
                bk = std::real(smu);
                a1 = fnu + 1.;
                ak = a1 * std::fabs(bk);
                if (ak > alim) {
                    kflag = 3;
                }
                p2 = s2 * css[kflag - 1];
                s2 = p2 * rz;
                s1 *= css[kflag - 1];
                if (koded != 1) {
                    f = std::exp(z);
                    s1 *= f;
                    s2 *= f;
                }
                goto L100;
            }
        }
        //
        // IFLAG=0 MEANS NO UNDERFLOW OCCURRED
        // IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
        // KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
        // RECURSION
        //
        coef = rthpi / std::sqrt(z);
        kflag = 2;
        if (koded != 2) {
            if (xx > alim) {
                koded = 2;
                iflag = 1;
                kflag = 2;
            } else {
                a1 = std::exp(-xx) * std::real(css[kflag - 1]);
                pt = a1 * std::complex<double>(std::cos(yy), -std::sin(yy));
                coef *= pt;
            }
        }

        if (std::fabs(dnu) == 0.5) {
            s1 = coef;
            s2 = coef;
            goto L100;
        }
        //
        //    MILLER ALGORITHM FOR ABS(Z) > R1
        //
        ak = std::fabs(std::cos(pi * dnu));
        if (ak == 0.) {
            s1 = coef;
            s2 = coef;
            goto L100;
        }
        fhs = std::fabs(0.25 - dnu2);
        if (fhs == 0.) {
            s1 = coef;
            s2 = coef;
            goto L100;
        }
        //
        // COMPUTE R2=F(E). IF ABS(Z) >= R2, USE FORWARD RECURRENCE TO
        // DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
        // 12 <= E <= 60. E IS COMPUTED FROM 2**(-E)=B**(1-DIGITS(0.0_dp))=
        // TOL WHERE B IS THE BASE OF THE ARITHMETIC.
        //
        t1 = (i1mach[13] - 1) * d1mach[4] * (std::log(10) / std::log(2));
        t1 = std::fmin(std::fmax(t1, 12.0), 60.0);
        t2 = tth * t1 - 6.0;
        t1 = std::fabs(std::atan2(yy, xx));
        if (t2 <= caz) {
            //
            // FORWARD RECURRENCE LOOP WHEN ABS(Z) >= R2
            //
            etest = ak / (pi * caz * tol);
            fk = 1.0;
            if (etest < 1.0) {
                goto L80;
            }
            fks = 2.0;
            rk = caz + caz + 2.0;
            a1 = 0.0;
            a2 = 1.0;
            for (i = 1; i < (kmax + 1); i++) {
                ak = fhs / fks;
                bk = rk / (fk + 1.0);
                tm = a2;
                a2 = bk * a2 - ak * a1;
                a1 = tm;
                rk += 2.;
                fks += fk + fk + 2.0;
                fhs += fk + fk;
                fk += 1.0;
                tm = std::fabs(a2) * fk;
                if (etest < tm) {
                    /* goto 160 */
                    break;
                }
                if (i == kmax) {
                    /* Didn't break so goes to 310 */
                    return -2;
                }
            }

            /* 160 */
            fk += spi * t1 * std::sqrt(t2 / caz);
            fhs = std::fabs(0.25 - dnu2);
        } else {
            //
            // COMPUTE BACKWARD INDEX K FOR ABS(Z) < R2
            //
            a2 = std::sqrt(caz);
            ak *= fpi / (tol * std::sqrt(a2));
            aa = 3.0 * t1 / (1.0 + caz);
            bb = 14.7 * t1 / (28.0 + caz);
            ak = (std::log(ak) + caz * std::cos(aa) / (1.0 + 0.008 * caz)) / std::cos(bb);
            fk = 0.12125 * ak * ak / caz + 1.5;
        }
    L80:
        //
        // BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
        //
        k = (int)fk;
        fk = (double)k;
        fks = fk * fk;
        p1 = 0.0;
        p2 = tol;
        cs = p2;
        for (i = 1; i < (k + 1); i++) {
            a1 = fks - fk;
            a2 = (fks + fk) / (a1 + fhs);
            rk = 2.0 / (fk + 1.);
            t1 = (fk + xx) * rk;
            t2 = yy * rk;
            pt = p2;
            p2 = (p2 * std::complex<double>(t1, t2) - p1) * a2;
            p1 = pt;
            cs += p2;
            fks = a1 - fk + 1.0;
            fk -= 1.0;
        }

        //
        // COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER SCALING
        //
        tm = std::abs(cs);
        pt = 1.0 / tm;
        s1 = pt * p2;
        cs = std::conj(cs) * pt;
        s1 *= coef * cs;
        if ((inu <= 0) && (n <= 1)) {
            zd = z;
            if (iflag == 1) {
                goto L190;
            }
            goto L130;
        }
        //
        // COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING
        //
        tm = std::abs(p2);
        pt = 1.0 / tm;
        p1 = pt * p1;
        p2 = std::conj(p2) * pt;
        pt = p1 * p2;
        s2 = s1 * (1. + (dnu + 0.5 - pt) / z);
        //
        // FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
        // SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
        //
    L100:
        ck = (dnu + 1.) * rz;
        if (n == 1) {
            inu -= 1;
        }
        if (inu <= 0) {
            if (n <= 1) {
                s1 = s2;
            }
            zd = z;
            if (iflag == 1) {
                goto L190;
            }
            goto L130;
        }
        inub = 1;
        if (iflag == 1) {
            goto L160;
        }
    L110:
        p1 = csr[kflag - 1];
        ascle = bry[kflag - 1];
        for (i = inub; i < inu + 1; i++) {
            st = s2;
            s2 = ck * s2 + s1;
            s1 = st;
            ck += rz;
            if (kflag < 3) {
                p2 = s2 * p1;
                p2m = std::fmax(std::fabs(std::real(p2)), std::fabs(std::imag(p2)));
                if (p2m > ascle) {
                    kflag += 1;
                    ascle = bry[kflag - 1];
                    s1 *= p1;
                    s2 = p2;
                    s1 *= css[kflag - 1];
                    s2 *= css[kflag - 1];
                    p1 = csr[kflag - 1];
                }
            }
        }
        if (n == 1) {
            s1 = s2;
        }

    L130:
        y[0] = s1 * csr[kflag - 1];
        if (n == 1) {
            return nz;
        }
        y[1] = s2 * csr[kflag - 1];
        if (n == 2) {
            return nz;
        }
        kk = 2;
    L140:
        kk += 1;
        if (kk > n) {
            return nz;
        }
        p1 = csr[kflag - 1];
        ascle = bry[kflag - 1];
        for (i = kk; i < (n + 1); i++) {
            p2 = s2;
            s2 = ck * s2 + s1;
            s1 = p2;
            ck += rz;
            p2 = s2 * p1;
            y[i - 1] = p2;
            if (kflag < 3) {
                p2m = std::fmax(std::fabs(std::real(p2)), std::fabs(std::imag(p2)));
                if (p2m > ascle) {
                    kflag += 1;
                    ascle = bry[kflag - 1];
                    s1 *= p1;
                    s2 = p2;
                    s1 *= css[kflag - 1];
                    s2 *= css[kflag - 1];
                    p1 = csr[kflag - 1];
                }
            }
        }
        return nz;
    //
    // IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
    //
    L160:
        elm = std::exp(-elim);
        ascle = bry[0];
        zd = z;
        xd = xx;
        yd = yy;
        ic = -1;
        j = 2;
        for (i = 1; i < (inu + 1); i++) {
            st = s2;
            s2 = ck * s2 + s1;
            s1 = st;
            ck += rz;
            as = std::abs(s2);
            alas = std::log(as);
            p2r = alas - xd;
            if (p2r >= -elim) {
                p2 = -zd + std::log(s2);
                p2r = std::real(p2);
                p2i = std::imag(p2);
                p2m = std::exp(p2r) / tol;
                p1 = p2m * std::complex<double>(std::cos(p2i), std::sin(p2i));
                if (!(uchk(p1, ascle, tol))) {
                    j = 3 - j;
                    cy[j - 1] = p1;
                    if (ic == i - 1) {
                        goto L180;
                    }
                    ic = i;
                    continue;
                }
            }
            if (alas >= 0.5 * elim) {
                xd -= elim;
                s1 *= elm;
                s2 *= elm;
                zd = std::complex<double>(xd, yd);
            }
        }
        if (n == 1) {
            s1 = s2;
        }
        goto L190;
    L180:
        kflag = 1;
        inub = i + 1;
        s2 = cy[j - 1];
        j = 3 - j;
        s1 = cy[j - 1];
        if (inub <= inu) {
            goto L110;
        }
        if (n == 1) {
            s1 = s2;
        }
        goto L130;
    L190:
        y[0] = s1;
        if (n != 1) {
            y[1] = s2;
        }
        ascle = bry[0];
        nz = kscl(zd, fnu, n, &y[0], rz, &ascle, tol, elim);
        inu = n - nz;
        if (inu <= 0) {
            return nz;
        }
        kk = nz + 1;
        s1 = y[kk - 1];
        y[kk - 1] = s1 * csr[0];
        if (inu == 1) {
            return nz;
        }
        kk = nz + 2;
        s2 = y[kk - 1];
        y[kk - 1] = s2 * csr[0];
        if (inu == 2) {
            return nz;
        }
        t2 = fnu + (kk - 1);
        ck = t2 * rz;
        kflag = 1;
        goto L140;
    }

    inline int buni(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, int nui, int *nlast, double fnul,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZBUNI
        //***REFER TO  ZBESI,ZBESK
        //
        //     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
        //
        //***ROUTINES CALLED  ZUNI1,ZUNI2,AZABS,D1MACH
        //***END PROLOGUE  ZBUNI

        std::complex<double> cscl, cscr, rz, st, s1, s2;
        double ax, ay, dfnu, fnui, gnu, xx, yy, ascle, str, sti, stm;
        int i, iflag, iform, k, nl, nw, nz;
        std::complex<double> cy[2] = {0.0};
        double bry[3] = {0.0};

        nz = 0;
        xx = std::real(z);
        yy = std::imag(z);
        ax = std::fabs(xx) * std::sqrt(3.);
        ay = std::fabs(yy);
        iform = 1;
        if (ay > ax) {
            iform = 2;
        }
        if (nui == 0) {
            if (iform != 2) {
                uni1(z, fnu, kode, n, y, &nw, nlast, fnul, tol, elim, alim);
            } else {
                uni2(z, fnu, kode, n, y, &nw, nlast, fnul, tol, elim, alim);
            }
            if (nw < 0) {
                nz = -1;
                if (nw == -2) {
                    nz = -2;
                }
                return nz;
            }
            return nw;
        }

        fnui = nui;
        dfnu = fnu + (n - 1);
        gnu = dfnu + fnui;
        if (iform != 2) {
            //
            // ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
            // -PI/3 <= ARG(Z) <= PI/3
            //
            uni1(z, gnu, kode, 2, cy, &nw, nlast, fnul, tol, elim, alim);
        } else {
            uni2(z, gnu, kode, 2, cy, &nw, nlast, fnul, tol, elim, alim);
        }
        if (nw >= 0) {
            if (nw != 0) {
                *nlast = n;
                return nz;
            }
            ay = std::abs(cy[0]);
            //
            // SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
            //
            bry[0] = 1e3 * d1mach[0] / tol;
            bry[1] = tol / 1e3 * d1mach[0];
            bry[2] = bry[1];
            iflag = 2;
            ascle = bry[1];
            ax = 1.0;
            cscl = ax;
            if (ay <= bry[0]) {
                iflag = 1;
                ascle = bry[0];
                ax = 1.0 / tol;
                cscl = ax;
            } else {
                if (ay >= bry[1]) {
                    iflag = 3;
                    ascle = bry[2];
                    ax = tol;
                    cscl = ax;
                }
            }
            ay = 1.0 / ax;
            cscr = ay;
            s1 = cy[1] * cscl;
            s2 = cy[0] * cscl;
            rz = 2.0 / z;
            for (i = 1; i < (nui + 1); i++) {
                st = s2;
                s2 = (dfnu + fnui) * rz * st + s1;
                s1 = st;
                fnui -= 1.0;
                if (iflag < 3) {
                    st = s2 * cscr;
                    str = std::fabs(std::real(st));
                    sti = std::fabs(std::imag(st));
                    stm = std::fmax(str, sti);
                    if (stm > ascle) {
                        iflag += 1;
                        ascle = bry[iflag - 1];
                        s1 *= cscr;
                        s2 = st;
                        ax *= tol;
                        ay = 1.0 / ax;
                        cscl = ax;
                        cscr = ay;
                        s1 *= cscl;
                        s2 *= cscl;
                    }
                }
            }
            y[n - 1] = s2 * cscr;
            if (n == 1) {
                return nz;
            }
            nl = n - 1;
            fnui = nl;
            k = nl;
            for (i = 0; i < nl; i++) {
                st = s2;
                s2 = (fnu + fnui) * rz * s2 + s1;
                s1 = st;
                st = s2 * cscr;
                y[k - 1] = st;
                fnui -= 1.0;
                k -= 1;
                if (iflag < 3) {
                    st = s2 * cscr;
                    str = std::fabs(std::real(st));
                    sti = std::fabs(std::imag(st));
                    stm = std::fmax(str, sti);
                    if (stm > ascle) {
                        iflag += 1;
                        ascle = bry[iflag - 1];
                        s1 *= cscr;
                        s2 = st;
                        ax *= tol;
                        ay = 1.0 / ax;
                        cscl = ax;
                        cscr = ay;
                        s1 *= cscl;
                        s2 *= cscl;
                    }
                }
            }
            return nz;
        }
        nz = -1;
        if (nw == -2) {
            nz = -2;
        }
        return nz;
    }

    inline int bunk(
        std::complex<double> z, double fnu, int kode, int mr, int n, std::complex<double> *y, double tol, double elim,
        double alim
    ) {

        //***BEGIN PROLOGUE  ZBUNK
        //***REFER TO  ZBESK,ZBESH
        //
        //     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
        //
        //***ROUTINES CALLED  ZUNK1,ZUNK2
        //***END PROLOGUE  ZBUNK

        double ax, ay;

        int nz = 0;
        ax = std::fabs(std::real(z)) * 1.7321;
        ay = std::fabs(std::imag(z));

        if (ay <= ax) {
            //
            // Asymptotic expansion for K(FNU,Z) for large FNU applied in
            // -PI/3 <= ARG(Z) <= PI/3
            //
            nz = unk1(z, fnu, kode, mr, n, y, tol, elim, alim);
        } else {
            //
            // Asymptotic expansion for H(2, FNU, Z*EXP(M*HPI)) for large FNU
            // applied in PI/3 < ABS(ARG(Z)) <= PI/2 where M = +I or -I and HPI = PI/2
            //
            nz = unk2(z, fnu, kode, mr, n, y, tol, elim, alim);
        }
        return nz;
    }

    inline double gamln(double z) {

        //***BEGIN PROLOGUE  DGAMLN
        //***DATE WRITTEN   830501   (YYMMDD)
        //***REVISION DATE  830501   (YYMMDD)
        //***CATEGORY NO.  B5F
        //***KEYWORDS  GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION
        //***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
        //***PURPOSE  TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION
        //***DESCRIPTION
        //
        //               **** A DOUBLE PRECISION ROUTINE ****
        //         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.
        //
        //     DESCRIPTION OF ARGUMENTS
        //
        //         INPUT      Z IS D0UBLE PRECISION
        //           Z      - ARGUMENT, Z.GT.0.0D0
        //
        //         OUTPUT      DGAMLN IS DOUBLE PRECISION
        //           DGAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0
        //           IERR    - ERROR FLAG
        //                     IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
        //                     IERR=1, Z.LE.0.0D0,    NO COMPUTATION
        //
        //
        //***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
        //                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
        //***ROUTINES CALLED  I1MACH,D1MACH
        //***END PROLOGUE  DGAMLN

        int i1m, mz;
        double fln, fz, rln, s, tlg, trm, tst, t1, wdtol, zdmy, zinc, zm, zmin, zp, zsq;
        const double con = 1.83787706640934548; /* LN(2*PI) */
        int nz = 0;
        if (z > 0.0) {
            if (z <= 101.0) {
                nz = (int)z;
                fz = z - nz;
                if (fz <= 0.0) {
                    if (nz <= 100) {
                        return dgamln_gln[nz - 1];
                    }
                }
            }
            wdtol = std::fmax(d1mach[3], 1e-18);
            i1m = i1mach[13];
            rln = d1mach[4] * i1m;
            fln = std::fmax(std::fmin(rln, 20.), 3.0) - 3.0;
            zm = 1.8 + 0.3875 * fln;
            mz = ((int)zm) + 1;
            zmin = mz;
            zdmy = z;
            zinc = 0.0;
            if (z < zmin) {
                zinc = zmin - nz;
                zdmy = z + zinc;
            }
            zp = 1. / zdmy;
            t1 = dgamln_cf[0] * zp;
            s = t1;
            if (zp >= wdtol) {
                zsq = zp * zp;
                tst = t1 * wdtol;
                for (int i = 2; i < 23; i++) {
                    zp *= zsq;
                    trm = dgamln_cf[i - 1] * zp;
                    if (std::fabs(trm) < tst) {
                        break;
                    }
                    s += trm;
                }
            }

            if (zinc == 0.) {
                tlg = std::log(z);
                return z * (tlg - 1.0) + 0.5 * (con - tlg) + s;
            }
            zp = 1.0;
            nz = (int)zinc;
            for (int i = 0; i < nz; i++) {
                zp *= (z + i);
            }
            tlg = std::log(zdmy);
            return zdmy * (tlg - 1.0) - std::log(zp) + 0.5 * (con - tlg) + s;
        }
        // Zero or negative argument
        return NAN;
    }

    inline int mlri(std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, double tol) {

        //***BEGIN PROLOGUE  ZMLRI
        //***REFER TO  ZBESI,ZBESK
        //
        //     ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
        //     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
        //
        //***ROUTINES CALLED  DGAMLN,D1MACH,AZABS,AZEXP,AZLOG,ZMLT
        //***END PROLOGUE  ZMLRI

        std::complex<double> ck, cnorm, pt, p1, p2, rz, sum;
        double ack, ak, ap, at, az, bk, fkap, fkk, flam, fnf, rho, rho2, scle, tfnf, tst, x;
        int i, iaz, ifnu, inu, itime, k, kk, km, m, nz;
        scle = d1mach[0] / tol;
        nz = 0;
        az = std::abs(z);
        x = std::real(z);
        iaz = (int)az;
        ifnu = (int)fnu;
        inu = ifnu + n - 1;
        at = iaz + 1;
        ck = at / z;
        rz = 2. / z;
        p1 = 0.;
        p2 = 1.;
        ack = (at + 1.0) / az;
        rho = ack + std::sqrt(ack * ack - 1.);
        rho2 = rho * rho;
        tst = (rho2 + rho2) / ((rho2 - 1.0) * (rho - 1.0));
        tst /= tol;
        //
        // COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
        //
        ak = at;
        i = 1;
        for (i = 1; i < 81; i++) {
            pt = p2;
            p2 = p1 - ck * p2;
            p1 = pt;
            ck += rz;
            ap = std::abs(p2);
            if (ap > tst * ak * ak) {
                break;
            }
            ak += 1.0;
            if (i == 80) {
                /* Exhausted loop without break */
                return -2;
            }
        }
        i += 1;
        k = 0;
        if (inu >= iaz) {
            //
            // COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
            //
            p1 = 0.0;
            p2 = 1.0;
            at = inu + 1;
            ck = at / z;
            ack = at / az;
            tst = std::sqrt(ack / tol);
            itime = 1;
            k = 1;
            for (k = 1; k < 81; k++) {
                pt = p2;
                p2 = p1 - ck * p2;
                p1 = pt;
                ck += rz;
                ap = std::abs(p2);
                if (ap >= tst) {
                    if (itime == 2) {
                        break;
                    }
                    ack = std::abs(ck);
                    flam = ack + std::sqrt(ack * ack - 1.0);
                    fkap = ap / std::abs(p1);
                    rho = std::fmin(flam, fkap);
                    tst *= std::sqrt(rho / (rho * rho - 1.0));
                    itime = 2;
                }
                if (k == 80) {
                    /* Exhausted loop without break */
                    return -2;
                }
            }
        }
        //
        // BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
        //
        k += 1;
        kk = std::fmax(i + iaz, k + inu);
        fkk = kk;
        p1 = 0.0;
        //
        // SCALE P2 AND SUM BY SCLE
        //
        p2 = scle;
        fnf = fnu - ifnu;
        tfnf = fnf + fnf;
        bk = gamln(fkk + tfnf + 1.0) - gamln(fkk + 1.0) - gamln(tfnf + 1.0);
        bk = std::exp(bk);
        sum = 0.;
        km = kk - inu;
        for (i = 1; i < (km + 1); i++) {
            pt = p2;
            p2 = p1 + (fkk + fnf) * rz * p2;
            p1 = pt;
            ak = 1. - tfnf / (fkk + tfnf);
            ack = bk * ak;
            sum += (ack + bk) * p1;
            bk = ack;
            fkk -= 1.;
        }
        y[n - 1] = p2;
        if (n != 1) {
            for (i = 2; i < (n + 1); i++) {
                pt = p2;
                p2 = p1 + (fkk + fnf) * rz * p2;
                p1 = pt;
                ak = 1. - tfnf / (fkk + tfnf);
                ack = bk * ak;
                sum += (ack + bk) * p1;
                bk = ack;
                fkk -= 1.;
                m = n - i + 1;
                y[m - 1] = p2;
            }
        }
        if (ifnu > 0) {
            for (i = 1; i < (ifnu + 1); i++) {
                pt = p2;
                p2 = p1 + (fkk + fnf) * rz * p2;
                p1 = pt;
                ak = 1. - tfnf / (fkk + tfnf);
                ack = bk * ak;
                sum += (ack + bk) * p1;
                bk = ack;
                fkk -= 1.;
            }
        }
        pt = z;
        if (kode == 2) {
            pt -= x;
        }
        p1 = -fnf * std::log(rz) + pt;
        ap = gamln(1. + fnf);
        pt = p1 - ap;
        //
        // THE DIVISION EXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
        // IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
        //
        p2 += sum;
        ap = std::abs(p2);
        p1 = 1. / ap;
        ck = std::exp(pt) * p1;
        pt = std::conj(p2) * p1;
        cnorm = ck * pt;
        for (int i = 0; i < n; i++) {
            y[i] *= cnorm;
        }
        return nz;
    }

    inline int kscl(
        std::complex<double> zr, double fnu, int n, std::complex<double> *y, std::complex<double> rz, double *ascle,
        double tol, double elim
    ) {

        //***BEGIN PROLOGUE  ZKSCL
        //***REFER TO  ZBESK
        //
        //     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
        //     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
        //     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
        //
        //***ROUTINES CALLED  ZUCHK,AZABS,AZLOG
        //***END PROLOGUE  ZKSCL

        std::complex<double> cy[2] = {0.};
        double as, acs, alas, fn, zri, xx;
        std::complex<double> s1, s2, cs, ck, zd;
        int nz = 0;
        int ic = 0;
        int nn = (n > 2 ? 2 : n);
        int kk = 0;
        int i;
        double elm = std::exp(-elim);
        xx = std::real(zr);

        for (i = 1; i < (nn + 1); i++) {
            s1 = y[i - 1];
            cy[i - 1] = s1;
            as = std::abs(s1);
            acs = -std::real(zr) + std::log(as);
            nz += 1;
            y[i - 1] = 0.;
            if (acs < -elim) {
                continue;
            }
            cs = -zr + std::log(s1);
            cs = (std::exp(std::real(cs)) / tol) *
                 std::complex<double>(std::cos(std::imag(cs)), std::sin(std::imag(cs)));
            if (!uchk(cs, *ascle, tol)) {
                y[i - 1] = cs;
                nz -= 1;
                ic = i;
            }
        }
        if (n == 1) {
            return nz;
        }
        if (ic <= 1) {
            y[0] = 0.;
            nz = 2;
        }
        if (n == 2) {
            return nz;
        }
        if (nz == 0) {
            return nz;
        }

        fn = fnu + 1.;
        ck = fn * rz;
        s1 = cy[0];
        s2 = cy[1];
        zri = std::imag(zr);
        zd = zr;
        for (i = 3; i < (n + 1); i++) {
            kk = i;
            cs = s2;
            s2 *= ck;
            s2 += s1;
            s1 = cs;
            ck += rz;
            as = std::abs(s2);
            alas = std::log(as);
            acs = alas - xx;
            nz += 1;
            y[i - 1] = 0.;
            if (acs >= -elim) {
                cs = std::log(s2);
                cs -= zd;
                cs = (std::exp(std::real(cs)) / tol) *
                     std::complex<double>(std::cos(std::imag(cs)), std::sin(std::imag(cs)));
                if (!uchk(cs, *ascle, tol)) {
                    y[i - 1] = cs;
                    nz -= 1;
                    if (ic == kk - 1) {
                        nz = kk - 2;
                        for (int i = 0; i < nz; i++) {
                            y[i] = 0.;
                        }
                        return nz;
                    }
                    ic = kk;
                    continue;
                }
            }
            if (alas >= 0.5 * elim) {
                xx -= elim;
                zd = std::complex<double>(xx, zri);
                s1 *= elm;
                s2 *= elm;
            }
        }
        nz = n;
        if (ic == n) {
            nz = n - 1;
        }

        for (int i = 0; i < nz; i++) {
            y[i] = 0.;
        }
        return nz;
    }

    inline void rati(std::complex<double> z, double fnu, int n, std::complex<double> *cy, double tol) {

        //***BEGIN PROLOGUE  ZRATI
        //***REFER TO  ZBESI,ZBESK,ZBESH
        //
        //     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.
        //
        //***ROUTINES CALLED  AZABS,ZDIV
        //***END PROLOGUE  ZRATI

        std::complex<double> cdfnu, pt, p1, p2, rz, t1;
        double ak, amagz, ap1, ap2, arg, az, dfnu, fdnu, flam, fnup, rap1, rho, test, test1;
        int i, id, idnu, inu, itime, k, kk, magz;

        az = std::abs(z);
        inu = (int)fnu;
        idnu = inu + n - 1;
        fdnu = idnu;
        magz = az;
        amagz = magz + 1;
        fnup = std::fmax(amagz, fdnu);
        id = idnu - magz - 1;
        itime = 1;
        k = 1;
        rz = 2.0 / z;
        t1 = fnup * rz;
        p2 = -t1;
        p1 = 1.0;
        t1 += rz;
        if (id > 0) {
            id = 0;
        }
        ap2 = std::abs(p2);
        ap1 = std::abs(p1);

        //
        // THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX
        // GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT P2
        // VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR PREMATURELY.
        //
        arg = (ap2 + ap2) / (ap1 * tol);
        test1 = std::sqrt(arg);
        test = test1;
        rap1 = 1.0 / ap1;
        p1 *= rap1;
        p2 *= rap1;
        ap2 *= rap1;

        while (1) {
            k += 1;
            ap1 = ap2;
            pt = p2;
            p2 = p1 - t1 * p2;
            p1 = pt;
            t1 += rz;
            ap2 = std::abs(p2);
            if (ap1 > test) {
                break;
            }
            if (itime != 2) {
                ak = std::abs(t1) * 0.5;
                flam = ak + std::sqrt(ak * ak - 1.0);
                rho = std::fmin(ap2 / ap1, flam);
                test = test1 * std::sqrt(rho / (rho * rho - 1.0));
                itime = 2;
            }
        }
        kk = k + 1 - id;
        ak = kk;
        dfnu = fnu + n - 1;
        cdfnu = dfnu;
        t1 = ak;
        p1 = 1.0 / ap2;
        p2 = 0.0;

        for (i = 1; i < (kk + 1); i++) {
            pt = p1;
            p1 = rz * (cdfnu + t1) * p1 + p2;
            p2 = pt;
            t1 -= 1.0;
        }

        if (p1 == 0.) {
            p1 = std::complex<double>(tol, tol);
        }

        cy[n - 1] = p2 / p1;
        if (n == 1) {
            return;
        }
        k = n - 1;
        ak = k;
        t1 = ak;
        cdfnu = fnu * rz;

        for (i = 2; i < (n + 1); i++) {
            pt = cdfnu + t1 * rz + cy[k];
            if (pt == 0.0) {
                pt = std::complex<double>(tol, tol);
            }
            cy[k - 1] = 1.0 / pt;
            t1 -= 1.0;
            k -= 1;
        }
        return;
    }

    inline int seri(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, double tol, double elim,
        double alim
    ) {

        //***BEGIN PROLOGUE  ZSERI
        //***REFER TO  ZBESI,ZBESK
        //
        //     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).
        //
        //***ROUTINES CALLED  DGAMLN,D1MACH,ZUCHK,AZABS,ZDIV,AZLOG,ZMLT
        //***END PROLOGUE  ZSERI

        std::complex<double> ak1, ck, coef, crsc, cz, half_z, rz, s1, s2, w[2];
        double aa, acz, ak, arm, ascle, atol, az, dfnu, fnup, rak1, rs, rtr1, s, ss, x;
        int ib, iflag, il, k, m, nn;

        int nz = 0;
        az = std::abs(z);
        if (az == 0.0) {
            y[0] = 0.0;
            if (fnu == 0.) {
                y[0] = 1.0;
            }
            if (n == 1) {
                return nz;
            }
            for (int i = 1; i < n; i++) {
                y[i] = 0.0;
            }
            return nz;
        }
        x = std::real(z);
        arm = 1e3 * d1mach[0];
        rtr1 = std::sqrt(arm);
        crsc = 1.0;
        iflag = 0;
        if (az >= arm) {
            half_z = 0.5 * z;
            cz = 0.0;
            if (az > rtr1) {
                cz = half_z * half_z;
            }
            acz = std::abs(cz);
            nn = n;
            ck = std::log(half_z);
        L10:
            dfnu = fnu + (nn - 1);
            fnup = dfnu + 1.0;
            //
            // UNDERFLOW TEST
            //
            ak1 = ck * dfnu;
            ak = gamln(fnup);
            ak1 -= ak;
            if (kode == 2) {
                ak1 -= x;
            }
            rak1 = std::real(ak1);
            if (rak1 > -elim) {
                goto L30;
            }
        L20:
            nz += 1;
            y[nn - 1] = 0.0;
            //
            // RETURN WITH NZ < 0 IF ABS(Z*Z/4) > FNU+N-NZ-1 COMPLETE
            // THE CALCULATION IN CBINU WITH N=N-ABS(NZ)
            //
            if (acz > dfnu) {
                return -nz;
            }
            nn -= 1;
            if (nn == 0) {
                return nz;
            }
            goto L10;
        L30:
            if (rak1 <= -alim) {
                iflag = 1;
                ss = 1.0 / tol;
                crsc = tol;
                ascle = arm * ss;
            }
            ak = std::imag(ak1);
            aa = std::exp(rak1);
            if (iflag == 1) {
                aa *= ss;
            }
            coef = aa * std::complex<double>(std::cos(ak), std::sin(ak));
            atol = tol * acz / fnup;
            il = (nn > 2 ? 2 : nn);
            for (int i = 1; i < (il + 1); i++) {
                dfnu = fnu + (nn - i);
                fnup = dfnu + 1.0;
                s1 = 1.0;
                if (acz >= tol * fnup) {
                    ak1 = 1.0;
                    ak = fnup + 2.0;
                    s = fnup;
                    aa = 2.0;
                    while (1) {
                        rs = 1.0 / s;
                        ak1 *= cz;
                        ak1 *= rs;
                        s1 += ak1;
                        s += ak;
                        ak += 2.0;
                        aa *= acz;
                        aa *= rs;
                        if (aa <= atol) {
                            break;
                        }
                    }
                }
                s2 = s1 * coef;
                w[i - 1] = s2;
                if (iflag != 0) {
                    if (uchk(s2, ascle, tol)) {
                        goto L20;
                    }
                }
                m = nn - i + 1;
                y[m - 1] = s2 * crsc;
                if (i != il) {
                    coef *= dfnu / half_z;
                }
            }
            if (nn <= 2) {
                return nz;
            }
            k = nn - 2;
            ak = k;
            rz = 2.0 / z;
            if (iflag == 1) {
                goto L80;
            }
            ib = 3;
        L60:
            for (int i = ib; i < (nn + 1); i++) {
                y[k - 1] = (ak + fnu) * rz * y[k] + y[k + 1];
                ak -= 1.0;
                k -= 1;
            }
            return nz;
        L80:
            //
            // RECUR BACKWARD WITH SCALED VALUES
            //
            //
            // EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
            // UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1000
            //
            s1 = w[0];
            s2 = w[1];
            int l;
            for (l = 3; l < (nn + 1); l++) {
                ck = s2;
                s2 = s1 + (ak + fnu) * rz * s2;
                s1 = ck;
                ck = s2 * crsc;
                y[k - 1] = ck;
                ak -= 1.0;
                k -= 1;
                if (std::abs(ck) > ascle) {
                    goto L100;
                }
            }
            return nz;
        L100:
            ib = l + 1;
            if (ib > nn) {
                return nz;
            }
            goto L60;
        }
        nz = n;
        if (fnu == 0.0) {
            nz -= 1;
        }
        y[0] = 0.0;
        if (fnu == 0.) {
            y[0] = 1.0;
        }
        if (n == 1) {
            return nz;
        }
        for (int i = 1; i < n; i++) {
            y[i] = 0.0;
        }
        return nz;
    }

    inline int s1s2(
        std::complex<double> zr, std::complex<double> *s1, std::complex<double> *s2, double ascle, double alim, int *iuf
    ) {

        //***BEGIN PROLOGUE  ZS1S2
        //***REFER TO  ZBESK,ZAIRY
        //
        //     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.
        //
        //***ROUTINES CALLED  AZABS,AZEXP,AZLOG
        //***END PROLOGUE  ZS1S2

        std::complex<double> c1, s1d;
        double aa, aln, as1, as2, xx;
        int nz = 0;
        as1 = std::abs(*s1);
        as2 = std::abs(*s2);
        aa = std::real(*s1);
        aln = std::imag(*s1);

        if ((aa != 0.) || (aln != 0.)) {
            if (as1 != 0.) {
                xx = std::real(zr);
                aln = -xx - xx + std::log(as1);
                s1d = *s1;
                *s1 = 0.;
                as1 = 0.;
                if (aln >= -alim) {
                    c1 = std::log(s1d) - zr - zr;
                    *s1 = std::exp(c1);
                    as1 = std::abs(*s1);
                    *iuf += 1;
                }
            }
        }
        aa = std::fmax(as1, as2);
        if (aa > ascle) {
            return nz;
        }
        *s1 = 0.;
        *s2 = 0.;
        *iuf = 0;
        return 1;
    }

    inline int uchk(std::complex<double> y, double ascle, double tol) {

        //***BEGIN PROLOGUE  ZUCHK
        //***REFER TO ZSERI,ZUOIK,ZUNK1,ZUNK2,ZUNI1,ZUNI2,ZKSCL
        //
        //      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.
        //
        //***ROUTINES CALLED  (NONE)
        //***END PROLOGUE  ZUCHK

        double yr = std::fabs(std::real(y));
        double yi = std::fabs(std::imag(y));
        double ss = std::fmax(yr, yi);
        double st = std::fmin(yr, yi);
        if (st > ascle) {
            return 0;
        } else {
            st /= tol;
            if (ss < st) {
                return 1;
            } else {
                return 0;
            }
        }
    }

    inline void unhj(
        std::complex<double> z, double fnu, int ipmtr, double tol, std::complex<double> *phi, std::complex<double> *arg,
        std::complex<double> *zeta1, std::complex<double> *zeta2, std::complex<double> *asum, std::complex<double> *bsum
    ) {

        //***BEGIN PROLOGUE  ZUNHJ
        //***REFER TO  ZBESI,ZBESK
        //
        //     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.
        //
        //***ROUTINES CALLED  AZABS,ZDIV,AZLOG,AZSQRT,D1MACH
        //***END PROLOGUE  ZUNHJ

        std::complex<double> cfnu, przth, ptfn, rtzta, rzth, suma, sumb;
        std::complex<double> tfn, t2, w, w2, za, zb, zc, zeta, zth;
        double ang, atol, aw2, azth, btol, fn13, fn23, pp, rfn13;
        double rfnu, rfnu2, wi, wr, zci, zcr, zetai, zetar, zthi;
        double zthr, asumr, asumi, bsumr, bsumi, test, ac;
        double ex1 = 1. / 3.;
        double ex2 = 2. / 3.;
        double hpi = 1.57079632679489662;
        double pi = 3.14159265358979324;
        double thpi = 4.71238898038468986;
        int ias, ibs, j, ju, k, kmax, kp1, ks, l, lrp1, l1, l2, m;
        /* array vars */
        std::complex<double> cr[14] = {0.};
        std::complex<double> dr[14] = {0.};
        std::complex<double> up[14] = {0.};
        std::complex<double> p[30] = {0.};
        double ap[30] = {0.};

        rfnu = 1. / fnu;
        //
        // OVERFLOW TEST (Z/FNU TOO SMALL)
        //
        test = d1mach[0] * 1e3;
        ac = fnu * test;
        if ((std::fabs(std::real(z)) <= ac) && (std::fabs(std::imag(z)) <= ac)) {
            *zeta1 = 2.0 * std::fabs(std::log(test)) + fnu;
            *zeta2 = fnu;
            *phi = 1.;
            *arg = 1.;
            return;
        }
        zb = z * rfnu;
        rfnu2 = rfnu * rfnu;
        //
        // COMPUTE IN THE FOURTH QUADRANT
        //
        fn13 = std::pow(fnu, ex1);
        fn23 = fn13 * fn13;
        rfn13 = 1.0 / fn13;
        w2 = 1.0 - zb * zb;
        /* AMOS AZSQRT and C CSQRT differs when imaginary 0.0 swaps sign */
        w2 = 1.0 - zb * zb;
        if (std::imag(w2) == -0.0) {
            w2 = std::real(w2);
        }
        aw2 = std::abs(w2);
        if (aw2 <= 0.25) {
            //
            // POWER SERIES FOR ABS(W2) <= 0.25
            //
            k = 1;
            p[0] = 1.;
            suma = zunhj_gama[0];
            ap[0] = 1.;
            if (aw2 >= tol) {
                for (k = 2; k < 31; k++) {
                    p[k - 1] = p[k - 2] * w2;
                    suma += p[k - 1] * zunhj_gama[k - 1];
                    ap[k - 1] = ap[k - 2] * aw2;
                    if (ap[k - 1] < tol) {
                        break;
                    }
                }
            }
            /* Check for exhausted loop */
            if (k == 31) {
                k = 30;
            }

            kmax = k;
            zeta = w2 * suma;
            *arg = zeta * fn23;
            za = std::sqrt(suma);
            *zeta2 = std::sqrt(w2) * fnu;
            *zeta1 = (*zeta2) * (1. + zeta * za * ex2);
            za = za + za;
            *phi = std::sqrt(za) * rfn13;
            if (ipmtr == 1) {
                return;
            }
            //
            // SUM SERIES FOR ASUM AND BSUM
            //
            sumb = 0.0;
            for (k = 1; k < (kmax + 1); k++) {
                sumb += p[k - 1] * zunhj_beta[k - 1];
            }
            *asum = 0.0;
            *bsum = sumb;
            l1 = 0;
            l2 = 30;
            btol = tol * (std::fabs(std::real(*bsum)) + std::fabs(std::imag(*bsum)));
            atol = tol;
            pp = 1.0;
            ias = 0;
            ibs = 0;
            if (rfnu2 < tol) {
                /* 110 */
                *asum += 1.;
                *bsum *= rfnu * rfn13;
                /* 120 */
                return;
            }
            for (int is = 2; is < 8; is++) {
                atol /= rfnu2;
                pp *= rfnu2;
                if (ias != 1) {
                    suma = 0.0;
                    for (k = 1; k < (kmax + 1); k++) {
                        m = l1 + k;
                        suma += p[k - 1] * zunhj_alfa[m - 1];
                        if (ap[k - 1] < atol) {
                            break;
                        }
                    }
                    *asum += suma * pp;
                    if (pp < tol) {
                        ias = 1;
                    }
                }
                if (ibs != 1) {
                    sumb = 0.0;
                    for (k = 1; k < (kmax + 1); k++) {
                        m = l2 + k;
                        sumb += p[k - 1] * zunhj_beta[m - 1];
                        if (ap[k - 1] < atol) {
                            break;
                        }
                    }
                    *bsum += sumb * pp;
                    if (pp < btol) {
                        ibs = 1;
                    }
                }
                if ((ias == 1) && (ibs == 1)) {
                    break;
                }
                l1 += 30;
                l2 += 30;
            }
            *asum += 1.;
            *bsum *= rfnu * rfn13;
            return;
        } else {
            //
            // ABS(W2) > 0.25
            w = std::sqrt(w2);
            wr = std::real(w);
            wi = std::imag(w);
            if (wr < 0) {
                wr = 0.;
            }
            if (wi < 0) {
                wi = 0.;
            }
            za = (1. + w) / zb;
            zc = std::log(za);
            zcr = std::real(zc);
            zci = std::imag(zc);
            if (zci < 0) {
                zci = 0.;
            }
            if (zci > hpi) {
                zci = hpi;
            }
            if (zcr < 0) {
                zcr = 0.;
            }
            zc = std::complex<double>(zcr, zci);
            zth = (zc - w) * 1.5;
            cfnu = fnu;
            *zeta1 = zc * cfnu;
            *zeta2 = w * cfnu;
            azth = std::abs(zth);
            zthr = std::real(zth);
            zthi = std::imag(zth);
            ang = thpi;
            if ((zthr < 0.) || (zthi >= 0.)) {
                ang = hpi;
                if (zthr != 0.) {
                    ang = std::atan(zthi / zthr);
                    if (zthr < 0.) {
                        ang += pi;
                    }
                }
            }
            pp = std::pow(azth, ex2);
            ang *= ex2;
            zetar = pp * std::cos(ang);
            zetai = pp * std::sin(ang);
            if (zetai < 0.) {
                zetai = 0.;
            }
            zeta = std::complex<double>(zetar, zetai);
            *arg = zeta * fn23;
            rtzta = zth / zeta;
            za = rtzta / w;
            *phi = std::sqrt(za + za) * rfn13;
            if (ipmtr == 1) {
                return;
            }
            tfn = rfnu / w;
            rzth = rfnu / zth;
            zc = rzth * zunhj_ar[1];
            t2 = 1. / w2;
            up[1] = (t2 * zunhj_c[1] + zunhj_c[2]) * tfn;
            *bsum = up[1] + zc;
            *asum = 0.;

            if (rfnu < tol) {
                *asum += 1.;
                *bsum *= -rfn13 / rtzta;
                return;
            }

            przth = rzth;
            ptfn = tfn;
            up[0] = 1.0;
            pp = 1.0;
            bsumr = std::real(*bsum);
            bsumi = std::imag(*bsum);
            btol = tol * (std::fabs(bsumr) + std::fabs(bsumi));
            ks = 0;
            kp1 = 2;
            l = 3;
            ias = 0;
            ibs = 0;

            for (int lr = 2; lr < 13; lr += 2) {
                lrp1 = lr + 1;
                //
                // COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN
                // NEXT SUMA AND SUMB
                //
                for (k = lr; k < (lrp1 + 1); k++) {
                    ks += 1;
                    kp1 += 1;
                    l += 1;
                    za = zunhj_c[l - 1];
                    for (j = 2; j < (kp1 + 1); j++) {
                        l += 1;
                        za = za * t2 + zunhj_c[l - 1];
                    }
                    ptfn *= tfn;
                    up[kp1 - 1] = ptfn * za;
                    cr[ks - 1] = przth * zunhj_br[ks];
                    przth *= rzth;
                    dr[ks - 1] = przth * zunhj_ar[ks + 1];
                }
                pp *= rfnu2;
                if (ias != 1) {
                    suma = up[lr];
                    ju = lrp1;
                    for (int jr = 1; jr < lrp1; jr++) {
                        ju -= 1;
                        suma += cr[jr - 1] * up[ju - 1];
                    }
                    *asum += suma;
                    asumr = std::real(*asum);
                    asumi = std::imag(*asum);
                    test = std::fabs(asumr) + std::fabs(asumi);
                    if ((pp < tol) && (test < tol)) {
                        ias = 1;
                    }
                }
                if (ibs != 1) {
                    sumb = up[lr + 1] + up[lr] * zc;
                    ju = lrp1;
                    for (int jr = 1; jr < lrp1; jr++) {
                        ju -= 1;
                        sumb += dr[jr - 1] * up[ju - 1];
                    }
                    *bsum += sumb;
                    bsumr = std::real(*bsum);
                    bsumi = std::imag(*bsum);
                    test = std::fabs(bsumr) + std::fabs(bsumi);
                    if ((pp < tol) && (test < tol)) {
                        ibs = 1;
                    }
                }
                if ((ias == 1) && (ibs == 1)) {
                    break;
                }
            }
            *asum += 1.;
            *bsum *= -rfn13 / rtzta;
            return;
        }
    }

    inline void uni1(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, int *nz, int *nlast, double fnul,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZUNI1
        //***REFER TO  ZBESI,ZBESK
        //
        //     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.
        //
        //     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
        //
        //***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,AZABS
        //***END PROLOGUE  ZUNI1

        std::complex<double> c2, phi, rz, sum, s1, s2, zeta1 = 0, zeta2 = 0;
        double aphi, ascle, c1r, crsc, cscl, fn, rs1;
        int i, iflag, init, k, m, nd, nn, nuf;
        std::complex<double> cwrk[16] = {0.};
        std::complex<double> cy[2] = {0.};
        *nz = 0;
        nd = n;
        *nlast = 0;
        //
        // COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN
        // MAGNITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
        // EXP(ALIM)=EXP(ELIM)*TOL
        //
        cscl = 1.0 / tol;
        crsc = tol;
        double css[3] = {cscl, 1., crsc};
        double csr[3] = {crsc, 1., cscl};
        double bry[3] = {1e3 * d1mach[0] / tol, 0., 0.};
        bry[1] = 1.0 / bry[0];
        bry[2] = d1mach[1];
        //
        // CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
        //
        fn = std::fmax(fnu, 1.0);
        init = 0;
        unik(z, fn, 1, 1, tol, &init, &phi, &zeta1, &zeta2, &sum, &cwrk[0]);
        if (kode != 1) {
            s1 = -zeta1 + fn * (fn / (z + zeta2));
        } else {
            s1 = -zeta1 + zeta2;
        }

        rs1 = std::real(s1);
        if (std::fabs(rs1) > elim) {
            if (rs1 > 0) {
                *nz = -1;
                return;
            }
            *nz = n;
            for (i = 0; i < n; i++) {
                y[i] = 0.0;
            }
        }
    L30:
        nn = (nd > 2 ? 2 : nd);
        for (i = 1; i < (nn + 1); i++) {
            fn = fnu + nd - i;
            init = 0;
            unik(z, fn, 1, 0, tol, &init, &phi, &zeta1, &zeta2, &sum, &cwrk[0]);
            if (kode != 1) {
                s1 = -zeta1 + fn * (fn / (z + zeta2)) + std::complex<double>(0.0, std::imag(z));
            } else {
                s1 = -zeta1 + zeta2;
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) > elim) {
                goto L110;
            }
            if (i == 1) {
                iflag = 2;
            }
            if (std::fabs(rs1) >= alim) {
                //
                // REFINE TEST AND SCALE
                //
                aphi = std::abs(phi);
                rs1 += std::log(aphi);
                if (std::fabs(rs1) > elim) {
                    goto L110;
                }
                if (i == 1) {
                    iflag = 1;
                }
                if (rs1 >= 0.0) {
                    if (i == 1) {
                        iflag = 3;
                    }
                }
            }
            /* 60 */
            //
            // SCALE S1 IF CABS(S1) < ASCLE
            //
            s2 = phi * sum;
            s1 = std::exp(std::real(s1)) * css[iflag - 1] *
                 std::complex<double>(std::cos(std::imag(s1)), std::sin(std::imag(s1)));
            s2 *= s1;
            if (iflag == 1) {
                if (uchk(s2, bry[0], tol)) {
                    goto L110;
                }
            }
            /* 70 */
            cy[i - 1] = s2;
            m = nd - i + 1;
            y[m - 1] = s2 * csr[iflag - 1];
        }
        /* 80 */
        if (nd > 2) {
            rz = 2.0 / z;
            s1 = cy[0];
            s2 = cy[1];
            c1r = csr[iflag - 1];
            ascle = bry[iflag - 1];
            k = nd - 2;
            fn = k;
            for (i = 3; i < (nd + 1); i++) {
                c2 = s2;
                s2 = s1 + (fnu + fn) * rz * c2;
                s1 = c2;
                c2 = s2 * c1r;
                y[k - 1] = c2;
                k -= 1;
                fn -= 1.0;
                if (iflag >= 3) {
                    continue;
                }
                if (std::fmax(std::fabs(std::real(c2)), std::fabs(std::imag(c2))) <= ascle) {
                    continue;
                }
                iflag += 1;
                ascle = bry[iflag - 1];
                s1 *= c1r;
                s2 = c2;
                s1 *= css[iflag - 1];
                s2 *= css[iflag - 1];
                c1r = csr[iflag - 1];
            }
            /* 90 */
        }
        /* 100 */
        return;
    L110:
        if (rs1 > 0.0) {
            *nz = -1;
            return;
        }
        y[nd - 1] = 0.0;
        *nz += 1;
        nd -= 1;
        if (nd == 0) {
            return;
        }
        nuf = uoik(z, fnu, kode, 1, nd, y, tol, elim, alim);
        if (nuf < 0) {
            *nz = -1;
            return;
        }
        nd -= nuf;
        *nz += nuf;
        if (nd == 0) {
            return;
        }
        fn = fnu + nd - 1;
        if (fn >= fnul) {
            goto L30;
        }
        *nlast = nd;
        return;
    }

    inline void uni2(
        std::complex<double> z, double fnu, int kode, int n, std::complex<double> *y, int *nz, int *nlast, double fnul,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZUNI2
        //***REFER TO  ZBESI,ZBESK
        //
        //     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
        //
        //***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,AZABS
        //***END PROLOGUE  ZUNI2

        std::complex<double> ai, arg, asum, bsum, cfn, cid, crsc, cscl, c1, c2, dai, phi, rz, s1, s2, zb, zeta1, zeta2,
            zn, zar;
        double aarg, ang, aphi, ascle, ay, c2i, c2m, c2r, fn, rs1, yy;
        int i, iflag, in, inu, j, k, nai, nd, ndai, nn, nuf, idum;
        double hpi = 1.57079632679489662;  /* 0.5 pi */
        double aic = 1.265512123484645396; /* log(2 sqrt(pi)) */
        std::complex<double> cip[4] = {1.0, std::complex<double>(0, 1), -1.0, -std::complex<double>(0, 1)};
        std::complex<double> ci = std::complex<double>(0, 1);
        //
        // COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAGNITUDE
        // ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
        // EXP(ALIM) = EXP(ELIM)*TOL
        //
        cscl = 1.0 / tol;
        crsc = tol;
        std::complex<double> csr[3] = {crsc, 1.0, cscl};
        std::complex<double> css[3] = {cscl, 1.0, crsc};
        double bry[3] = {1e3 * d1mach[0] / tol, 0.0, 0.0};
        std::complex<double> cy[2] = {0.0};
        yy = std::imag(z);
        *nz = 0;
        nd = n;
        *nlast = 0;
        //
        // ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
        //
        zn = -z * ci;
        zb = z;
        cid = -ci;
        inu = (int)fnu;
        ang = hpi * (fnu - inu);
        c2 = std::complex<double>(std::cos(ang), std::sin(ang));
        zar = c2;
        in = inu + n - 1;
        in = in % 4;
        c2 *= cip[in];
        if (yy <= 0.0) {
            zn = std::conj(-zn);
            zb = std::conj(zb);
            cid = -cid;
            c2 = std::conj(c2);
        }
        //
        // CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
        //
        fn = std::fmax(fnu, 1.0);
        unhj(zn, fn, 0, tol, &phi, &arg, &zeta1, &zeta2, &asum, &bsum);
        if (kode != 1) {
            cfn = fnu;
            s1 = -zeta1 + cfn * (cfn / (zb + zeta2));
        } else {
            s1 = -zeta1 + zeta2;
        }
        rs1 = std::real(s1);
        if (std::fabs(rs1) > elim) {
            if (rs1 > 0.) {
                *nz = -1;
                return;
            }
            for (i = 0; i < n; i++) {
                y[i] = 0.0;
            }
            return;
        }
    L10:
        nn = (nd > 2 ? 2 : nd);
        i = 1;
        for (i = 1; i < (nn + 1); i++) {
            fn = fnu + (nd - i);
            unhj(zn, fn, 0, tol, &phi, &arg, &zeta1, &zeta2, &asum, &bsum);
            if (kode != 1) {
                cfn = fn;
                ay = std::fabs(yy);
                s1 = -zeta1 + cfn * (cfn / (zb + zeta2)) + ay * std::complex<double>(0, 1);
            } else {
                s1 = -zeta1 + zeta2;
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) > elim) {
                goto L50;
            }
            if (i == 1) {
                iflag = 2;
            }
            if (std::fabs(rs1) >= alim) {
                //
                // REFINE TEST AND SCALE
                //
                aphi = std::abs(phi);
                aarg = std::abs(arg);
                rs1 += std::log(aphi) - 0.25 * std::log(aarg) - aic;
                if (std::fabs(rs1) > elim) {
                    goto L50;
                }
                if (i == 1) {
                    iflag = 1;
                }
                if (rs1 >= 0.0) {
                    if (i == 1) {
                        iflag = 3;
                    }
                }
            }
            //
            // SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
            // EXPONENT EXTREMES
            //
            ai = airy(arg, 0, 2, &nai, &idum);
            dai = airy(arg, 1, 2, &ndai, &idum);
            s2 = phi * (ai * asum + dai * bsum);
            c2r = std::exp(std::real(s1)) * std::real(css[iflag - 1]);
            c2i = std::imag(s1);
            s1 = c2r * std::complex<double>(std::cos(c2i), std::sin(c2i));
            s2 *= s1;
            if (iflag == 1) {
                if (uchk(s2, bry[0], tol)) {
                    goto L50;
                }
            }
            if (yy <= 0.0) {
                s2 = std::conj(s2);
            }
            j = nd - i + 1;
            s2 *= c2;
            cy[i - 1] = s2;
            y[j - 1] = s2 * csr[iflag - 1];
            c2 *= cid;
        }
        if (nd > 2) {
            rz = 2.0 / z;
            bry[1] = 1.0 / bry[0];
            bry[2] = d1mach[1];
            s1 = cy[0];
            s2 = cy[1];
            c1 = csr[iflag - 1];
            ascle = bry[iflag - 1];
            k = nd - 2;
            fn = k;
            for (i = 3; i < (nd + 1); i++) {
                c2 = s2;
                s2 = s1 + (fnu + fn) * rz * s2;
                s1 = c2;
                c2 = s2 * c1;
                y[k - 1] = c2;
                k -= 1;
                fn -= 1.0;
                if (iflag < 3) {
                    c2r = std::fabs(std::real(c2));
                    c2i = std::fabs(std::imag(c2));
                    c2m = std::fmax(c2r, c2i);
                    if (c2m > ascle) {
                        iflag += 1;
                        ascle = bry[iflag - 1];
                        s1 *= c1;
                        s2 = c2;
                        s1 *= css[iflag - 1];
                        s2 *= css[iflag - 1];
                        c1 = csr[iflag - 1];
                    }
                }
            }
        }
        return;
    L50:
        if (rs1 <= 0.0) {
            //
            // SET UNDERFLOW AND UPDATE PARAMETERS
            //
            y[nd - 1] = 0.0;
            *nz += 1;
            nd -= 1;
            if (nd == 0) {
                return;
            }
            nuf = uoik(z, fnu, kode, 1, nd, y, tol, elim, alim);
            if (nuf >= 0) {
                nd -= nuf;
                *nz += nuf;
                if (nd == 0) {
                    return;
                }
                fn = fnu + nd - 1;
                if (fn >= fnul) {
                    // The following was commented out in the original F77 code
                    // C      FN = CIDI
                    // C      J = NUF + 1
                    // C      K = MOD(J,4) + 1
                    // C      S1R = CIPR(K)
                    // C      S1I = CIPI(K)
                    // C      IF (FN.LT.0.0D0) S1I = -S1I
                    // C      STR = C2R*S1R - C2I*S1I
                    // C      C2I = C2R*S1I + C2I*S1R
                    // C      C2R = STR
                    in = (inu + nd - 1) % 4;
                    c2 = zar * cip[in];
                    if (yy <= 0.0) {
                        c2 = std::conj(c2);
                    }
                    goto L10;
                }
                *nlast = nd;
                return;
            }
        }
        *nz = -1;
        return;
    }

    inline void unik(
        std::complex<double> zr, double fnu, int ikflg, int ipmtr, double tol, int *init, std::complex<double> *phi,
        std::complex<double> *zeta1, std::complex<double> *zeta2, std::complex<double> *total,
        std::complex<double> *cwrk
    ) {

        //***BEGIN PROLOGUE  ZUNIK
        //***REFER TO  ZBESI,ZBESK
        //
        //        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.
        //
        //***ROUTINES CALLED  ZDIV,AZLOG,AZSQRT,D1MACH
        //***END PROLOGUE  ZUNIK

        std::complex<double> cfn, crfn, s, sr, t, t2, zn;
        double ac, rfn, test, tstr, tsti;
        int i, j, k, l;
        /* ( 1/sqrt(2 PI), sqrt(PI/2) ) */
        double con[2] = {3.98942280401432678e-01, 1.25331413731550025};

        if (*init == 0) {
            rfn = 1. / fnu;
            crfn = rfn;

            tstr = std::real(zr);
            tsti = std::imag(zr);
            test = d1mach[0] * 1e3;
            ac = fnu * test;
            if ((std::fabs(tstr) <= ac) && (std::fabs(tsti) <= ac)) {
                ac = 2.0 * std::fabs(std::log(test)) + fnu;
                *zeta1 = ac;
                *zeta2 = fnu;
                *phi = 1.0;
            }
            t = zr * crfn;
            s = 1.0 + t * t;
            sr = std::sqrt(s);
            cfn = fnu;
            zn = (1. + sr) / t;
            *zeta1 = cfn * std::log(zn);
            *zeta2 = cfn * sr;
            t = 1.0 / sr;
            sr = t * crfn;
            cwrk[15] = std::sqrt(sr);
            *phi = cwrk[15] * con[ikflg - 1];
            if (ipmtr != 0) {
                return;
            }
            t2 = 1. / s;
            cwrk[0] = 1.;
            crfn = 1.;
            ac = 1.;
            l = 1;
            k = 2;
            for (k = 2; k < 16; k++) {
                s = 0.0;
                for (j = 1; j < (k + 1); j++) {
                    l += 1;
                    s = s * t2 + zunik_c[l - 1];
                }
                crfn *= sr;
                cwrk[k - 1] = crfn * s;
                ac *= rfn;
                tstr = std::real(cwrk[k - 1]);
                tsti = std::imag(cwrk[k - 1]);
                test = std::fabs(tstr) + std::fabs(tsti);
                if ((ac < tol) && (test < tol)) {
                    break;
                }
            }
            /* Guard against exhausted loop */
            if (k == 16) {
                k -= 1;
            }
            *init = k;
        }

        *total = 0.0;
        t = 1.0;
        if (ikflg != 2) {

            for (i = 0; i < (*init); i++) {
                *total += cwrk[i];
            }
            *phi = cwrk[15] * con[0];

        } else {

            for (i = 1; i < (*init + 1); i++) {
                *total += t * cwrk[i - 1];
                t = -t;
            }
            *phi = cwrk[15] * con[1];
        }
        return;
    }

    inline int unk1(
        std::complex<double> z, double fnu, int kode, int mr, int n, std::complex<double> *y, double tol, double elim,
        double alim
    ) {

        //***BEGIN PROLOGUE  ZUNK1
        //***REFER TO  ZBESK
        //
        //     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
        //
        //***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,AZABS
        //***END PROLOGUE  ZUNK1

        std::complex<double> cfn, ck, crsc, cs, cscl, csgn, cspn, c1, c2, rz, s1, s2, zr, phid, zeta1d = 0.0,
                                                                                                zeta2d = 0.0, sumd;
        double ang, aphi, asc, ascle, c2i, c2m, c2r, fmr, fn, fnf, rs1, sgn, x;
        int i, ib, iflag = 0, ifn, il, inu, iuf, k, kdflg, kflag, kk, m, nw, nz, j, jc, ipard, initd = 0, ic;

        cscl = 1.0 / tol;
        crsc = tol;
        std::complex<double> css[3] = {cscl, 1.0, crsc};
        std::complex<double> csr[3] = {crsc, 1.0, cscl};
        std::complex<double> cwrk[3][16] = {{0.0}};
        std::complex<double> phi[2] = {0.0};
        std::complex<double> sum[2] = {0.0};
        std::complex<double> zeta1[2] = {0.0};
        std::complex<double> zeta2[2] = {0.0};
        std::complex<double> cy[2] = {0.0};
        double bry[3] = {1e3 * d1mach[0] / tol, tol / 1e3 * d1mach[0], d1mach[1]};
        int init[2] = {0};
        double pi = 3.14159265358979324;

        kdflg = 1;
        kflag = 1;
        fn = fnu;
        nz = 0;
        x = std::real(z);
        zr = z;
        if (x < 0.0) {
            zr = -z;
        }
        j = 2;
        for (i = 1; i < (n + 1); i++) {
            j = 3 - j;  /* j flip flops between 1, 2 */
            jc = j - 1; /* dummy index for 0-indexing */
            fn = fnu + (i - 1);
            init[jc] = 0;
            unik(zr, fn, 2, 0, tol, &init[jc], &phi[jc], &zeta1[jc], &zeta2[jc], &sum[jc], &cwrk[jc][0]);
            if (kode != 1) {
                cfn = fn;
                s1 = zeta1[jc] - cfn * (cfn / (zr + zeta2[jc]));
            } else {
                s1 = zeta1[jc] - zeta2[jc];
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) <= elim) {
                if (kdflg == 1) {
                    kflag = 2;
                }
                if (std::fabs(rs1) >= alim) {
                    //
                    // REFINE TEST AND SCALE
                    //
                    aphi = std::abs(phi[jc]);
                    rs1 += std::log(aphi);
                    if (std::fabs(rs1) > elim) {
                        goto L10;
                    }
                    if (kdflg == 1) {
                        kflag = 1;
                    }
                    if (rs1 >= 0.0) {
                        if (kdflg == 1) {
                            kflag = 3;
                        }
                    }
                }

                //
                // SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
                // EXPONENT EXTREMES
                //
                s2 = phi[jc] * sum[jc];
                c2r = std::real(s1);
                c2i = std::imag(s1);
                c2m = std::exp(c2r) * std::real(css[kflag - 1]);
                s1 = c2m * std::complex<double>(std::cos(c2i), std::sin(c2i));
                s2 *= s1;
                if (!((kflag == 1) && (uchk(s2, bry[0], tol)))) {
                    cy[kdflg - 1] = s2;
                    y[i - 1] = s2 * csr[kflag - 1];
                    if (kdflg == 2) {
                        break;
                    }
                    kdflg = 2;
                    continue;
                }
            }
        L10:
            if (rs1 > 0.0) {
                return -1;
            }
            if (x < 0.0) {
                return -1;
            }
            kdflg = 1;
            y[i - 1] = 0.0;
            nz += 1;
            if (i > 1) {
                if (y[i - 2] != 0.0) {
                    y[i - 2] = 0.0;
                    nz += 1;
                }
            }
        }
        /* Check for exhausted loop */
        if (i == (n + 1)) {
            i = n;
        }

        rz = 2.0 / zr;
        ck = fn * rz;
        ib = i + 1;
        if (n >= ib) {
            //
            // TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
            // ON UNDERFLOW
            //
            fn = fnu + (n - 1);
            ipard = 1;
            if (mr != 0) {
                ipard = 0;
            }
            initd = 0;
            unik(zr, fn, 2, ipard, tol, &initd, &phid, &zeta1d, &zeta2d, &sumd, &cwrk[2][0]);
            if (kode != 1) {
                cfn = fn;
                s1 = zeta1d - cfn * (cfn / (zr + zeta2d));
            } else {
                s1 = zeta1d - zeta2d;
            }
            rs1 = std::real(s1);
            if (std::fabs(rs1) <= elim) {
                if (std::fabs(rs1) < alim) {
                    goto L50;
                }
                //
                // REFINE ESTIMATE AND TEST
                //
                aphi = std::abs(phid);
                rs1 += std::log(aphi);
                if (std::fabs(rs1) < elim) {
                    goto L50;
                }
            }
            if (rs1 > 0.0) {
                return -1;
            }
            //
            // FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
            //
            if (x < 0.0) {
                return -1;
            }
            nz = n;
            for (i = 0; i < n; i++) {
                y[i] = 0.0;
            }
            return nz;
        L50:
            //
            // RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
            //
            s1 = cy[0];
            s2 = cy[1];
            c1 = csr[kflag - 1];
            ascle = bry[kflag - 1];
            for (i = ib; i < (n + 1); i++) {
                c2 = s2;
                s2 = ck * s2 + s1;
                s1 = c2;
                ck += rz;
                c2 = s2 * c1;
                y[i - 1] = c2;
                if (kflag < 3) {
                    c2m = std::fmax(std::fabs(std::real(c2)), std::fabs(std::imag(c2)));
                    if (c2m > ascle) {
                        kflag += 1;
                        ascle = bry[kflag - 1];
                        s1 *= c1;
                        s2 = c2;
                        s1 *= css[kflag - 1];
                        s2 *= css[kflag - 1];
                        c1 = csr[kflag - 1];
                    }
                }
            }
        }
        if (mr == 0) {
            return nz;
        }
        //
        // ANALYTIC CONTINUATION FOR RE(Z) < 0.0
        //
        nz = 0;
        fmr = mr;
        sgn = (fmr < 0.0 ? pi : -pi);
        //
        // CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
        //
        csgn = std::complex<double>(0.0, sgn);
        inu = (int)fnu;
        fnf = fnu - inu;
        ifn = inu + n - 1;
        ang = fnf * sgn;
        cspn = std::complex<double>(std::cos(ang), std::sin(ang));
        if (ifn % 2 == 1) {
            cspn = -cspn;
        }
        asc = bry[0];
        kk = n;
        iuf = 0;
        kdflg = 1;
        ib -= 1;
        ic = ib - 1;
        k = 1;
        for (k = 1; k < (n + 1); k++) {
            fn = fnu + (kk - 1);
            //
            // LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
            // FUNCTION ABOVE
            //
            m = 3;
            if (n > 2) {
                goto L80;
            }
        L70:
            initd = init[j - 1];
            phid = phi[j - 1];
            zeta1d = zeta1[j - 1];
            zeta2d = zeta2[j - 1];
            sumd = sum[j - 1];
            m = j;
            j = 3 - j;
            goto L90;
        L80:
            if (!((kk == n) && (ib < n))) {
                if ((kk == ib) || (kk == ic)) {
                    goto L70;
                }
                initd = 0;
            }
        L90:
            unik(zr, fn, 1, 0, tol, &initd, &phid, &zeta1d, &zeta2d, &sumd, &cwrk[m - 1][0]);
            if (kode != 1) {
                cfn = fn;
                s1 = -zeta1d + cfn * (cfn / (zr + zeta2d));
            } else {
                s1 = -zeta1d + zeta2d;
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) > elim) {
                goto L110;
            }
            if (kdflg == 1) {
                iflag = 2;
            }
            if (std::fabs(rs1) >= alim) {
                //
                // REFINE TEST AND SCALE
                //
                aphi = std::abs(phid);
                rs1 += std::log(aphi);
                if (std::fabs(rs1) > elim) {
                    goto L110;
                }
                if (kdflg == 1) {
                    iflag = 1;
                }
                if (rs1 >= 0.0) {
                    if (kdflg == 1) {
                        iflag = 3;
                    }
                }
            }

            s2 = csgn * phid * sumd;
            c2r = std::real(s1);
            c2i = std::imag(s1);
            c2m = std::exp(c2r) * std::real(css[iflag - 1]);
            s1 = c2m * std::complex<double>(std::cos(c2i), std::sin(c2i));
            s2 *= s1;
            if (iflag == 1) {
                if (uchk(s2, bry[0], tol)) {
                    s2 = 0.0;
                }
            }
        L100:
            cy[kdflg - 1] = s2;
            c2 = s2;
            s2 *= csr[iflag - 1];
            //
            // ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
            //
            s1 = y[kk - 1];
            if (kode != 1) {
                nw = s1s2(zr, &s1, &s2, asc, alim, &iuf);
                nz += nw;
            }
            y[kk - 1] = s1 * cspn + s2;
            kk -= 1;
            cspn = -cspn;
            if (c2 == 0.0) {
                kdflg = 1;
                continue;
            }
            if (kdflg == 2) {
                goto L130;
            }
            kdflg = 2;
            continue;
        L110:
            if (rs1 > 0.0) {
                return -1;
            }
            s2 = 0.0;
            goto L100;
        }
        /* If loop is exhausted */
        if (k == n + 1) {
            k -= 1;
        }
    L130:
        il = n - k;
        if (il == 0) {
            return nz;
        };
        s1 = cy[0];
        s2 = cy[1];
        cs = csr[iflag - 1];
        ascle = bry[iflag - 1];
        fn = inu + il;
        for (i = 1; i < (il + 1); i++) {
            c2 = s2;
            s2 = s1 + (fn + fnf) * rz * s2;
            s1 = c2;
            fn -= 1.0;
            c2 = s2 * cs;
            ck = c2;
            c1 = y[kk - 1];
            if (kode != 1) {
                nw = s1s2(zr, &c1, &c2, asc, alim, &iuf);
                nz = nz + nw;
            }
            y[kk - 1] = c1 * cspn + c2;
            kk -= 1;
            cspn = -cspn;
            if (iflag < 3) {
                c2m = std::fmax(std::fabs(std::real(c2)), std::fabs(std::imag(c2)));
                if (c2m > ascle) {
                    iflag += 1;
                    ascle = bry[iflag - 1];
                    s1 *= cs;
                    s2 = ck;
                    s1 *= css[iflag - 1];
                    s2 *= css[iflag - 1];
                    cs = csr[iflag - 1];
                }
            }
        }
        return nz;
    }

    inline int unk2(
        std::complex<double> z, double fnu, int kode, int mr, int n, std::complex<double> *y, double tol, double elim,
        double alim
    ) {

        //***BEGIN PROLOGUE  ZUNK2
        //***REFER TO  ZBESK
        //
        //     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
        //
        //***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,AZABS
        //***END PROLOGUE  ZUNK2

        std::complex<double> ai, cfn, ck, cs, csgn, cspn, c1, c2, dai, rz, s1, s2, zb, zn, zr, phid, argd, zeta1d,
            zeta2d, asumd, bsumd;
        double aarg, ang, aphi, asc, ascle, car, cpn, c2i, c2m, c2r, crsc, cscl, fmr, fn, fnf, rs1, sar, sgn, spn, x,
            yy;
        int i, ib, iflag = 0, ifn, il, in, inu, iuf, k, kdflg, kflag, kk, nai, ndai, nw, nz, idum, j, ipard, ic;

        std::complex<double> cr1 = std::complex<double>(1.0, 1.73205080756887729);      /*   1 + sqrt(3)i  */
        std::complex<double> cr2 = std::complex<double>(-0.5, -8.66025403784438647e-1); /*      0.5 cr1    */
        double hpi = 1.57079632679489662;                                               /*      0.5 pi     */
        double pi = 3.14159265358979324;
        double aic = 1.26551212348464539; /* log(2 sqrt(pi)) */
        std::complex<double> cip[4] = {1.0, -std::complex<double>(0, 1), -1.0, std::complex<double>(0, 1)};
        cscl = 1.0 / tol;
        crsc = tol;
        std::complex<double> css[3] = {cscl, 1.0, crsc};
        std::complex<double> csr[3] = {crsc, 1.0, cscl};
        std::complex<double> phi[2] = {0.0};
        std::complex<double> arg[2] = {0.0};
        std::complex<double> zeta1[2] = {0.0};
        std::complex<double> zeta2[2] = {0.0};
        std::complex<double> asum[2] = {0.0};
        std::complex<double> bsum[2] = {0.0};
        std::complex<double> cy[2] = {0.0};
        double bry[3] = {(1.0 + 1e3 * d1mach[0] / tol), 1.0 / (1.0 + 1e3 * d1mach[0] / tol), d1mach[1]};

        kdflg = 1;
        kflag = 1;
        fn = fnu;
        nz = 0;
        //
        // EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
        // THE UNDERFLOW LIMIT
        //
        x = std::real(z);
        zr = z;
        if (x < 0.0) {
            zr = -z;
        }
        yy = std::imag(zr);
        zn = -zr * std::complex<double>(0, 1);
        zb = zr;
        inu = (int)fnu;
        fnf = fnu - inu;
        ang = -hpi * fnf;
        car = std::cos(ang);
        sar = std::sin(ang);
        cpn = hpi * car;
        spn = hpi * sar;
        c2 = std::complex<double>(spn, -cpn);
        kk = (inu % 4) + 1;
        cs = cr1 * c2 * cip[kk - 1];
        if (yy <= 0.0) {
            zn = std::conj(-zn);
            zb = std::conj(zb);
        }
        //
        // K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
        // QUADRANT.  FOURTH QUADRANT VALUES (YY <= 0.0_dp) ARE COMPUTED BY
        // CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
        //
        j = 2;
        for (i = 1; i < (n + 1); i++) {
            j = 3 - j;
            fn = fnu + (i - 1);
            unhj(zn, fn, 0, tol, &phi[j - 1], &arg[j - 1], &zeta1[j - 1], &zeta2[j - 1], &asum[j - 1], &bsum[j - 1]);
            if (kode != 1) {
                cfn = fn;
                s1 = zeta1[j - 1] - cfn * (cfn / (zb + zeta2[j - 1]));
            } else {
                s1 = zeta1[j - 1] - zeta2[j - 1];
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) <= elim) {
                if (kdflg == 1) {
                    kflag = 2;
                }
                if (std::fabs(rs1) >= alim) {
                    //
                    // REFINE TEST AND SCALE
                    //
                    aphi = std::abs(phi[j - 1]);
                    aarg = std::abs(arg[j - 1]);
                    rs1 += std::log(aphi) - 0.25 * std::log(aarg) - aic;
                    if (std::fabs(rs1) > elim) {
                        /* GO TO 70 */
                        if (rs1 > 0.0) {
                            return -1;
                        }
                        /* FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW */
                        if (x < 0.0) {
                            return -1;
                        }
                        kdflg = 1;
                        y[i - 1] = 0.0;
                        cs *= -std::complex<double>(0, 1);
                        nz += 1;
                        if (i != 1) {
                            if (y[i - 2] != 0.0) {
                                y[i - 2] = 0.0;
                                nz += 1;
                            }
                        }
                        continue;
                    }
                    if (kdflg == 1) {
                        kflag = 1;
                    }
                    if (rs1 >= 0.0) {
                        if (kdflg == 1) {
                            kflag = 3;
                        }
                    }
                }
                //
                // SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
                // EXPONENT EXTREMES
                //
                c2 = arg[j - 1] * cr2;
                ai = airy(c2, 0, 2, &nai, &idum);
                dai = airy(c2, 1, 2, &ndai, &idum);
                s2 = cs * phi[j - 1] * (ai * asum[j - 1] + cr2 * dai * bsum[j - 1]);
                c2r = std::real(s1);
                c2i = std::imag(s1);
                c2m = std::exp(c2r) * std::real(css[kflag - 1]);
                s1 = c2m * std::complex<double>(std::cos(c2i), std::sin(c2i));
                s2 *= s1;
                if (kflag == 1) {
                    if (uchk(s2, bry[0], tol)) {
                        /* GO TO 70 */
                        if (rs1 > 0.0) {
                            return -1;
                        }
                        /* FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW */
                        if (x < 0.0) {
                            return -1;
                        }
                        kdflg = 1;
                        y[i - 1] = 0.0;
                        cs *= -std::complex<double>(0, 1);
                        nz += 1;
                        if (i != 1) {
                            if (y[i - 2] != 0.0) {
                                y[i - 2] = 0.0;
                                nz += 1;
                            }
                        }
                        continue;
                    }
                }
                if (yy <= 0.0) {
                    s2 = std::conj(s2);
                }
                cy[kdflg - 1] = s2;
                y[i - 1] = s2 * csr[kflag - 1];
                cs *= -std::complex<double>(0, 1);
                if (kdflg == 2) {
                    break;
                }
                kdflg = 2;
                continue;
            }
            /* GO TO 70 */
            if (rs1 > 0.0) {
                return -1;
            }
            /* FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW */
            if (x < 0.0) {
                return -1;
            }
            kdflg = 1;
            y[i - 1] = 0.0;
            cs *= -std::complex<double>(0, 1);
            nz += 1;
            if (i != 1) {
                if (y[i - 2] != 0.0) {
                    y[i - 2] = 0.0;
                    nz += 1;
                }
            }
            continue;
        }
        /* Check for exhausted loop */
        if (i == n + 1) {
            i = n;
        }

        rz = 2.0 / zr;
        ck = fn * rz;
        ib = i + 1;
        if (n >= ib) {
            fn = fnu + (n - 1);
            ipard = 1;
            if (mr != 0) {
                ipard = 0;
            }
            unhj(zn, fn, ipard, tol, &phid, &argd, &zeta1d, &zeta2d, &asumd, &bsumd);
            if (kode != 1) {
                cfn = fn;
                s1 = zeta1d - cfn * (cfn / (zb + zeta2d));
            } else {
                s1 = zeta1d - zeta2d;
            }
            rs1 = std::real(s1);
            if (std::fabs(rs1) <= elim) {
                if (std::fabs(rs1) < alim) {
                    goto L120;
                }
                //
                // REFINE ESTIMATE AND TEST
                //
                aphi = std::abs(phid);
                rs1 += std::log(aphi);
                if (std::fabs(rs1) < elim) {
                    goto L120;
                }
            }
            if (rs1 > 0.0) {
                return -1;
            }
            //
            //  FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
            //
            if (x < 0.0) {
                return -1;
            }
            nz = n;
            for (i = 0; i < n; i++) {
                y[i] = 0.0;
            }
            return nz;
        L120:
            //
            // SCALED FORWARD RECURRENCE FOR REMAINDER OF THE SEQUENCE
            //
            s1 = cy[0];
            s2 = cy[1];
            c1 = csr[kflag - 1];
            ascle = bry[kflag - 1];
            for (i = ib; i < (n + 1); i++) {
                c2 = s2;
                s2 = ck * s2 + s1;
                s1 = c2;
                ck += rz;
                c2 = s2 * c1;
                y[i - 1] = c2;
                if (kflag < 3) {
                    c2m = std::fmax(std::fabs(std::real(c2)), std::fabs(std::imag(c2)));
                    if (c2m > ascle) {
                        kflag += 1;
                        ascle = bry[kflag - 1];
                        s1 *= c1;
                        s2 = c2;
                        s1 *= css[kflag - 1];
                        s2 *= css[kflag - 1];
                        c1 = csr[kflag - 1];
                    }
                }
            }
        }
        if (mr == 0) {
            return nz;
        }
        //
        // ANALYTIC CONTINUATION FOR RE(Z) < 0.0
        //
        nz = 0;
        fmr = mr;
        sgn = (fmr < 0.0 ? pi : -pi);
        //
        // CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
        //
        csgn = std::complex<double>(0.0, sgn);
        if (yy <= 0.0) {
            csgn = -csgn;
        }
        ifn = inu + n - 1;
        ang = fnf * sgn;
        cspn = std::complex<double>(std::cos(ang), std::sin(ang));
        if (ifn % 2 == 1) {
            cspn = -cspn;
        }
        //
        // CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION.  I(FNU,Z) IS
        // COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
        // QUADRANT.  FOURTH QUADRANT VALUES (YY <= 0.0_dp) ARE COMPUTED BY
        // CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
        //
        cs = std::complex<double>(car, -sar) * csgn;
        in = (ifn % 4) + 1;
        c2 = cip[in - 1];
        cs *= std::conj(c2);
        asc = bry[0];
        iuf = 0;
        kk = n;
        kdflg = 1;
        ib -= 1;
        ic = ib - 1;
        for (k = 1; k < (n + 1); k++) {
            fn = fnu + (kk - 1);
            if (n > 2) {
                goto L175;
            }
        L172:
            phid = phi[j - 1];
            argd = arg[j - 1];
            zeta1d = zeta1[j - 1];
            zeta2d = zeta2[j - 1];
            asumd = asum[j - 1];
            bsumd = bsum[j - 1];
            j = 3 - j;
            goto L210;
        L175:
            if (!((kk == n) && (ib < n))) {
                if ((kk == ib) || (kk == ic)) {
                    goto L172;
                }
                unhj(zn, fn, 0, tol, &phid, &argd, &zeta1d, &zeta2d, &asumd, &bsumd);
            }
        L210:
            if (kode != 1) {
                cfn = fn;
                s1 = -zeta1d + cfn * (cfn / (zb + zeta2d));
            } else {
                s1 = -zeta1d + zeta2d;
            }
            //
            // TEST FOR UNDERFLOW AND OVERFLOW
            //
            rs1 = std::real(s1);
            if (std::fabs(rs1) > elim) {
                if (rs1 > 0.0) {
                    return -1;
                }
                s2 = 0.0;
                goto L250;
            }
            if (kdflg == 1) {
                iflag = 2;
            }
            if (std::fabs(rs1) >= alim) {
                //
                // REFINE  TEST AND SCALE
                //
                aphi = std::abs(phid);
                aarg = std::abs(argd);
                rs1 += std::log(aphi) - 0.25f * std::log(aarg) - aic;
                if (std::fabs(rs1) > elim) {
                    if (rs1 > 0.0) {
                        return -1;
                    }
                    s2 = 0.0;
                    goto L250;
                }
                if (kdflg == 1) {
                    iflag = 1;
                }
                if (rs1 >= 0.0) {
                    if (kdflg == 1) {
                        iflag = 3;
                    }
                }
            }

            ai = airy(argd, 0, 2, &nai, &idum);
            dai = airy(argd, 1, 2, &ndai, &idum);
            s2 = cs * phid * (ai * asumd + dai * bsumd);
            c2r = std::real(s1);
            c2i = std::imag(s1);
            c2m = std::exp(c2r) * std::real(css[iflag - 1]);
            s1 = c2m * std::complex<double>(std::cos(c2i), std::sin(c2i));
            s2 *= s1;
            if (iflag == 1) {
                if (uchk(s2, bry[0], tol)) {
                    s2 = 0.0;
                }
            }
        L250:
            if (yy <= 0.0) {
                s2 = std::conj(s2);
            }
            cy[kdflg - 1] = s2;
            c2 = s2;
            s2 *= csr[iflag - 1];
            //
            // ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
            //
            s1 = y[kk - 1];
            if (kode != 1) {
                nw = s1s2(zr, &s1, &s2, asc, alim, &iuf);
                nz += nw;
            }
            y[kk - 1] = s1 * cspn + s2;
            kk -= 1;
            cspn = -cspn;
            cs *= -std::complex<double>(0, 1);
            if (c2 == 0.0) {
                kdflg = 1;
                continue;
            }
            if (kdflg == 2) {
                break;
            }
            kdflg = 2;
            continue;
        }
        /* Check for exhausted loop */
        if (k == n + 1) {
            k = n;
        }

        il = n - k;
        if (il == 0) {
            return nz;
        }
        //
        // RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
        // K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
        // INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
        //
        s1 = cy[0];
        s2 = cy[1];
        cs = csr[iflag - 1];
        ascle = bry[iflag - 1];
        fn = inu + il;
        for (i = 1; i < (il + 1); i++) {
            c2 = s2;
            s2 = s1 + (fn + fnf) * rz * c2;
            s1 = c2;
            fn -= 1.0;
            c2 = s2 * cs;
            ck = c2;
            c1 = y[kk - 1];
            if (kode != 1) {
                nw = s1s2(zr, &c1, &c2, asc, alim, &iuf);
                nz = nz + nw;
            }
            y[kk - 1] = c1 * cspn + c2;
            kk -= 1;
            cspn = -cspn;
            if (iflag < 3) {
                c2m = std::fmax(std::fabs(std::real(ck)), std::fabs(std::imag(ck)));
                if (c2m > ascle) {
                    iflag += 1;
                    ascle = bry[iflag - 1];
                    s1 *= cs;
                    s2 = ck;
                    s1 *= css[iflag - 1];
                    s2 *= css[iflag - 1];
                    cs = csr[iflag - 1];
                }
            }
        }
        return nz;
    }

    inline int uoik(
        std::complex<double> z, double fnu, int kode, int ikflg, int n, std::complex<double> *y, double tol,
        double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZUOIK
        //***REFER TO  ZBESI,ZBESK,ZBESH
        //
        //     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
        //
        //***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,AZABS,AZLOG
        //***END PROLOGUE  ZUOIK

        std::complex<double> arg, asum, bsum, cz, phi, sum, zb, zeta1;
        std::complex<double> zeta2, zn, zr;
        double aarg, aphi, ascle, ax, ay, fnn, gnn, gnu, rcz, x, yy;
        int iform, init, nn;
        double aic = 1.265512123484645396;
        std::complex<double> cwrk[16] = {0.};

        int nuf = 0;
        nn = n;
        x = std::real(z);
        zr = z;
        if (x < 0.) {
            zr = -z;
        }
        zb = zr;
        yy = std::imag(zr);
        ax = std::fabs(x) * std::sqrt(3.);
        ay = std::fabs(yy);
        iform = 1;
        if (ay > ax) {
            iform = 2;
        }
        gnu = std::fmax(fnu, 1.);
        if (ikflg != 1) {
            fnn = nn;
            gnn = fnu + fnn - 1;
            gnu = std::fmax(gnn, fnn);
        }

        if (iform != 2) {
            init = 0;
            unik(zr, gnu, ikflg, 1, tol, &init, &phi, &zeta1, &zeta2, &sum, &cwrk[0]);
            cz = -zeta1 + zeta2;
        } else {
            zn = -zr * std::complex<double>(0, 1);
            if (yy <= 0.) {
                zn.real(-zn.real());
            }
            unhj(zn, gnu, 1, tol, &phi, &arg, &zeta1, &zeta2, &asum, &bsum);
            cz = zeta2 - zeta1;
            aarg = std::abs(arg);
        }
        if (kode == 2) {
            cz -= zb;
        }
        if (ikflg == 2) {
            cz = -cz;
        }
        aphi = std::abs(phi);
        rcz = std::real(cz);

        /*  OVERFLOW TEST  */
        if (rcz > elim) {
            return -1;
        }
        if (rcz >= alim) {
            rcz += std::log(aphi);
            if (iform == 2) {
                rcz -= 0.25 * std::log(aarg) + aic;
            }
            if (rcz > elim) {
                return -1;
            }
        } else {
            /*  UNDERFLOW TEST  */
            if (rcz >= -elim) {
                if (rcz > -alim) {
                    /* pass */
                } else {
                    rcz += std::log(aphi);
                    if (iform == 2) {
                        rcz -= 0.25 * std::log(aarg) + aic;
                    }
                    if (rcz > -elim) {
                        /* goto 30 */
                        ascle = 1e3 * d1mach[0] / tol;
                        cz += std::log(phi);
                        if (iform != 1) {
                            cz -= 0.25 * std::log(arg) + aic;
                        }
                        ax = std::exp(rcz) / tol;
                        ay = std::imag(cz);
                        cz = ax * std::exp(ay);
                        if (uchk(cz, ascle, tol)) {
                            for (int i = 0; i < nn; i++) {
                                y[i] = 0.;
                            }
                            return nn;
                        }
                    } else {
                        for (int i = 0; i < nn; i++) {
                            y[i] = 0.;
                        }
                        return nn;
                    }
                }
            } else {
                for (int i = 0; i < nn; i++) {
                    y[i] = 0.;
                }
                return nn;
            }
        }
        if ((ikflg == 2) || (n == 1)) {
            return nuf;
        }
        /* 140 */
        while (1) {
            gnu = fnu + (nn - 1);
            if (iform != 2) {
                init = 0;
                unik(zr, gnu, ikflg, 1, tol, &init, &phi, &zeta1, &zeta2, &sum, &cwrk[0]);
                cz = zeta2 - zeta1;
            } else {
                unhj(zn, gnu, 1, tol, &phi, &arg, &zeta1, &zeta2, &asum, &bsum);
                cz = zeta2 - zeta1;
                aarg = std::abs(phi);
            }
            if (kode == 2) {
                cz -= zb;
            }

            /* 170 */
            aphi = std::abs(phi);
            rcz = std::real(cz);

            if (rcz >= -elim) {
                if (rcz > -alim) {
                    return nuf;
                }
                rcz += std::log(aphi);
                if (iform == 2) {
                    rcz -= 0.25 * std::log(aarg) + aic;
                }
                if (rcz > -elim) {
                    ascle = 1e3 * d1mach[0] / tol;
                    cz += std::log(phi);
                    if (iform != 1) {
                        cz -= 0.25 * std::log(arg) + aic;
                    }
                    ax = std::exp(rcz) / tol;
                    ay = std::imag(cz);
                    cz = ax * (std::cos(ay) + std::sin(ay) * std::complex<double>(0, 1));
                    if (!(uchk(cz, ascle, tol))) {
                        return nuf;
                    }
                }
            }

            y[nn - 1] = 0.;
            nn -= 1;
            nuf += 1;
            if (nn == 0) {
                return nuf;
            }
        }
        return -1;
    }

    inline int wrsk(
        std::complex<double> zr, double fnu, int kode, int n, std::complex<double> *y, std::complex<double> *cw,
        double tol, double elim, double alim
    ) {

        //***BEGIN PROLOGUE  ZWRSK
        //***REFER TO  ZBESI,ZBESK
        //
        //     ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
        //     NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
        //
        //***ROUTINES CALLED  D1MACH,ZBKNU,ZRATI,AZABS
        //***END PROLOGUE  ZWRSK

        std::complex<double> cinu, cscl, ct, c1, c2, rct, st;
        double act, acw, ascle, yy;
        int i, nw, nz;

        //
        // I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
        // Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
        // WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
        //
        nz = 0;
        nw = bknu(zr, fnu, kode, 2, cw, tol, elim, alim);
        if (nw != 0) {
            /* 50 */
            nz = -1;
            if (nw == -2) {
                nz = -2;
            }
            return nz;
        }
        rati(zr, fnu, n, y, tol);
        //
        // RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
        // R(FNU+J-1,Z)=Y(J),  J=1,...,N
        //
        cinu = 1.0;
        if (kode != 1) {
            yy = std::imag(zr);
            cinu = std::complex<double>(std::cos(yy), std::sin(yy));
        }
        //
        // ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH THE
        // UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE SCALED TO
        // PREVENT OVER OR UNDERFLOW.  CUOIK HAS DETERMINED THAT THE RESULT
        // IS ON SCALE.
        //
        acw = std::abs(cw[1]);
        ascle = 1e3 * d1mach[0] / tol;
        cscl = 1.0;

        if (acw <= ascle) {
            cscl = 1.0 / tol;
        } else {
            ascle = 1.0 / ascle;
            if (acw >= ascle) {
                cscl = tol;
            }
        }

        c1 = cw[0] * cscl;
        c2 = cw[1] * cscl;
        st = y[0];
        //
        // CINU=CINU*(CONJG(CT)/ABS(CT))*(1.0_dp/ABS(CT) PREVENTS
        // UNDER- OR OVERFLOW PREMATURELY BY SQUARING ABS(CT)
        //
        ct = zr * (c2 + st * c1);
        act = std::abs(ct);
        rct = 1.0 / act;
        ct = std::conj(ct) * rct;
        cinu *= ct * rct;
        y[0] = cinu * cscl;
        if (n == 1) {
            return nz;
        }
        for (i = 2; i < (n + 1); i++) {
            cinu *= st;
            st = y[i - 1];
            y[i - 1] = cinu * cscl;
        }
        return nz;
    }
} // namespace amos
} // namespace xsf