cbessel.f90

Module complex_bessel
 
!      REMARK ON ALGORITHM 644, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 21, NO. 4, December, 1995, P.  388--393.
 
! Code converted using TO_F90 by Alan Miller
! Date: 2002-02-08  Time: 17:53:05
! Latest revision - 16 April 2002
 
IMPLICIT NONE
INTEGER, PARAMETER, PUBLIC  :: dp = kind(1.0d0)
 
PRIVATE
PUBLIC  :: cbesh, cbesi, cbesj, cbesk, cbesy, cairy, cbiry, gamln
 
 
CONTAINS
 
 
SUBROUTINE cbesh(z, fnu, kode, m, n, cy, nz, ierr)
!***BEGIN PROLOGUE  CBESH
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!   ON KODE=1, CBESH 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.CMPLX(0.0E0,0.0E0) IN THE CUT PLANE -PI < ARG(Z) <= PI.
!   ON KODE=2, CBESH COMPUTES THE SCALED HANKEL FUNCTIONS
 
!   CY(I)=H(M,FNU+J-1,Z)*EXP(-MM*Z*I)       MM=3-2M,      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
!     Z      - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI < ARG(Z) <= PI
!     FNU    - ORDER OF INITIAL H FUNCTION, FNU >= 0.0E0
!     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 OF THE SEQUENCE, N >= 1
 
!   OUTPUT
!     CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
!              VALUES 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 > 0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE TO UNDERFLOW,
!                        CY(J)=CMPLX(0.0,0.0) J=1,...,NZ WHEN Y > 0.0 AND M=1
!                        OR Y < 0.0 AND M=2. FOR THE COMPLEMENTARY 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+N-1 TOO
!                      LARGE OR ABS(Z) TOO SMALL OR BOTH
!              IERR=3, ABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                      BUT LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION
!                      PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!              IERR=4, ABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTATION BECAUSE OF
!                      COMPLETE LOSSES OF SIGNIFICANCE 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) >= 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) > 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)*EXP( PI*FNU*I)
!          H(2,-FNU,Z) = H(2,FNU,Z)*EXP(-PI*FNU*I)
!                    I**2=-1
 
!    CAN BE USED.
 
!    IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELEMENTARY
!    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=R1MACH(4)=UNIT ROUNDOFF.  ALSO 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.8D+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 REPRESENTS THE INCREASE IN ERROR DUE
!    TO ARGUMENT REDUCTION IN THE ELEMENTARY FUNCTIONS.
!    HERE, S=MAX(1, ABS(LOG10(ABS(Z))), ABS(LOG10(FNU))) APPROXIMATELY
!    (I.E. S=MAX(1,ABS(EXPONENT OF ABS(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, ACM TRANS. MATH. SOFTWARE,
!        VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH
!***END PROLOGUE  CBESH
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: zn, zt, csgn
real(dp)     :: 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
INTEGER       :: i, inu, inuh, ir, k, k1, k2, mm, mr, nn, nuf, nw
 
real(dp), PARAMETER  :: hpi = 1.57079632679489662_dp
 
!***FIRST EXECUTABLE STATEMENT  CBESH
nz = 0
xx = real(z, kind=dp)
yy = aimag(z)
ierr = 0
IF (xx == 0.0_dp .AND. yy == 0.0_dp) ierr = 1
IF (fnu < 0.0_dp) ierr = 1
IF (m < 1 .OR. m > 2) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (n < 1) ierr = 1
IF (ierr /= 0) RETURN
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) < 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
!-----------------------------------------------------------------------
tol = max(epsilon(0.0_dp), 1.0d-18)
k1 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1), abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa, 18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa, -41.45_dp)
fnul = 10.0_dp + 6.0_dp * (dig - 3.0_dp)
rl = 1.2_dp * dig + 3.0_dp
fn = fnu + (nn-1)
mm = 3 - m - m
fmm = mm
zn = z * cmplx(0.0_dp, -fmm, kind=dp)
xn = real(zn, kind=dp)
yn = aimag(zn)
az = abs(z)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
IF (az <= aa) THEN
  IF (fn <= aa) THEN
    aa = sqrt(aa)
    IF (az > aa) ierr = 3
    IF (fn > aa) ierr = 3
!-----------------------------------------------------------------------
!     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
!-----------------------------------------------------------------------
    ufl = tiny(0.0_dp) * 1.0d+3
    IF (az >= ufl) THEN
      IF (fnu <= fnul) THEN
        IF (fn > 1.0_dp) THEN
          IF (fn <= 2.0_dp) THEN
            IF (az > tol) GO TO 10
            arg = 0.5_dp * az
            aln = -fn * log(arg)
            IF (aln > elim) GO TO 50
          ELSE
            CALL cuoik(zn, fnu, kode, 2, nn, cy, nuf, tol, elim, alim)
            IF (nuf < 0) GO TO 50
            nz = nz + nuf
            nn = 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 40
          END IF
        END IF
 
        10 IF (.NOT.(xn < 0.0_dp .OR. (xn == 0.0_dp .AND. yn < 0.0_dp  &
               .AND. m == 2))) THEN
!-----------------------------------------------------------------------
!     RIGHT HALF PLANE COMPUTATION, XN >= 0.  .AND.  (XN.NE.0.  .OR.
!     YN >= 0.  .OR.  M=1)
!-----------------------------------------------------------------------
          CALL cbknu(zn, fnu, kode, nn, cy, nz, tol, elim, alim)
          GO TO 20
        END IF
!-----------------------------------------------------------------------
!     LEFT HALF PLANE COMPUTATION
!-----------------------------------------------------------------------
        mr = -mm
        CALL cacon(zn, fnu, kode, mr, nn, cy, nw, rl, fnul, tol, elim, alim)
        IF (nw < 0) GO TO 60
        nz = nw
      ELSE
!-----------------------------------------------------------------------
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU > FNUL
!-----------------------------------------------------------------------
        mr = 0
        IF (.NOT.(xn >= 0.0_dp .AND. (xn /= 0.0_dp .OR. yn >= 0.0_dp  &
               .OR. m /= 2))) THEN
          mr = -mm
          IF (xn == 0.0_dp .AND. yn < 0.0_dp) zn = -zn
        END IF
        CALL cbunk(zn, fnu, kode, mr, nn, cy, nw, tol, elim, alim)
        IF (nw < 0) GO TO 60
        nz = nz + nw
      END IF
!-----------------------------------------------------------------------
!     H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
 
!     ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
!-----------------------------------------------------------------------
      20 sgn = sign(hpi,-fmm)
!-----------------------------------------------------------------------
!     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_dp / sgn
      cpn = rhpi * cos(arg)
      spn = rhpi * sin(arg)
!     ZN = CMPLX(-SPN,CPN)
      csgn = cmplx(-spn, cpn, kind=dp)
!     IF (MOD(INUH,2).EQ.1) ZN = -ZN
      IF (mod(inuh,2) == 1) csgn = -csgn
      zt = cmplx(0.0_dp, -fmm, kind=dp)
      rtol = 1.0_dp / tol
      ascle = ufl * rtol
      DO  i = 1, nn
!       CY(I) = CY(I)*ZN
!       ZN = ZN*ZT
        zn = cy(i)
        aa = real(zn, kind=dp)
        bb = aimag(zn)
        atol = 1.0_dp
        IF (max(abs(aa),abs(bb)) <= ascle) THEN
          zn = zn * rtol
          atol = tol
        END IF
        zn = zn * csgn
        cy(i) = zn * atol
        csgn = csgn * zt
      END DO
      RETURN
 
      40 IF (xn >= 0.0_dp) RETURN
    END IF
 
    50 ierr = 2
    nz = 0
    RETURN
 
    60 IF (nw == -1) GO TO 50
    nz = 0
    ierr = 5
    RETURN
  END IF
END IF
nz = 0
ierr = 4
RETURN
END SUBROUTINE cbesh
 
 
 
SUBROUTINE cbesi(z, fnu, kode, n, cy, nz, ierr)
!***BEGIN PROLOGUE  CBESI
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!   ON KODE=1, CBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
!   BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL (dp), NONNEGATIVE
!   ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE
!   -PI < ARG(Z) <= PI. ON KODE=2, CBESI 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
!     Z      - Z=CMPLX(X,Y),  -PI < ARG(Z) <= PI
!     FNU    - ORDER OF INITIAL I FUNCTION, FNU >= 0.0_dp
!     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 >= 1
 
!   OUTPUT
!     CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
!              VALUES 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 > 0 , LAST NZ COMPONENTS OF CY SET TO ZERO
!                        DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0),
!                        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, ABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                      BUT LOSSES OF SIGNIFICANCE BY ARGUMENT
!                      REDUCTION PRODUCE LESS THAN HALF OF MACHINE
!                      ACCURACY
!              IERR=4, ABS(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 ABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE ABS(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) > 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 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 > ABS(Z).
 
!   IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELEMENTARY
!   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=EPSILON(0.0_dp)=UNIT ROUNDOFF.  ALSO
!   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=HUGE(0). 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.8D+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 REPRESENTS THE INCREASE IN ERROR DUE TO
!   ARGUMENT REDUCTION IN THE ELEMENTARY FUNCTIONS.  HERE, S =
!   MAX(1,ABS(LOG10(ABS(Z))), ABS(LOG10(FNU))) APPROXIMATELY
!   (I.E. S = MAX(1,ABS(EXPONENT OF ABS(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, ACM TRANS. MATH. SOFTWARE,
!         VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CBINU,I1MACH,R1MACH
!***END PROLOGUE  CBESI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: csgn, zn
real(dp)     :: aa, alim, arg, dig, elim, fnul, rl, r1m5, s1, s2,  &
                 tol, xx, yy, az, fn, bb, ascle, rtol, atol
INTEGER       :: i, inu, k, k1, k2, nn
 
real(dp), PARAMETER     :: pi = 3.14159265358979324_dp
COMPLEX (dp), PARAMETER  :: cone = (1.0_dp, 0.0_dp)
 
!***FIRST EXECUTABLE STATEMENT  CBESI
ierr = 0
nz = 0
IF (fnu < 0.0_dp) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (n < 1) ierr = 1
IF (ierr /= 0) RETURN
xx = real(z, kind=dp)
yy = aimag(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) < 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.
!-----------------------------------------------------------------------
tol = max(epsilon(0.0_dp), 1.0d-18)
k1 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1),abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa, 18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa, -41.45_dp)
rl = 1.2_dp * dig + 3.0_dp
fnul = 10.0_dp + 6.0_dp * (dig - 3.0_dp)
az = abs(z)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
IF (az <= aa) THEN
  fn = fnu + (n-1)
  IF (fn <= aa) THEN
    aa = sqrt(aa)
    IF (az > aa) ierr = 3
    IF (fn > aa) ierr = 3
    zn = z
    csgn = cone
    IF (xx < 0.0_dp) THEN
      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_dp) arg = -arg
      s1 = cos(arg)
      s2 = sin(arg)
      csgn = cmplx(s1, s2, kind=dp)
      IF (mod(inu,2) == 1) csgn = -csgn
    END IF
    CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
    IF (nz >= 0) THEN
      IF (xx >= 0.0_dp) RETURN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
!-----------------------------------------------------------------------
      nn = n - nz
      IF (nn == 0) RETURN
      rtol = 1.0_dp / tol
      ascle = tiny(0.0_dp) * rtol * 1.0e+3
      DO  i = 1, nn
!       CY(I) = CY(I)*CSGN
        zn = cy(i)
        aa = real(zn, kind=dp)
        bb = aimag(zn)
        atol = 1.0_dp
        IF (max(abs(aa),abs(bb)) <= ascle) THEN
          zn = zn * rtol
          atol = tol
        END IF
        zn = zn * csgn
        cy(i) = zn * atol
        csgn = -csgn
      END DO
      RETURN
    END IF
    IF (nz /= -2) THEN
      nz = 0
      ierr = 2
      RETURN
    END IF
    nz = 0
    ierr = 5
    RETURN
  END IF
END IF
nz = 0
ierr = 4
RETURN
END SUBROUTINE cbesi
 
 
 
SUBROUTINE cbesj(z, fnu, kode, n, cy, nz, ierr)
!***BEGIN PROLOGUE  CBESJ
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!    ON KODE=1, CBESJ COMPUTES AN N MEMBER  SEQUENCE OF COMPLEX
!    BESSEL FUNCTIONS CY(I) = J(FNU+I-1,Z) FOR REAL (dp), NONNEGATIVE
!    ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
!    -PI < ARG(Z) <= 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
!      Z      - Z=CMPLX(X,Y),  -PI < ARG(Z) <= PI
!      FNU    - ORDER OF INITIAL J FUNCTION, FNU >= 0.0_dp
!      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      - NUMBER OF MEMBERS OF THE SEQUENCE, N >= 1
 
!    OUTPUT
!      CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
!               VALUES 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 > 0 , LAST NZ COMPONENTS OF CY SET TO ZERO
!                         DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0),
!                         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, ABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                       BUT LOSSES OF SIGNIFICANCE BY ARGUMENT
!                       REDUCTION PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!               IERR=4, ABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTATION BECAUSE
!                       OF COMPLETE LOSSES OF SIGNIFICANCE 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) >= 0.0
 
!    J(FNU,Z)=EXP(-FNU*PI*I/2)*I(FNU, I*Z)    AIMAG(Z) < 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 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 > ABS(Z).
 
!    IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELEMENTARY
!    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 = EPSILON(0.0_dp) = UNIT
!    ROUNDOFF.  ALSO 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 = HUGE(0).
!    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.8D+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 REPRESENTS THE INCREASE IN ERROR DUE
!    TO ARGUMENT REDUCTION IN THE ELEMENTARY FUNCTIONS.  HERE,
!    S = MAX(1,ABS(LOG10(ABS(Z))), ABS(LOG10(FNU))) APPROXIMATELY
!    (I.E. S = MAX(1, ABS(EXPONENT OF ABS(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, ACM TRANS. MATH. SOFTWARE,
!         VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CBINU,I1MACH,R1MACH
!***END PROLOGUE  CBESJ
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: ci, csgn, zn
real(dp)     :: aa, alim, arg, dig, elim, fnul, rl, r1, r1m5, r2,  &
                 tol, yy, az, fn, bb, ascle, rtol, atol
INTEGER       :: i, inu, inuh, ir, k1, k2, nl, k
 
real(dp), PARAMETER  :: hpi = 1.570796326794896619_dp
 
!***FIRST EXECUTABLE STATEMENT  CBESJ
ierr = 0
nz = 0
IF (fnu < 0.0_dp) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (n < 1) ierr = 1
IF (ierr /= 0) RETURN
!-----------------------------------------------------------------------
!     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.
!-----------------------------------------------------------------------
tol = max(epsilon(0.0_dp), 1.0d-18)
k1 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1),abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa, 18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa, -41.45_dp)
rl = 1.2_dp * dig + 3.0_dp
fnul = 10.0_dp + 6.0_dp * (dig - 3.0_dp)
ci = cmplx(0.0_dp, 1.0_dp, kind=dp)
yy = aimag(z)
az = abs(z)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
fn = fnu + (n-1)
IF (az <= aa) THEN
  IF (fn <= aa) THEN
    aa = 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
!-----------------------------------------------------------------------
    inu = int(fnu)
    inuh = inu / 2
    ir = inu - 2 * inuh
    arg = (fnu - (inu-ir)) * hpi
    r1 = cos(arg)
    r2 = sin(arg)
    csgn = cmplx(r1, r2, kind=dp)
    IF (mod(inuh,2) == 1) csgn = -csgn
!-----------------------------------------------------------------------
!     ZN IS IN THE RIGHT HALF PLANE
!-----------------------------------------------------------------------
    zn = -z * ci
    IF (yy < 0.0_dp) THEN
      zn = -zn
      csgn = conjg(csgn)
      ci = conjg(ci)
    END IF
    CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
    IF (nz >= 0) THEN
      nl = n - nz
      IF (nl == 0) RETURN
      rtol = 1.0_dp / tol
      ascle = tiny(0.0_dp) * rtol * 1.0e+3
      DO  i = 1, nl
!       CY(I)=CY(I)*CSGN
        zn = cy(i)
        aa = real(zn, kind=dp)
        bb = aimag(zn)
        atol = 1.0_dp
        IF (max(abs(aa),abs(bb)) <= ascle) THEN
          zn = zn * rtol
          atol = tol
        END IF
        zn = zn * csgn
        cy(i) = zn * atol
        csgn = csgn * ci
      END DO
      RETURN
    END IF
    IF (nz /= -2) THEN
      nz = 0
      ierr = 2
      RETURN
    END IF
    nz = 0
    ierr = 5
    RETURN
  END IF
END IF
nz = 0
ierr = 4
RETURN
END SUBROUTINE cbesj
 
 
 
SUBROUTINE cbesk(z, fnu, kode, n, cy, nz, ierr)
!***BEGIN PROLOGUE  CBESK
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!   ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX BESSEL FUNCTIONS
!   CY(J)=K(FNU+J-1,Z) FOR REAL (dp), NONNEGATIVE ORDERS FNU+J-1, J=1,...,N
!   AND COMPLEX Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI < ARG(Z) <= 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
!     Z      - Z=CMPLX(X,Y),Z.NE.CMPLX(0.,0.),-PI < ARG(Z) <= PI
!     FNU    - ORDER OF INITIAL K FUNCTION, FNU >= 0.0_dp
!     N      - NUMBER OF MEMBERS OF THE SEQUENCE, N >= 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
!     CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
!              VALUES 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 > 0 , FIRST NZ COMPONENTS OF CY SET TO ZERO
!                        DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0),
!                        I=1,...,N WHEN X >= 0.0.  WHEN X < 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+N-1 IS
!                      TOO LARGE OR ABS(Z) IS TOO SMALL OR BOTH
!              IERR=3, ABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE, BUT
!                      LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION PRODUCE
!                      LESS THAN HALF OF MACHINE ACCURACY
!              IERR=4, ABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTATION BECAUSE OF
!                      COMPLETE LOSSES OF SIGNIFICANCE 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 >= 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) > 0, I**2=-1
 
!   WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
 
!   FOR LARGE ORDERS, FNU > 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 ELEMENTARY
!   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 = EPSILON(0.0_dp) = UNIT ROUNDOFF.  ALSO 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=HUGE(0).  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.8D+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 REPRESENTS THE INCREASE IN ERROR DUE TO
!   ARGUMENT REDUCTION IN THE ELEMENTARY FUNCTIONS.  HERE, S =
!   MAX(1,ABS(LOG10(ABS(Z))), ABS(LOG10(FNU))) APPROXIMATELY (I.E. S =
!   MAX(1,ABS(EXPONENT OF ABS(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, ACM TRANS. MATH. SOFTWARE,
!           VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH
!***END PROLOGUE  CBESK
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
real(dp)  :: aa, alim, aln, arg, az, dig, elim, fn, fnul, rl, r1m5,  &
              tol, ufl, xx, yy, bb
INTEGER    :: k, k1, k2, mr, nn, nuf, nw
 
!***FIRST EXECUTABLE STATEMENT  CBESK
ierr = 0
nz = 0
xx = real(z, kind=dp)
yy = aimag(z)
IF (yy == 0.0_dp .AND. xx == 0.0_dp) ierr = 1
IF (fnu < 0.0_dp) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (n < 1) ierr = 1
IF (ierr /= 0) RETURN
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) < 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
!-----------------------------------------------------------------------
tol = max(epsilon(0.0_dp), 1.0d-18)
k1 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1),abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa, 18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa, -41.45_dp)
fnul = 10.0_dp + 6.0_dp * (dig - 3.0_dp)
rl = 1.2_dp * dig + 3.0_dp
az = abs(z)
fn = fnu + (nn-1)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
IF (az <= aa) THEN
  IF (fn <= aa) THEN
    aa = sqrt(aa)
    IF (az > aa) ierr = 3
    IF (fn > aa) ierr = 3
!-----------------------------------------------------------------------
!     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
!-----------------------------------------------------------------------
!     UFL = EXP(-ELIM)
    ufl = tiny(0.0_dp) * 1.0e+3
    IF (az >= ufl) THEN
      IF (fnu <= fnul) THEN
        IF (fn > 1.0_dp) THEN
          IF (fn <= 2.0_dp) THEN
            IF (az > tol) GO TO 10
            arg = 0.5_dp * az
            aln = -fn * log(arg)
            IF (aln > elim) GO TO 30
          ELSE
            CALL cuoik(z, fnu, kode, 2, nn, cy, nuf, tol, elim, alim)
            IF (nuf < 0) GO TO 30
            nz = nz + nuf
            nn = 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 20
          END IF
        END IF
 
        10 IF (xx >= 0.0_dp) THEN
!-----------------------------------------------------------------------
!     RIGHT HALF PLANE COMPUTATION, REAL(Z) >= 0.
!-----------------------------------------------------------------------
          CALL cbknu(z, fnu, kode, nn, cy, nw, tol, elim, alim)
          IF (nw < 0) GO TO 40
          nz = nw
          RETURN
        END IF
!-----------------------------------------------------------------------
!     LEFT HALF PLANE COMPUTATION
!     PI/2 < ARG(Z) <= PI AND -PI < ARG(Z) < -PI/2.
!-----------------------------------------------------------------------
        IF (nz /= 0) GO TO 30
        mr = 1
        IF (yy < 0.0_dp) mr = -1
        CALL cacon(z, fnu, kode, mr, nn, cy, nw, rl, fnul, tol, elim, alim)
        IF (nw < 0) GO TO 40
        nz = nw
        RETURN
      END IF
!-----------------------------------------------------------------------
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU > FNUL
!-----------------------------------------------------------------------
      mr = 0
      IF (xx < 0.0_dp) THEN
        mr = 1
        IF (yy < 0.0_dp) mr = -1
      END IF
      CALL cbunk(z, fnu, kode, mr, nn, cy, nw, tol, elim, alim)
      IF (nw < 0) GO TO 40
      nz = nz + nw
      RETURN
 
      20 IF (xx >= 0.0_dp) RETURN
    END IF
 
    30 nz = 0
    ierr = 2
    RETURN
 
    40 IF (nw == -1) GO TO 30
    nz = 0
    ierr = 5
    RETURN
  END IF
END IF
nz = 0
ierr = 4
RETURN
END SUBROUTINE cbesk
 
 
 
SUBROUTINE cbesy(z, fnu, kode, n, cy, nz, ierr)
 
! N.B. Argument CWRK has been removed.
 
!***BEGIN PROLOGUE  CBESY
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101  (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
 
!   ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
!   BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL (dp), NONNEGATIVE
!   ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
!   -PI < ARG(Z) <= 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
!     Z      - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI < ARG(Z) <= PI
!     FNU    - ORDER OF INITIAL Y FUNCTION, FNU >= 0.0_dp
!     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 >= 1
!     CWRK   - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N
 
!   OUTPUT
!     CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
!              VALUES 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 > 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+N-1 IS
!                      TOO LARGE OR ABS(Z) IS TOO SMALL OR BOTH
!              IERR=3, ABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
!                      BUT LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION
!                      PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!              IERR=4, ABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTATION
!                      BECAUSE OF COMPLETE LOSSES OF SIGNIFICANCE
!                      BY ARGUMENT REDUCTION
!              IERR=5, ERROR              - NO COMPUTATION,
!                      ALGORITHM TERMINATION CONDITION NOT MET
 
!***LONG DESCRIPTION
 
!   THE COMPUTATION IS CARRIED OUT IN TERMS OF THE I(FNU,Z) AND
!   K(FNU,Z) BESSEL FUNCTIONS IN THE RIGHT HALF PLANE BY
 
!       Y(FNU,Z) = I*CC*I(FNU,ARG) - (2/PI)*CONJG(CC)*K(FNU,ARG)
 
!       Y(FNU,Z) = CONJG(Y(FNU,CONJG(Z)))
 
!   FOR AIMAG(Z) >= 0 AND AIMAG(Z) < 0 RESPECTIVELY, WHERE
!   CC=EXP(I*PI*FNU/2), ARG=Z*EXP(-I*PI/2) AND I**2=-1.
 
!   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 > ABS(Z).
 
!   IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELEMENTARY
!   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 = EPSILON(0.0_dp) = UNIT
!   ROUNDOFF.  ALSO 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 = HUGE(0).
!   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.8D+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 REPRESENTS THE INCREASE IN ERROR DUE TO
!   ARGUMENT REDUCTION IN THE ELEMENTARY FUNCTIONS.  HERE, S =
!   MAX(1,ABS(LOG10(ABS(Z))), ABS(LOG10(FNU))) APPROXIMATELY (I.E. S =
!   MAX(1,ABS(EXPONENT OF ABS(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, ACM TRANS. MATH. SOFTWARE,
!        VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CBESI,CBESK,I1MACH,R1MACH
!***END PROLOGUE  CBESY
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: ci, csgn, cspn, cwrk(n), ex, zu, zv, zz, zn
real(dp)     :: arg, elim, ey, r1, r2, tay, xx, yy, ascle, rtol,  &
                 atol, tol, aa, bb, ffnu, rhpi, r1m5
INTEGER       :: i, ifnu, k, k1, k2, nz1, nz2, i4
COMPLEX (dp), PARAMETER  :: cip(4) = (/ (1.0_dp, 0.0_dp), (0.0_dp, 1.0_dp),  &
                                        (-1.0_dp, 0.0_dp), (0.0_dp, -1.0_dp) /)
real(dp), PARAMETER     :: hpi = 1.57079632679489662_dp
 
!***FIRST EXECUTABLE STATEMENT  CBESY
xx = real(z, kind=dp)
yy = aimag(z)
ierr = 0
nz = 0
IF (xx == 0.0_dp .AND. yy == 0.0_dp) ierr = 1
IF (fnu < 0.0_dp) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (n < 1) ierr = 1
IF (ierr /= 0) RETURN
ci = cmplx(0.0_dp, 1.0_dp, kind=dp)
zz = z
IF (yy < 0.0_dp) zz = conjg(z)
zn = -ci * zz
CALL cbesi(zn, fnu, kode, n, cy, nz1, ierr)
IF (ierr == 0 .OR. ierr == 3) THEN
  CALL cbesk(zn, fnu, kode, n, cwrk, nz2, ierr)
  IF (ierr == 0 .OR. ierr == 3) THEN
    nz = min(nz1, nz2)
    ifnu = int(fnu)
    ffnu = fnu - ifnu
    arg = hpi * ffnu
    csgn = cmplx(cos(arg), sin(arg), kind=dp)
    i4 = mod(ifnu, 4) + 1
    csgn = csgn * cip(i4)
    rhpi = 1.0_dp / hpi
    cspn = conjg(csgn) * rhpi
    csgn = csgn * ci
    IF (kode /= 2) THEN
      DO  i = 1, n
        cy(i) = csgn * cy(i) - cspn * cwrk(i)
        csgn = ci * csgn
        cspn = -ci * cspn
      END DO
      IF (yy < 0.0_dp) cy(1:n) = conjg(cy(1:n))
      RETURN
    END IF
 
    r1 = cos(xx)
    r2 = sin(xx)
    ex = cmplx(r1, r2, kind=dp)
    tol = max(epsilon(0.0_dp), 1.0d-18)
    k1 = minexponent(0.0_dp)
    k2 = maxexponent(0.0_dp)
    k = min(abs(k1),abs(k2))
    r1m5 = log10( real( radix(0.0_dp), kind=dp) )
!-----------------------------------------------------------------------
!     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
!-----------------------------------------------------------------------
    elim = 2.303_dp * (k*r1m5 - 3.0_dp)
    ey = 0.0_dp
    tay = abs(yy+yy)
    IF (tay < elim) ey = exp(-tay)
    cspn = ex * ey * cspn
    nz = 0
    rtol = 1.0_dp / tol
    ascle = tiny(0.0_dp) * rtol * 1.0e+3
    DO  i = 1, n
!----------------------------------------------------------------------
!       CY(I) = CSGN*CY(I)-CSPN*CWRK(I): PRODUCTS ARE COMPUTED IN
!       SCALED MODE IF CY(I) OR CWRK(I) ARE CLOSE TO UNDERFLOW TO
!       PREVENT UNDERFLOW IN AN INTERMEDIATE COMPUTATION.
!----------------------------------------------------------------------
      zv = cwrk(i)
      aa = real(zv, kind=dp)
      bb = aimag(zv)
      atol = 1.0_dp
      IF (max(abs(aa),abs(bb)) <= ascle) THEN
        zv = zv * rtol
        atol = tol
      END IF
      zv = zv * cspn
      zv = zv * atol
      zu = cy(i)
      aa = real(zu, kind=dp)
      bb = aimag(zu)
      atol = 1.0_dp
      IF (max(abs(aa),abs(bb)) <= ascle) THEN
        zu = zu * rtol
        atol = tol
      END IF
      zu = zu * csgn
      zu = zu * atol
      cy(i) = zu - zv
      IF (yy < 0.0_dp) cy(i) = conjg(cy(i))
      IF (cy(i) == cmplx(0.0_dp, 0.0_dp, kind=dp) .AND. ey == 0.0_dp) nz = nz + 1
      csgn = ci * csgn
      cspn = -ci * cspn
    END DO
    RETURN
  END IF
END IF
nz = 0
RETURN
END SUBROUTINE cbesy
 
 
 
SUBROUTINE cairy(z, id, kode, ai, nz, ierr)
!***BEGIN PROLOGUE  CAIRY
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!   ON KODE=1, CAIRY 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 EXP(ZTA)*AI(Z) OR EXP(ZTA)*
!   DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
!   -PI/3 < ARG(Z) < PI/3 AND THE EXPONENTIAL GROWTH IN
!   PI/3 < ABS(ARG(Z)) < PI WHERE ZTA=(2/3)*Z*SQRT(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
!     Z      - Z=CMPLX(X,Y)
!     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=EXP(ZTA)*AI(Z)       ON ID=0 OR
!                       AI=EXP(ZTA)*DAI(Z)/DZ   ON ID=1 WHERE
!                       ZTA=(2/3)*Z*SQRT(Z)
 
!   OUTPUT
!     AI     - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND KODE
!     NZ     - UNDERFLOW INDICATOR
!              NZ= 0   , NORMAL RETURN
!              NZ= 1   , AI=CMPLX(0.0,0.0) DUE TO UNDERFLOW IN
!                        -PI/3 < ARG(Z) < 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 WITH KODE=1.
!              IERR=3, ABS(Z) LARGE      - COMPUTATION COMPLETED
!                      LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION
!                      PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!              IERR=4, ABS(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 ABS(Z) > 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 ABS(Z) <= 1.0.
 
!   IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELEMENTARY
!   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 =
!   EPSILON(0.0_dp) = UNIT ROUNDOFF.  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 = HUGE(0).  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.8D+16, 2.1E+9 IN DOUBLE PRECISION ARITHMETIC RESPECTIVELY.
!   THIS MAKES U2 AND U3 LIMITING IN THEIR RESPECTIVE ARITHMETICS.  THIS MEANS
!   THAT THE MAGNITUDE 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 REPRESENTS
!   THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
!   ELEMENTARY FUNCTIONS.  HERE, S = MAX(1,ABS(LOG10(ABS(Z))),
!   ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
!   ABS(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, ACM TRANS. MATH. SOFTWARE,
!      VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CACAI,CBKNU,I1MACH,R1MACH
!***END PROLOGUE  CAIRY
 
COMPLEX (dp), INTENT(IN)   :: z
INTEGER, INTENT(IN)        :: id
INTEGER, INTENT(IN)        :: kode
COMPLEX (dp), INTENT(OUT)  :: ai
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: csq, cy(1), s1, s2, trm1, trm2, zta, z3
real(dp)     :: aa, ad, ak, alim, atrm, az, az3, bk, ck, dig, dk, d1, d2, &
                 elim, fid, fnu, rl, r1m5, sfac, tol, zi, zr, z3i, z3r, bb, &
                 alaz
INTEGER       :: iflag, k, k1, k2, mr, nn
real(dp), PARAMETER  :: tth = 6.66666666666666667d-01,  &
           c1 = 3.55028053887817240d-01, c2 = 2.58819403792806799d-01,  &
           coef = 1.83776298473930683d-01
COMPLEX (dp), PARAMETER  :: cone = (1.0_dp, 0.0_dp)
 
!***FIRST EXECUTABLE STATEMENT  CAIRY
ierr = 0
nz = 0
IF (id < 0 .OR. id > 1) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (ierr /= 0) RETURN
az = abs(z)
tol = max(epsilon(0.0_dp), 1.0d-18)
fid = id
IF (az <= 1.0_dp) THEN
!-----------------------------------------------------------------------
!     POWER SERIES FOR ABS(Z) <= 1.
!-----------------------------------------------------------------------
  s1 = cone
  s2 = cone
  IF (az < tol) GO TO 30
  aa = az * az
  IF (aa >= tol/az) THEN
    trm1 = cone
    trm2 = cone
    atrm = 1.0_dp
    z3 = z * z * z
    az3 = az * aa
    ak = 2.0_dp + fid
    bk = 3.0_dp - fid - fid
    ck = 4.0_dp - fid
    dk = 3.0_dp + fid + fid
    d1 = ak * dk
    d2 = bk * ck
    ad = min(d1,d2)
    ak = 24.0_dp + 9.0_dp * fid
    bk = 30.0_dp - 9.0_dp * fid
    z3r = real(z3, kind=dp)
    z3i = aimag(z3)
    DO  k = 1, 25
      trm1 = trm1 * cmplx(z3r/d1, z3i/d1, kind=dp)
      s1 = s1 + trm1
      trm2 = trm2 * cmplx(z3r/d2, z3i/d2, kind=dp)
      s2 = s2 + trm2
      atrm = atrm * az3 / ad
      d1 = d1 + ak
      d2 = d2 + bk
      ad = min(d1,d2)
      IF (atrm < tol*ad) EXIT
      ak = ak + 18.0_dp
      bk = bk + 18.0_dp
    END DO
  END IF
 
  IF (id /= 1) THEN
    ai = s1 * c1 - z * s2 * c2
    IF (kode == 1) RETURN
    zta = z * sqrt(z) * tth
    ai = ai * exp(zta)
    RETURN
  END IF
  ai = -s2 * c2
  IF (az > tol) ai = ai + z * z * s1 * c1/(1.0_dp + fid)
  IF (kode == 1) RETURN
  zta = z * sqrt(z) * tth
  ai = ai * exp(zta)
  RETURN
END IF
!-----------------------------------------------------------------------
!     CASE FOR ABS(Z) > 1.0
!-----------------------------------------------------------------------
fnu = (1.0_dp + fid) / 3.0_dp
!-----------------------------------------------------------------------
!     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 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1),abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa,18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa,-41.45_dp)
rl = 1.2_dp * dig + 3.0_dp
alaz = log(az)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
aa = aa ** tth
IF (az > aa) GO TO 70
aa = sqrt(aa)
IF (az > aa) ierr = 3
csq = sqrt(z)
zta = z * csq * tth
!-----------------------------------------------------------------------
!     RE(ZTA) <= 0 WHEN RE(Z) < 0, ESPECIALLY WHEN IM(Z) IS SMALL
!-----------------------------------------------------------------------
iflag = 0
sfac = 1.0_dp
zi = aimag(z)
zr = real(z, kind=dp)
ak = aimag(zta)
IF (zr < 0.0_dp) THEN
  bk = real(zta, kind=dp)
  ck = -abs(bk)
  zta = cmplx(ck, ak, kind=dp)
END IF
IF (zi == 0.0_dp) THEN
  IF (zr <= 0.0_dp) THEN
    zta = cmplx(0.0_dp, ak, kind=dp)
  END IF
END IF
aa = real(zta, kind=dp)
IF (aa < 0.0_dp .OR. zr <= 0.0_dp) THEN
  IF (kode /= 2) THEN
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
    IF (aa <= -alim) THEN
      aa = -aa + 0.25_dp * alaz
      iflag = 1
      sfac = tol
      IF (aa > elim) GO TO 50
    END IF
  END IF
!-----------------------------------------------------------------------
!     CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
!-----------------------------------------------------------------------
  mr = 1
  IF (zi < 0.0_dp) mr = -1
  CALL cacai(zta, fnu, kode, mr, 1, cy, nn, rl, tol, elim, alim)
  IF (nn < 0) GO TO 60
  nz = nz + nn
ELSE
  IF (kode /= 2) THEN
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
    IF (aa >= alim) THEN
      aa = -aa - 0.25_dp * alaz
      iflag = 2
      sfac = 1.0_dp / tol
      IF (aa < -elim) GO TO 40
    END IF
  END IF
  CALL cbknu(zta, fnu, kode, 1, cy, nz, tol, elim, alim)
END IF
s1 = cy(1) * coef
IF (iflag == 0) THEN
  IF (id /= 1) THEN
    ai = csq * s1
    RETURN
  END IF
  ai = -z * s1
  RETURN
END IF
s1 = s1 * sfac
IF (id /= 1) THEN
  s1 = s1 * csq
  ai = s1 / sfac
  RETURN
END IF
s1 = -s1 * z
ai = s1 / sfac
RETURN
 
30 aa = 1.0e+3 * tiny(0.0_dp)
s1 = cmplx(0.0_dp, 0.0_dp, kind=dp)
IF (id /= 1) THEN
  IF (az > aa) s1 = c2 * z
  ai = c1 - s1
  RETURN
END IF
ai = -c2
aa = sqrt(aa)
IF (az > aa) s1 = z * z * 0.5_dp
ai = ai + s1 * c1
RETURN
 
40 nz = 1
ai = cmplx(0.0_dp, 0.0_dp, kind=dp)
RETURN
 
50 nz = 0
ierr = 2
RETURN
 
60 IF (nn == -1) GO TO 50
nz = 0
ierr = 5
RETURN
 
70 ierr = 4
nz = 0
RETURN
END SUBROUTINE cairy
 
 
 
SUBROUTINE cbiry(z, id, kode, bi, ierr)
!***BEGIN PROLOGUE  CBIRY
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  890801, 930101   (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
 
!   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 EXP(-AXZTA)*BI(Z) OR EXP(-AXZTA)*DBI(Z)/DZ
!   IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND
!   RIGHT HALF PLANES WHERE ZTA = (2/3)*Z*SQRT(Z) = CMPLX(XZTA,YZTA)
!   AND AXZTA=ABS(XZTA).
!   DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
!   FUNCTIONS (REF. 1).
 
!   INPUT
!     Z      - Z=CMPLX(X,Y)
!     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=EXP(-AXZTA)*BI(Z)     ON ID=0 OR
!                       BI=EXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE
!                       ZTA=(2/3)*Z*SQRT(Z)=CMPLX(XZTA,YZTA)
!                       AND AXZTA=ABS(XZTA)
 
!   OUTPUT
!     BI     - 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 WITH KODE=1
!              IERR=3, ABS(Z) LARGE      - COMPUTATION COMPLETED
!                      LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION
!                      PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!              IERR=4, ABS(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 ABS(Z) > 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 ABS(Z) <= 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=EPSILON(0.0_dp)=UNIT ROUNDOFF.
!      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=HUGE(0). 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.8D+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.
 
!      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(ABS(Z))),
!      ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
!      ABS(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, ACM TRANS. MATH. SOFTWARE,
!           VOL. 12, NO. 3, SEPTEMBER 1986, PP 265-273.
 
!***ROUTINES CALLED  CBINU,I1MACH,R1MACH
!***END PROLOGUE  CBIRY
 
COMPLEX (dp), INTENT(IN)   :: z
INTEGER, INTENT(IN)        :: id
INTEGER, INTENT(IN)        :: kode
COMPLEX (dp), INTENT(OUT)  :: bi
INTEGER, INTENT(OUT)       :: ierr
 
COMPLEX (dp)  :: csq, cy(2), s1, s2, trm1, trm2, zta, z3
real(dp)     :: 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, z3i, z3r
INTEGER       :: k, k1, k2, nz
real(dp), PARAMETER  :: tth = 6.66666666666666667d-01,  &
    c1 = 6.14926627446000736d-01, c2 = 4.48288357353826359d-01,  &
    coef = 5.77350269189625765d-01, pi = 3.141592653589793238_dp
real(dp), PARAMETER  :: cone = (1.0_dp,0.0_dp)
 
!***FIRST EXECUTABLE STATEMENT  CBIRY
ierr = 0
nz = 0
IF (id < 0 .OR. id > 1) ierr = 1
IF (kode < 1 .OR. kode > 2) ierr = 1
IF (ierr /= 0) RETURN
az = abs(z)
tol = max(epsilon(0.0_dp), 1.0d-18)
fid = id
IF (az <= 1.0_dp) THEN
!-----------------------------------------------------------------------
!     POWER SERIES FOR ABS(Z) <= 1.
!-----------------------------------------------------------------------
  s1 = cone
  s2 = cone
  IF (az < tol) GO TO 30
  aa = az * az
  IF (aa >= tol/az) THEN
    trm1 = cone
    trm2 = cone
    atrm = 1.0_dp
    z3 = z * z * z
    az3 = az * aa
    ak = 2.0_dp + fid
    bk = 3.0_dp - fid - fid
    ck = 4.0_dp - fid
    dk = 3.0_dp + fid + fid
    d1 = ak * dk
    d2 = bk * ck
    ad = min(d1,d2)
    ak = 24.0_dp + 9.0_dp * fid
    bk = 30.0_dp - 9.0_dp * fid
    z3r = real(z3, kind=dp)
    z3i = aimag(z3)
    DO  k = 1, 25
      trm1 = trm1 * cmplx(z3r/d1, z3i/d1, kind=dp)
      s1 = s1 + trm1
      trm2 = trm2 * cmplx(z3r/d2, z3i/d2, kind=dp)
      s2 = s2 + trm2
      atrm = atrm * az3 / ad
      d1 = d1 + ak
      d2 = d2 + bk
      ad = min(d1,d2)
      IF (atrm < tol*ad) EXIT
      ak = ak + 18.0_dp
      bk = bk + 18.0_dp
    END DO
  END IF
 
  IF (id /= 1) THEN
    bi = s1 * c1 + z * s2 * c2
    IF (kode == 1) RETURN
    zta = z * sqrt(z) * tth
    aa = real(zta, kind=dp)
    aa = -abs(aa)
    bi = bi * exp(aa)
    RETURN
  END IF
  bi = s2 * c2
  IF (az > tol) bi = bi + z * z * s1 * c1/(1.0_dp+fid)
  IF (kode == 1) RETURN
  zta = z * sqrt(z) * tth
  aa = real(zta, kind=dp)
  aa = -abs(aa)
  bi = bi * exp(aa)
  RETURN
END IF
!-----------------------------------------------------------------------
!     CASE FOR ABS(Z) > 1.0
!-----------------------------------------------------------------------
fnu = (1.0_dp+fid) / 3.0_dp
!-----------------------------------------------------------------------
!     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 = minexponent(0.0_dp)
k2 = maxexponent(0.0_dp)
r1m5 = log10( real( radix(0.0_dp), kind=dp) )
k = min(abs(k1),abs(k2))
elim = 2.303_dp * (k*r1m5 - 3.0_dp)
k1 = digits(0.0_dp) - 1
aa = r1m5 * k1
dig = min(aa,18.0_dp)
aa = aa * 2.303_dp
alim = elim + max(-aa,-41.45_dp)
rl = 1.2_dp * dig + 3.0_dp
fnul = 10.0_dp + 6.0_dp * (dig - 3.0_dp)
!-----------------------------------------------------------------------
!     TEST FOR RANGE
!-----------------------------------------------------------------------
aa = 0.5_dp / tol
bb = huge(0) * 0.5_dp
aa = min(aa,bb)
aa = aa ** tth
IF (az > aa) GO TO 60
aa = sqrt(aa)
IF (az > aa) ierr = 3
csq = sqrt(z)
zta = z * csq * tth
!-----------------------------------------------------------------------
!     RE(ZTA) <= 0 WHEN RE(Z) < 0, ESPECIALLY WHEN IM(Z) IS SMALL
!-----------------------------------------------------------------------
sfac = 1.0_dp
zi = aimag(z)
zr = real(z, kind=dp)
ak = aimag(zta)
IF (zr < 0.0_dp) THEN
  bk = real(zta, kind=dp)
  ck = -abs(bk)
  zta = cmplx(ck, ak, kind=dp)
END IF
IF (zi == 0.0_dp .AND. zr <= 0.0_dp) zta = cmplx(0.0_dp, ak, kind=dp)
aa = real(zta, kind=dp)
IF (kode /= 2) THEN
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
  bb = abs(aa)
  IF (bb >= alim) THEN
    bb = bb + 0.25_dp * log(az)
    sfac = tol
    IF (bb > elim) GO TO 40
  END IF
END IF
fmr = 0.0_dp
IF (aa < 0.0_dp .OR. zr <= 0.0_dp) THEN
  fmr = pi
  IF (zi < 0.0_dp) fmr = -pi
  zta = -zta
END IF
!-----------------------------------------------------------------------
!     AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA)
!     KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBINU
!-----------------------------------------------------------------------
CALL cbinu(zta,fnu,kode,1,cy,nz,rl,fnul,tol,elim,alim)
IF (nz < 0) GO TO 50
aa = fmr * fnu
z3 = cmplx(sfac, 0.0_dp, kind=dp)
s1 = cy(1) * cmplx(cos(aa), sin(aa), kind=dp) * z3
fnu = (2.0_dp - fid) / 3.0_dp
CALL cbinu(zta, fnu, kode, 2, cy, nz, rl, fnul, tol, elim, alim)
cy(1) = cy(1) * z3
cy(2) = cy(2) * z3
!-----------------------------------------------------------------------
!     BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3
!-----------------------------------------------------------------------
s2 = cy(1) * cmplx(fnu+fnu, 0.0_dp, kind=dp) / zta + cy(2)
aa = fmr * (fnu-1.0_dp)
s1 = (s1 + s2*cmplx(cos(aa), sin(aa), kind=dp)) * coef
IF (id /= 1) THEN
  s1 = csq * s1
  bi = s1 / sfac
  RETURN
END IF
s1 = z * s1
bi = s1 / sfac
RETURN
 
30 aa = c1 * (1.0_dp-fid) + fid * c2
bi = cmplx(aa, 0.0_dp, kind=dp)
RETURN
 
40 nz = 0
ierr = 2
RETURN
 
50 IF (nz == -1) GO TO 40
nz = 0
ierr = 5
RETURN
 
60 ierr = 4
nz = 0
RETURN
END SUBROUTINE cbiry
 
 
 
SUBROUTINE cunik(zr, fnu, ikflg, ipmtr, tol, init, phi, zeta1, zeta2, total, &
                 cwrk)
 
!***BEGIN PROLOGUE  CUNIK
!***REFER TO  CBESI,CBESK
 
!  CUNIK 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  (NONE)
!***END PROLOGUE  CUNIK
 
COMPLEX (dp), INTENT(IN)      :: zr
real(dp), INTENT(IN)         :: fnu
INTEGER, INTENT(IN)           :: ikflg
INTEGER, INTENT(IN)           :: ipmtr
real(dp), INTENT(IN)         :: tol
INTEGER, INTENT(IN OUT)       :: init
COMPLEX (dp), INTENT(OUT)     :: phi
COMPLEX (dp), INTENT(IN OUT)  :: zeta1
COMPLEX (dp), INTENT(IN OUT)  :: zeta2
COMPLEX (dp), INTENT(IN OUT)  :: total
COMPLEX (dp), INTENT(IN OUT)  :: cwrk(16)
 
COMPLEX (dp)  :: cfn, crfn, s, sr, t, t2,  zn
real(dp)     :: ac, rfn, test, tstr, tsti
INTEGER       :: i, j, k, l
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp, 0.0_dp)
real(dp), PARAMETER     :: con(2) = (/ 3.98942280401432678d-01, &
                                        1.25331413731550025_dp /)
real(dp), PARAMETER  :: c(120) = (/  &
  1.00000000000000000_dp, -2.08333333333333333d-01, 1.25000000000000000d-01, 3.34201388888888889d-01, &
 -4.01041666666666667d-01, 7.03125000000000000d-02, -1.02581259645061728_dp, 1.84646267361111111_dp,  &
 -8.91210937500000000d-01, 7.32421875000000000d-02, 4.66958442342624743_dp, -1.12070026162229938d+01,  &
  8.78912353515625000_dp, -2.36408691406250000_dp, 1.12152099609375000d-01, -2.82120725582002449d+01,  &
  8.46362176746007346d+01, -9.18182415432400174d+01, 4.25349987453884549d+01, -7.36879435947963170_dp,  &
  2.27108001708984375d-01, 2.12570130039217123d+02, -7.65252468141181642d+02, 1.05999045252799988d+03, &
 -6.99579627376132541d+02, 2.18190511744211590d+02, -2.64914304869515555d+01, 5.72501420974731445d-01, &
 -1.91945766231840700d+03, 8.06172218173730938d+03, -1.35865500064341374d+04, 1.16553933368645332d+04, &
 -5.30564697861340311d+03, 1.20090291321635246d+03, -1.08090919788394656d+02, 1.72772750258445740_dp,  &
  2.02042913309661486d+04, -9.69805983886375135d+04, 1.92547001232531532d+05, -2.03400177280415534d+05,  &
  1.22200464983017460d+05, -4.11926549688975513d+04, 7.10951430248936372d+03, -4.93915304773088012d+02,  &
  6.07404200127348304_dp, -2.42919187900551333d+05, 1.31176361466297720d+06, -2.99801591853810675d+06,  &
  3.76327129765640400d+06, -2.81356322658653411d+06, 1.26836527332162478d+06, -3.31645172484563578d+05,  &
  4.52187689813627263d+04, -2.49983048181120962d+03, 2.43805296995560639d+01, 3.28446985307203782d+06,  &
 -1.97068191184322269d+07, 5.09526024926646422d+07, -7.41051482115326577d+07, 6.63445122747290267d+07,  &
 -3.75671766607633513d+07, 1.32887671664218183d+07, -2.78561812808645469d+06, 3.08186404612662398d+05,  &
 -1.38860897537170405d+04, 1.10017140269246738d+02, -4.93292536645099620d+07, 3.25573074185765749d+08,  &
 -9.39462359681578403d+08, 1.55359689957058006d+09, -1.62108055210833708d+09, 1.10684281682301447d+09,  &
 -4.95889784275030309d+08, 1.42062907797533095d+08, -2.44740627257387285d+07, 2.24376817792244943d+06, &
 -8.40054336030240853d+04, 5.51335896122020586d+02, 8.14789096118312115d+08, -5.86648149205184723d+09,  &
  1.86882075092958249d+10, -3.46320433881587779d+10, 4.12801855797539740d+10, -3.30265997498007231d+10,  &
  1.79542137311556001d+10, -6.56329379261928433d+09, 1.55927986487925751d+09, -2.25105661889415278d+08,  &
  1.73951075539781645d+07, -5.49842327572288687d+05, 3.03809051092238427d+03, -1.46792612476956167d+10,  &
  1.14498237732025810d+11, -3.99096175224466498d+11, 8.19218669548577329d+11, -1.09837515608122331d+12,  &
  1.00815810686538209d+12, -6.45364869245376503d+11, 2.87900649906150589d+11, -8.78670721780232657d+10,  &
  1.76347306068349694d+10, -2.16716498322379509d+09, 1.43157876718888981d+08, -3.87183344257261262d+06,  &
  1.82577554742931747d+04, 2.86464035717679043d+11, -2.40629790002850396d+12, 9.10934118523989896d+12,  &
 -2.05168994109344374d+13, 3.05651255199353206d+13, -3.16670885847851584d+13, 2.33483640445818409d+13,  &
 -1.23204913055982872d+13, 4.61272578084913197d+12, -1.19655288019618160d+12, 2.05914503232410016d+11,  &
 -2.18229277575292237d+10, 1.24700929351271032d+09, -2.91883881222208134d+07, 1.18838426256783253d+05 /)
 
IF (init == 0) THEN
!-----------------------------------------------------------------------
!     INITIALIZE ALL VARIABLES
!-----------------------------------------------------------------------
  rfn = 1.0_dp / fnu
  crfn = rfn
  cwrk = czero
!     T = ZR*CRFN
!-----------------------------------------------------------------------
!     OVERFLOW TEST (ZR/FNU TOO SMALL)
!-----------------------------------------------------------------------
  tstr = real(zr, kind=dp)
  tsti = aimag(zr)
  test = tiny(0.0_dp) * 1.0e+3
  ac = fnu * test
  IF (abs(tstr) <= ac .AND. abs(tsti) <= ac) THEN
    ac = 2.0_dp * abs(log(test)) + fnu
    zeta1 = ac
    zeta2 = fnu
    phi = cone
    RETURN
  END IF
  t = zr * crfn
  s = cone + t * t
  sr = sqrt(s)
  cfn = fnu
  zn = (cone+sr) / t
  zeta1 = cfn * log(zn)
  zeta2 = cfn * sr
  t = cone / sr
  sr = t * crfn
  cwrk(16) = sqrt(sr)
  phi = cwrk(16) * con(ikflg)
  IF (ipmtr /= 0) RETURN
  t2 = cone / s
  cwrk(1) = cone
  crfn = cone
  ac = 1.0_dp
  l = 1
  DO  k = 2, 15
    s = czero
    DO  j = 1, k
      l = l + 1
      s = s * t2 + c(l)
    END DO
    crfn = crfn * sr
    cwrk(k) = crfn * s
    ac = ac * rfn
    tstr = real(cwrk(k), kind=dp)
    tsti = aimag(cwrk(k))
    test = abs(tstr) + abs(tsti)
    IF (ac < tol .AND. test < tol) GO TO 30
  END DO
  k = 15
 
  30 init = k
END IF
 
IF (ikflg /= 2) THEN
!-----------------------------------------------------------------------
!     COMPUTE SUM FOR THE I FUNCTION
!-----------------------------------------------------------------------
  total = sum( cwrk(1:init) )
  phi = cwrk(16) * con(1)
  RETURN
END IF
!-----------------------------------------------------------------------
!     COMPUTE SUM FOR THE K FUNCTION
!-----------------------------------------------------------------------
s = czero
t = cone
DO  i = 1, init
  s = s + t * cwrk(i)
  t = -t
END DO
total = s
phi = cwrk(16) * con(2)
RETURN
END SUBROUTINE cunik
 
 
 
SUBROUTINE cuoik(z, fnu, kode, ikflg, n, y, nuf, tol, elim, alim)
!***BEGIN PROLOGUE  CUOIK
!***REFER TO  CBESI,CBESK,CBESH
 
!   CUOIK 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 < 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*1000
!   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 > 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 = N MEANS ALL Y VALUES WERE SET TO ZERO
!   IKFLG=2 AND 0 < NUF < N NOT CONSIDERED.  Y MUST BE SET BY ANOTHER ROUTINE
 
!***ROUTINES CALLED  CUCHK,CUNHJ,CUNIK,R1MACH
!***END PROLOGUE  CUOIK
 
COMPLEX (dp), INTENT(IN)      :: z
real(dp), INTENT(IN)         :: fnu
INTEGER, INTENT(IN)           :: kode
INTEGER, INTENT(IN)           :: ikflg
INTEGER, INTENT(IN)           :: n
COMPLEX (dp), INTENT(IN OUT)  :: y(n)
INTEGER, INTENT(OUT)          :: nuf
real(dp), INTENT(IN)         :: tol
real(dp), INTENT(IN)         :: elim
real(dp), INTENT(IN)         :: alim
 
COMPLEX (dp)  :: arg, asum, bsum, cwrk(16), cz, phi, sum, zb, zeta1, zeta2, &
                 zn, zr
real(dp)     :: aarg, aphi, ascle, ax, ay, fnn, gnn, gnu, rcz, x, yy
INTEGER       :: iform, init, nn, nw
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp, 0.0_dp)
real(dp), PARAMETER     :: aic = 1.265512123484645396_dp
 
nuf = 0
nn = n
x = real(z, kind=dp)
zr = z
IF (x < 0.0_dp) zr = -z
zb = zr
yy = aimag(zr)
ax = abs(x) * 1.73205080756887_dp
ay = abs(yy)
iform = 1
IF (ay > ax) iform = 2
gnu = max(fnu, 1.0_dp)
IF (ikflg /= 1) THEN
  fnn = nn
  gnn = fnu + fnn - 1.0_dp
  gnu = max(gnn, fnn)
END IF
!-----------------------------------------------------------------------
!     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
!     REAL PARTS OF ZETA1, ZETA2 AND ZB.  NO ATTEMPT IS MADE TO GET
!     THE SIGN OF THE IMAGINARY PART CORRECT.
!-----------------------------------------------------------------------
IF (iform /= 2) THEN
  init = 0
  CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum, cwrk)
  cz = -zeta1 + zeta2
ELSE
  zn = -zr * cmplx(0.0_dp, 1.0_dp, kind=dp)
  IF (yy <= 0.0_dp) THEN
    zn = conjg(-zn)
  END IF
  CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
  cz = -zeta1 + zeta2
  aarg = abs(arg)
END IF
IF (kode == 2) cz = cz - zb
IF (ikflg == 2) cz = -cz
aphi = abs(phi)
rcz = real(cz, kind=dp)
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
IF (rcz <= elim) THEN
  IF (rcz >= alim) THEN
    rcz = rcz + log(aphi)
    IF (iform == 2) rcz = rcz - 0.25_dp * log(aarg) - aic
    IF (rcz > elim) GO TO 80
  ELSE
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
    IF (rcz >= -elim) THEN
      IF (rcz > -alim) GO TO 40
      rcz = rcz + log(aphi)
      IF (iform == 2) rcz = rcz - 0.25_dp * log(aarg) - aic
      IF (rcz > -elim) GO TO 30
    END IF
 
    10 y(1:nn) = czero
    nuf = nn
    RETURN
 
    30 ascle = 1.0e+3 * tiny(0.0_dp) / tol
    cz = cz + log(phi)
    IF (iform /= 1) THEN
      cz = cz - 0.25_dp * log(arg) - aic
    END IF
    ax = exp(rcz) / tol
    ay = aimag(cz)
    cz = ax * cmplx(cos(ay), sin(ay), kind=dp)
    CALL cuchk(cz, nw, ascle, tol)
    IF (nw == 1) GO TO 10
  END IF
 
  40 IF (ikflg == 2) RETURN
  IF (n == 1) RETURN
!-----------------------------------------------------------------------
!     SET UNDERFLOWS ON I SEQUENCE
!-----------------------------------------------------------------------
  50 gnu = fnu + (nn-1)
  IF (iform /= 2) THEN
    init = 0
    CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum, cwrk)
    cz = -zeta1 + zeta2
  ELSE
    CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
    cz = -zeta1 + zeta2
    aarg = abs(arg)
  END IF
  IF (kode == 2) cz = cz - zb
  aphi = abs(phi)
  rcz = real(cz, kind=dp)
  IF (rcz >= -elim) THEN
    IF (rcz > -alim) RETURN
    rcz = rcz + log(aphi)
    IF (iform == 2) rcz = rcz - 0.25_dp * log(aarg) - aic
    IF (rcz > -elim) GO TO 70
  END IF
 
  60 y(nn) = czero
  nn = nn - 1
  nuf = nuf + 1
  IF (nn == 0) RETURN
  GO TO 50
 
  70 ascle = 1.0e+3 * tiny(0.0_dp) / tol
  cz = cz + log(phi)
  IF (iform /= 1) THEN
    cz = cz - 0.25_dp * log(arg) - aic
  END IF
  ax = exp(rcz) / tol
  ay = aimag(cz)
  cz = ax * cmplx(cos(ay), sin(ay), kind=dp)
  CALL cuchk(cz, nw, ascle, tol)
  IF (nw == 1) GO TO 60
  RETURN
END IF
 
80 nuf = -1
RETURN
END SUBROUTINE cuoik
 
 
 
SUBROUTINE cwrsk(zr, fnu, kode, n, y, nz, cw, tol, elim, alim)
!***BEGIN PROLOGUE  CWRSK
!***REFER TO  CBESI,CBESK
 
!     CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z) >= 0.0 BY
!     NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
 
!***ROUTINES CALLED  CBKNU,CRATI,R1MACH
!***END PROLOGUE  CWRSK
 
COMPLEX (dp), INTENT(IN)   :: zr
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
COMPLEX (dp), INTENT(OUT)  :: cw(2)
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cinu, cscl, ct, c1, c2, rct, st
real(dp)     :: act, acw, ascle, s1, s2, yy
INTEGER       :: i, nw
 
!-----------------------------------------------------------------------
!     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
CALL cbknu(zr, fnu, kode, 2, cw, nw, tol, elim, alim)
IF (nw == 0) THEN
  CALL crati(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 = cmplx(1.0_dp, 0.0_dp, kind=dp)
  IF (kode /= 1) THEN
    yy = aimag(zr)
    s1 = cos(yy)
    s2 = sin(yy)
    cinu = cmplx(s1, s2, kind=dp)
  END IF
!-----------------------------------------------------------------------
!     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 = abs(cw(2))
  ascle = 1.0e+3 * tiny(0.0_dp) / tol
  cscl = cmplx(1.0_dp, 0.0_dp, kind=dp)
  IF (acw <= ascle) THEN
    cscl = cmplx(1.0_dp/tol, 0.0_dp, kind=dp)
  ELSE
    ascle = 1.0_dp / ascle
    IF (acw >= ascle) THEN
      cscl = cmplx(tol, 0.0_dp, kind=dp)
    END IF
  END IF
  c1 = cw(1) * cscl
  c2 = cw(2) * cscl
  st = y(1)
!-----------------------------------------------------------------------
!     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 = abs(ct)
  rct = cmplx(1.0_dp/act, 0.0_dp, kind=dp)
  ct = conjg(ct) * rct
  cinu = cinu * rct * ct
  y(1) = cinu * cscl
  IF (n == 1) RETURN
  DO  i = 2, n
    cinu = st * cinu
    st = y(i)
    y(i) = cinu * cscl
  END DO
  RETURN
END IF
nz = -1
IF (nw == -2) nz = -2
RETURN
END SUBROUTINE cwrsk
 
 
 
SUBROUTINE cmlri(z, fnu, kode, n, y, nz, tol)
!***BEGIN PROLOGUE  CMLRI
!***REFER TO  CBESI,CBESK
 
!     CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z) >= 0.0 BY THE
!     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
 
!***ROUTINES CALLED  GAMLN,R1MACH
!***END PROLOGUE  CMLRI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: tol
 
COMPLEX (dp)  :: ck, cnorm, pt, p1, p2, rz, sum
real(dp)     :: ack, ak, ap, at, az, bk, fkap, fkk, flam, fnf, rho,  &
                 rho2, scle, tfnf, tst, x
INTEGER       :: i, iaz, ifnu, inu, itime, k, kk, km, m
 
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp), &
                            ctwo = (2.0_dp,0.0_dp)
 
scle = 1.0e+3 * tiny(0.0_dp) / tol
nz = 0
az = abs(z)
x = real(z, kind=dp)
iaz = int(az)
ifnu = int(fnu)
inu = ifnu + n - 1
at = iaz + 1
ck = cmplx(at, 0.0_dp, kind=dp) / z
rz = ctwo / z
p1 = czero
p2 = cone
ack = (at + 1.0_dp) / az
rho = ack + sqrt(ack*ack - 1.0_dp)
rho2 = rho * rho
tst = (rho2+rho2) / ((rho2-1.0_dp)*(rho-1.0_dp))
tst = tst / tol
!-----------------------------------------------------------------------
!     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
!-----------------------------------------------------------------------
ak = at
DO  i = 1, 80
  pt = p2
  p2 = p1 - ck * p2
  p1 = pt
  ck = ck + rz
  ap = abs(p2)
  IF (ap > tst*ak*ak) GO TO 20
  ak = ak + 1.0_dp
END DO
GO TO 90
 
20 i = i + 1
k = 0
IF (inu >= iaz) THEN
!-----------------------------------------------------------------------
!     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
!-----------------------------------------------------------------------
  p1 = czero
  p2 = cone
  at = inu + 1
  ck = cmplx(at, 0.0_dp, kind=dp) / z
  ack = at / az
  tst = sqrt(ack/tol)
  itime = 1
  DO  k = 1, 80
    pt = p2
    p2 = p1 - ck * p2
    p1 = pt
    ck = ck + rz
    ap = abs(p2)
    IF (ap >= tst) THEN
      IF (itime == 2) GO TO 40
      ack = abs(ck)
      flam = ack + sqrt(ack*ack - 1.0_dp)
      fkap = ap / abs(p1)
      rho = min(flam,fkap)
      tst = tst * sqrt(rho/(rho*rho - 1.0_dp))
      itime = 2
    END IF
  END DO
  GO TO 90
END IF
!-----------------------------------------------------------------------
!     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
!-----------------------------------------------------------------------
40 k = k + 1
kk = max(i+iaz, k+inu)
fkk = kk
p1 = czero
!-----------------------------------------------------------------------
!     SCALE P2 AND SUM BY SCLE
!-----------------------------------------------------------------------
p2 = cmplx(scle, 0.0_dp, kind=dp)
fnf = fnu - ifnu
tfnf = fnf + fnf
bk = gamln(fkk+tfnf+1.0_dp) - gamln(fkk+1.0_dp) - gamln(tfnf+1.0_dp)
bk = exp(bk)
sum = czero
km = kk - inu
DO  i = 1, km
  pt = p2
  p2 = p1 + cmplx(fkk+fnf, 0.0_dp, kind=dp) * rz * p2
  p1 = pt
  ak = 1.0_dp - tfnf / (fkk+tfnf)
  ack = bk * ak
  sum = sum + cmplx(ack+bk, 0.0_dp, kind=dp) * p1
  bk = ack
  fkk = fkk - 1.0_dp
END DO
y(n) = p2
IF (n /= 1) THEN
  DO  i = 2, n
    pt = p2
    p2 = p1 + cmplx(fkk+fnf, 0.0_dp, kind=dp) * rz * p2
    p1 = pt
    ak = 1.0_dp - tfnf / (fkk+tfnf)
    ack = bk * ak
    sum = sum + cmplx(ack+bk, 0.0_dp, kind=dp) * p1
    bk = ack
    fkk = fkk - 1.0_dp
    m = n - i + 1
    y(m) = p2
  END DO
END IF
IF (ifnu > 0) THEN
  DO  i = 1, ifnu
    pt = p2
    p2 = p1 + cmplx(fkk+fnf, 0.0_dp, kind=dp) * rz * p2
    p1 = pt
    ak = 1.0_dp - tfnf / (fkk+tfnf)
    ack = bk * ak
    sum = sum + cmplx(ack+bk, 0.0_dp, kind=dp) * p1
    bk = ack
    fkk = fkk - 1.0_dp
  END DO
END IF
pt = z
IF (kode == 2) pt = pt - x
p1 = -fnf * log(rz) + pt
ap = gamln(1.0_dp+fnf)
pt = p1 - ap
!-----------------------------------------------------------------------
!     THE DIVISION EXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
!     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
!-----------------------------------------------------------------------
p2 = p2 + sum
ap = abs(p2)
p1 = cmplx(1.0_dp/ap, 0.0_dp, kind=dp)
ck = exp(pt) * p1
pt = conjg(p2) * p1
cnorm = ck * pt
y(1:n) = y(1:n) * cnorm
RETURN
 
90 nz = -2
RETURN
END SUBROUTINE cmlri
 
 
 
SUBROUTINE cunhj(z, fnu, ipmtr, tol, phi, arg, zeta1, zeta2, asum, bsum)
!***BEGIN PROLOGUE  CUNHJ
!***REFER TO  CBESI,CBESK
 
!  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
!      CUNHJ 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*LOG((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  (NONE)
!***END PROLOGUE  CUNHJ
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: ipmtr
real(dp), INTENT(IN)      :: tol
COMPLEX (dp), INTENT(OUT)  :: phi
COMPLEX (dp), INTENT(OUT)  :: arg
COMPLEX (dp), INTENT(OUT)  :: zeta1
COMPLEX (dp), INTENT(OUT)  :: zeta2
COMPLEX (dp), INTENT(OUT)  :: asum
COMPLEX (dp), INTENT(OUT)  :: bsum
 
COMPLEX (dp)  :: cfnu, cr(14), dr(14), p(30), przth, ptfn, rtzta, rzth,  &
                 suma, sumb, tfn, t2, up(14), w, w2, za, zb, zc, zeta, zth
real(dp)     :: ang, ap(30), atol, aw2, azth, btol, fn13, fn23, pp, rfn13, &
                 rfnu, rfnu2, wi, wr, zci, zcr, zetai, zetar, zthi, zthr,  &
                 asumr, asumi, bsumr, bsumi, test, tstr, tsti, ac
INTEGER       :: ias, ibs, is, j, jr, ju, k, kmax, kp1, ks, l, lr, lrp1,  &
                 l1, l2, m
real(dp), PARAMETER  :: ar(14) = (/  &
    1.00000000000000000_dp, 1.04166666666666667d-01,  &
    8.35503472222222222d-02, 1.28226574556327160d-01,  &
    2.91849026464140464d-01, 8.81627267443757652d-01,  &
    3.32140828186276754_dp, 1.49957629868625547d+01,  &
    7.89230130115865181d+01, 4.74451538868264323d+02,  &
    3.20749009089066193d+03, 2.40865496408740049d+04,  &
    1.98923119169509794d+05, 1.79190200777534383d+06 /)
real(dp), PARAMETER  :: br(14) = (/  &
    1.00000000000000000_dp, -1.45833333333333333d-01,  &
   -9.87413194444444444d-02, -1.43312053915895062d-01,  &
   -3.17227202678413548d-01, -9.42429147957120249d-01,  &
   -3.51120304082635426_dp, -1.57272636203680451d+01,  &
   -8.22814390971859444d+01, -4.92355370523670524d+02,  &
   -3.31621856854797251d+03, -2.48276742452085896d+04,  &
   -2.04526587315129788d+05, -1.83844491706820990d+06 /)
real(dp), PARAMETER  :: c(105) = (/  &
    1.00000000000000000_dp, -2.08333333333333333d-01, 1.25000000000000000d-01, &
    3.34201388888888889d-01, -4.01041666666666667d-01, 7.03125000000000000d-02, &
   -1.02581259645061728_dp, 1.84646267361111111_dp, -8.91210937500000000d-01,  &
    7.32421875000000000d-02, 4.66958442342624743_dp, -1.12070026162229938d+01,  &
    8.78912353515625000_dp, -2.36408691406250000_dp, 1.12152099609375000d-01,  &
   -2.82120725582002449d+01, 8.46362176746007346d+01, -9.18182415432400174d+01,  &
    4.25349987453884549d+01, -7.36879435947963170_dp, 2.27108001708984375d-01,  &
    2.12570130039217123d+02, -7.65252468141181642d+02, 1.05999045252799988d+03,  &
   -6.99579627376132541d+02, 2.18190511744211590d+02, -2.64914304869515555d+01, &
    5.72501420974731445d-01, -1.91945766231840700d+03, 8.06172218173730938d+03, &
   -1.35865500064341374d+04, 1.16553933368645332d+04, -5.30564697861340311d+03, &
    1.20090291321635246d+03, -1.08090919788394656d+02, 1.72772750258445740_dp,  &
    2.02042913309661486d+04, -9.69805983886375135d+04, 1.92547001232531532d+05, &
   -2.03400177280415534d+05, 1.22200464983017460d+05, -4.11926549688975513d+04,  &
    7.10951430248936372d+03, -4.93915304773088012d+02, 6.07404200127348304_dp,  &
   -2.42919187900551333d+05, 1.31176361466297720d+06, -2.99801591853810675d+06, &
    3.76327129765640400d+06, -2.81356322658653411d+06, 1.26836527332162478d+06, &
   -3.31645172484563578d+05, 4.52187689813627263d+04, -2.49983048181120962d+03,  &
    2.43805296995560639d+01, 3.28446985307203782d+06, -1.97068191184322269d+07, &
    5.09526024926646422d+07, -7.41051482115326577d+07, 6.63445122747290267d+07, &
   -3.75671766607633513d+07, 1.32887671664218183d+07, -2.78561812808645469d+06, &
    3.08186404612662398d+05, -1.38860897537170405d+04, 1.10017140269246738d+02, &
   -4.93292536645099620d+07, 3.25573074185765749d+08, -9.39462359681578403d+08, &
    1.55359689957058006d+09, -1.62108055210833708d+09, 1.10684281682301447d+09, &
   -4.95889784275030309d+08, 1.42062907797533095d+08, -2.44740627257387285d+07, &
    2.24376817792244943d+06, -8.40054336030240853d+04, 5.51335896122020586d+02,  &
    8.14789096118312115d+08, -5.86648149205184723d+09, 1.86882075092958249d+10, &
   -3.46320433881587779d+10, 4.12801855797539740d+10, -3.30265997498007231d+10,  &
    1.79542137311556001d+10, -6.56329379261928433d+09, 1.55927986487925751d+09, &
   -2.25105661889415278d+08, 1.73951075539781645d+07, -5.49842327572288687d+05,  &
    3.03809051092238427d+03, -1.46792612476956167d+10, 1.14498237732025810d+11, &
   -3.99096175224466498d+11, 8.19218669548577329d+11, -1.09837515608122331d+12, &
    1.00815810686538209d+12, -6.45364869245376503d+11, 2.87900649906150589d+11, &
   -8.78670721780232657d+10, 1.76347306068349694d+10, -2.16716498322379509d+09,  &
    1.43157876718888981d+08, -3.87183344257261262d+06, 1.82577554742931747d+04 /)
real(dp), PARAMETER  :: alfa1(30) = (/  &
 -4.44444444444444444d-03, -9.22077922077922078d-04, -8.84892884892884893d-05, &
  1.65927687832449737d-04, 2.46691372741792910d-04, 2.65995589346254780d-04,  &
  2.61824297061500945d-04, 2.48730437344655609d-04, 2.32721040083232098d-04,  &
  2.16362485712365082d-04, 2.00738858762752355d-04, 1.86267636637545172d-04,  &
  1.73060775917876493d-04, 1.61091705929015752d-04, 1.50274774160908134d-04,  &
  1.40503497391269794d-04, 1.31668816545922806d-04, 1.23667445598253261d-04,  &
  1.16405271474737902d-04, 1.09798298372713369d-04, 1.03772410422992823d-04,  &
  9.82626078369363448d-05, 9.32120517249503256d-05, 8.85710852478711718d-05,  &
  8.42963105715700223d-05, 8.03497548407791151d-05, 7.66981345359207388d-05,  &
  7.33122157481777809d-05, 7.01662625163141333d-05, 6.72375633790160292d-05 /)
real(dp), PARAMETER  :: alfa2(30) = (/  &
  6.93735541354588974d-04, 2.32241745182921654d-04, -1.41986273556691197d-05, &
 -1.16444931672048640d-04, -1.50803558053048762d-04, -1.55121924918096223d-04, &
 -1.46809756646465549d-04, -1.33815503867491367d-04, -1.19744975684254051d-04, &
 -1.06184319207974020d-04, -9.37699549891194492d-05, -8.26923045588193274d-05, &
 -7.29374348155221211d-05, -6.44042357721016283d-05, -5.69611566009369048d-05, &
 -5.04731044303561628d-05, -4.48134868008882786d-05, -3.98688727717598864d-05, &
 -3.55400532972042498d-05, -3.17414256609022480d-05, -2.83996793904174811d-05, &
 -2.54522720634870566d-05, -2.28459297164724555d-05, -2.05352753106480604d-05, &
 -1.84816217627666085d-05, -1.66519330021393806d-05, -1.50179412980119482d-05, &
 -1.35554031379040526d-05, -1.22434746473858131d-05, -1.10641884811308169d-05 /)
real(dp), PARAMETER  :: alfa3(30) = (/  &
 -3.54211971457743841d-04, -1.56161263945159416d-04, 3.04465503594936410d-05, &
  1.30198655773242693d-04, 1.67471106699712269d-04, 1.70222587683592569d-04,  &
  1.56501427608594704d-04, 1.36339170977445120d-04, 1.14886692029825128d-04,  &
  9.45869093034688111d-05, 7.64498419250898258d-05, 6.07570334965197354d-05,  &
  4.74394299290508799d-05, 3.62757512005344297d-05, 2.69939714979224901d-05,  &
  1.93210938247939253d-05, 1.30056674793963203d-05, 7.82620866744496661d-06,  &
  3.59257485819351583d-06, 1.44040049814251817d-07, -2.65396769697939116d-06, &
 -4.91346867098485910d-06, -6.72739296091248287d-06, -8.17269379678657923d-06, &
 -9.31304715093561232d-06, -1.02011418798016441d-05, -1.08805962510592880d-05, &
 -1.13875481509603555d-05, -1.17519675674556414d-05, -1.19987364870944141d-05 /)
real(dp), PARAMETER  :: alfa4(30) = (/  &
  3.78194199201772914d-04, 2.02471952761816167d-04, -6.37938506318862408d-05, &
 -2.38598230603005903d-04, -3.10916256027361568d-04, -3.13680115247576316d-04, &
 -2.78950273791323387d-04, -2.28564082619141374d-04, -1.75245280340846749d-04, &
 -1.25544063060690348d-04, -8.22982872820208365d-05, -4.62860730588116458d-05, &
 -1.72334302366962267d-05, 5.60690482304602267d-06, 2.31395443148286800d-05,  &
  3.62642745856793957d-05, 4.58006124490188752d-05, 5.24595294959114050d-05,  &
  5.68396208545815266d-05, 5.94349820393104052d-05, 6.06478527578421742d-05,  &
  6.08023907788436497d-05, 6.01577894539460388d-05, 5.89199657344698500d-05,  &
  5.72515823777593053d-05, 5.52804375585852577d-05, 5.31063773802880170d-05,  &
  5.08069302012325706d-05, 4.84418647620094842d-05, 4.60568581607475370d-05 /)
real(dp), PARAMETER  :: alfa5(30) = (/  &
 -6.91141397288294174d-04, -4.29976633058871912d-04, 1.83067735980039018d-04, &
  6.60088147542014144d-04, 8.75964969951185931d-04, 8.77335235958235514d-04,  &
  7.49369585378990637d-04, 5.63832329756980918d-04, 3.68059319971443156d-04,  &
  1.88464535514455599d-04, 3.70663057664904149d-05, -8.28520220232137023d-05, &
 -1.72751952869172998d-04, -2.36314873605872983d-04, -2.77966150694906658d-04, &
 -3.02079514155456919d-04, -3.12594712643820127d-04, -3.12872558758067163d-04, &
 -3.05678038466324377d-04, -2.93226470614557331d-04, -2.77255655582934777d-04, &
 -2.59103928467031709d-04, -2.39784014396480342d-04, -2.20048260045422848d-04, &
 -2.00443911094971498d-04, -1.81358692210970687d-04, -1.63057674478657464d-04, &
 -1.45712672175205844d-04, -1.29425421983924587d-04, -1.14245691942445952d-04 /)
real(dp), PARAMETER  :: alfa6(30) = (/  &
  1.92821964248775885d-03, 1.35592576302022234d-03, -7.17858090421302995d-04, &
 -2.58084802575270346d-03, -3.49271130826168475d-03, -3.46986299340960628d-03, &
 -2.82285233351310182d-03, -1.88103076404891354d-03, -8.89531718383947600d-04, &
  3.87912102631035228d-06, 7.28688540119691412d-04, 1.26566373053457758d-03,  &
  1.62518158372674427d-03, 1.83203153216373172d-03, 1.91588388990527909d-03,  &
  1.90588846755546138d-03, 1.82798982421825727d-03, 1.70389506421121530d-03,  &
  1.55097127171097686d-03, 1.38261421852276159d-03, 1.20881424230064774d-03,  &
  1.03676532638344962d-03, 8.71437918068619115d-04, 7.16080155297701002d-04,  &
  5.72637002558129372d-04, 4.42089819465802277d-04, 3.24724948503090564d-04,  &
  2.20342042730246599d-04, 1.28412898401353882d-04, 4.82005924552095464d-05 /)
real(dp)  :: alfa(180)
real(dp), PARAMETER  :: beta1(30) = (/  &
  1.79988721413553309d-02, 5.59964911064388073d-03, 2.88501402231132779d-03,  &
  1.80096606761053941d-03, 1.24753110589199202d-03, 9.22878876572938311d-04,  &
  7.14430421727287357d-04, 5.71787281789704872d-04, 4.69431007606481533d-04,  &
  3.93232835462916638d-04, 3.34818889318297664d-04, 2.88952148495751517d-04,  &
  2.52211615549573284d-04, 2.22280580798883327d-04, 1.97541838033062524d-04,  &
  1.76836855019718004d-04, 1.59316899661821081d-04, 1.44347930197333986d-04,  &
  1.31448068119965379d-04, 1.20245444949302884d-04, 1.10449144504599392d-04,  &
  1.01828770740567258d-04, 9.41998224204237509d-05, 8.74130545753834437d-05,  &
  8.13466262162801467d-05, 7.59002269646219339d-05, 7.09906300634153481d-05,  &
  6.65482874842468183d-05, 6.25146958969275078d-05, 5.88403394426251749d-05 /)
real(dp), PARAMETER  :: beta2(30) = (/  &
 -1.49282953213429172d-03, -8.78204709546389328d-04, -5.02916549572034614d-04, &
 -2.94822138512746025d-04, -1.75463996970782828d-04, -1.04008550460816434d-04, &
 -5.96141953046457895d-05, -3.12038929076098340d-05, -1.26089735980230047d-05, &
 -2.42892608575730389d-07, 8.05996165414273571d-06, 1.36507009262147391d-05,  &
  1.73964125472926261d-05, 1.98672978842133780d-05, 2.14463263790822639d-05,  &
  2.23954659232456514d-05, 2.28967783814712629d-05, 2.30785389811177817d-05,  &
  2.30321976080909144d-05, 2.28236073720348722d-05, 2.25005881105292418d-05,  &
  2.20981015361991429d-05, 2.16418427448103905d-05, 2.11507649256220843d-05,  &
  2.06388749782170737d-05, 2.01165241997081666d-05, 1.95913450141179244d-05,  &
  1.90689367910436740d-05, 1.85533719641636667d-05, 1.80475722259674218d-05 /)
real(dp), PARAMETER  :: beta3(30) = (/  &
  5.52213076721292790d-04, 4.47932581552384646d-04, 2.79520653992020589d-04,  &
  1.52468156198446602d-04, 6.93271105657043598d-05, 1.76258683069991397d-05,  &
 -1.35744996343269136d-05, -3.17972413350427135d-05, -4.18861861696693365d-05, &
 -4.69004889379141029d-05, -4.87665447413787352d-05, -4.87010031186735069d-05, &
 -4.74755620890086638d-05, -4.55813058138628452d-05, -4.33309644511266036d-05, &
 -4.09230193157750364d-05, -3.84822638603221274d-05, -3.60857167535410501d-05, &
 -3.37793306123367417d-05, -3.15888560772109621d-05, -2.95269561750807315d-05, &
 -2.75978914828335759d-05, -2.58006174666883713d-05, -2.41308356761280200d-05, &
 -2.25823509518346033d-05, -2.11479656768912971d-05, -1.98200638885294927d-05, &
 -1.85909870801065077d-05, -1.74532699844210224d-05, -1.63997823854497997d-05 /)
real(dp), PARAMETER  :: beta4(30) = (/  &
 -4.74617796559959808d-04, -4.77864567147321487d-04, -3.20390228067037603d-04, &
 -1.61105016119962282d-04, -4.25778101285435204d-05, 3.44571294294967503d-05, &
  7.97092684075674924d-05, 1.03138236708272200d-04, 1.12466775262204158d-04,  &
  1.13103642108481389d-04, 1.08651634848774268d-04, 1.01437951597661973d-04,  &
  9.29298396593363896d-05, 8.40293133016089978d-05, 7.52727991349134062d-05,  &
  6.69632521975730872d-05, 5.92564547323194704d-05, 5.22169308826975567d-05,  &
  4.58539485165360646d-05, 4.01445513891486808d-05, 3.50481730031328081d-05,  &
  3.05157995034346659d-05, 2.64956119950516039d-05, 2.29363633690998152d-05,  &
  1.97893056664021636d-05, 1.70091984636412623d-05, 1.45547428261524004d-05,  &
  1.23886640995878413d-05, 1.04775876076583236d-05, 8.79179954978479373d-06 /)
real(dp), PARAMETER  :: beta5(30) = (/  &
  7.36465810572578444d-04, 8.72790805146193976d-04, 6.22614862573135066d-04,  &
  2.85998154194304147d-04, 3.84737672879366102d-06, -1.87906003636971558d-04, &
 -2.97603646594554535d-04, -3.45998126832656348d-04, -3.53382470916037712d-04, &
 -3.35715635775048757d-04, -3.04321124789039809d-04, -2.66722723047612821d-04, &
 -2.27654214122819527d-04, -1.89922611854562356d-04, -1.55058918599093870d-04, &
 -1.23778240761873630d-04, -9.62926147717644187d-05, -7.25178327714425337d-05, &
 -5.22070028895633801d-05, -3.50347750511900522d-05, -2.06489761035551757d-05, &
 -8.70106096849767054d-06, 1.13698686675100290d-06, 9.16426474122778849d-06,  &
  1.56477785428872620d-05, 2.08223629482466847d-05, 2.48923381004595156d-05,  &
  2.80340509574146325d-05, 3.03987774629861915d-05, 3.21156731406700616d-05 /)
real(dp), PARAMETER  :: beta6(30) = (/  &
 -1.80182191963885708d-03, -2.43402962938042533d-03, -1.83422663549856802d-03, &
 -7.62204596354009765d-04, 2.39079475256927218d-04, 9.49266117176881141d-04,  &
  1.34467449701540359d-03, 1.48457495259449178d-03, 1.44732339830617591d-03,  &
  1.30268261285657186d-03, 1.10351597375642682d-03, 8.86047440419791759d-04,  &
  6.73073208165665473d-04, 4.77603872856582378d-04, 3.05991926358789362d-04,  &
  1.60315694594721630d-04, 4.00749555270613286d-05, -5.66607461635251611d-05, &
 -1.32506186772982638d-04, -1.90296187989614057d-04, -2.32811450376937408d-04, &
 -2.62628811464668841d-04, -2.82050469867598672d-04, -2.93081563192861167d-04, &
 -2.97435962176316616d-04, -2.96557334239348078d-04, -2.91647363312090861d-04, &
 -2.83696203837734166d-04, -2.73512317095673346d-04, -2.61750155806768580d-04 /)
real(dp), PARAMETER  :: beta7(30) = (/  &
  6.38585891212050914d-03, 9.62374215806377941d-03, 7.61878061207001043d-03,  &
  2.83219055545628054d-03, -2.09841352012720090d-03, -5.73826764216626498d-03, &
 -7.70804244495414620d-03, -8.21011692264844401d-03, -7.65824520346905413d-03, &
 -6.47209729391045177d-03, -4.99132412004966473d-03, -3.45612289713133280d-03, &
 -2.01785580014170775d-03, -7.59430686781961401d-04, 2.84173631523859138d-04, &
  1.10891667586337403d-03, 1.72901493872728771d-03, 2.16812590802684701d-03,  &
  2.45357710494539735d-03, 2.61281821058334862d-03, 2.67141039656276912d-03,  &
  2.65203073395980430d-03, 2.57411652877287315d-03, 2.45389126236094427d-03,  &
  2.30460058071795494d-03, 2.13684837686712662d-03, 1.95896528478870911d-03,  &
  1.77737008679454412d-03, 1.59690280765839059d-03, 1.42111975664438546d-03 /)
real(dp)  :: beta(210)
real(dp), PARAMETER  :: gama(30) = (/  &
  6.29960524947436582d-01, 2.51984209978974633d-01, 1.54790300415655846d-01, &
  1.10713062416159013d-01, 8.57309395527394825d-02, 6.97161316958684292d-02, &
  5.86085671893713576d-02, 5.04698873536310685d-02, 4.42600580689154809d-02, &
  3.93720661543509966d-02, 3.54283195924455368d-02, 3.21818857502098231d-02, &
  2.94646240791157679d-02, 2.71581677112934479d-02, 2.51768272973861779d-02, &
  2.34570755306078891d-02, 2.19508390134907203d-02, 2.06210828235646240d-02, &
  1.94388240897880846d-02, 1.83810633800683158d-02, 1.74293213231963172d-02, &
  1.65685837786612353d-02, 1.57865285987918445d-02, 1.50729501494095594d-02, &
  1.44193250839954639d-02, 1.38184805735341786d-02, 1.32643378994276568d-02, &
  1.27517121970498651d-02, 1.22761545318762767d-02, 1.18338262398482403d-02 /)
real(dp), PARAMETER  :: ex1 = 3.33333333333333333d-01,  &
    ex2 = 6.66666666666666667d-01, hpi = 1.57079632679489662_dp,  &
    pi = 3.14159265358979324_dp, thpi = 4.71238898038468986_dp
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
 
! Associate arrays alfa & beta
 
alfa(  1: 30) = alfa1
alfa( 31: 60) = alfa2
alfa( 61: 90) = alfa3
alfa( 91:120) = alfa4
alfa(121:150) = alfa5
alfa(151:180) = alfa6
beta(  1: 30) = beta1
beta( 31: 60) = beta2
beta( 61: 90) = beta3
beta( 91:120) = beta4
beta(121:150) = beta5
beta(151:180) = beta6
beta(181:210) = beta7
 
rfnu = 1.0_dp / fnu
!     ZB = Z*CMPLX(RFNU,0.0_dp)
!-----------------------------------------------------------------------
!     OVERFLOW TEST (Z/FNU TOO SMALL)
!-----------------------------------------------------------------------
tstr = real(z, kind=dp)
tsti = aimag(z)
test = tiny(0.0_dp) * 1.0e+3
ac = fnu * test
IF (abs(tstr) <= ac .AND. abs(tsti) <= ac) THEN
  ac = 2.0_dp * abs(log(test)) + fnu
  zeta1 = ac
  zeta2 = fnu
  phi = cone
  arg = cone
  RETURN
END IF
zb = z * rfnu
rfnu2 = rfnu * rfnu
!-----------------------------------------------------------------------
!     COMPUTE IN THE FOURTH QUADRANT
!-----------------------------------------------------------------------
fn13 = fnu ** ex1
fn23 = fn13 * fn13
rfn13 = 1.0_dp/fn13
w2 = cone - zb * zb
aw2 = abs(w2)
IF (aw2 > 0.25_dp) GO TO 110
!-----------------------------------------------------------------------
!     POWER SERIES FOR ABS(W2) <= 0.25_dp
!-----------------------------------------------------------------------
k = 1
 
p(1) = cone
suma = gama(1)
ap(1) = 1.0_dp
IF (aw2 >= tol) THEN
  DO  k = 2, 30
    p(k) = p(k-1) * w2
    suma = suma + p(k) * gama(k)
    ap(k) = ap(k-1) * aw2
    IF (ap(k) < tol) GO TO 20
  END DO
  k = 30
END IF
 
20 kmax = k
zeta = w2 * suma
arg = zeta * fn23
za = sqrt(suma)
zeta2 = sqrt(w2) * fnu
zeta1 = zeta2 * (cone + zeta*za*ex2)
za = za + za
phi = sqrt(za) * rfn13
IF (ipmtr /= 1) THEN
!-----------------------------------------------------------------------
!     SUM SERIES FOR ASUM AND BSUM
!-----------------------------------------------------------------------
  sumb = czero
  DO k = 1, kmax
    sumb = sumb + p(k)*beta(k)
  END DO
  asum = czero
  bsum = sumb
  l1 = 0
  l2 = 30
  btol = tol * (abs(real(bsum)) + abs(aimag(bsum)))
  atol = tol
  pp = 1.0_dp
  ias = 0
  ibs = 0
  IF (rfnu2 >= tol) THEN
    DO  is = 2, 7
      atol = atol / rfnu2
      pp = pp * rfnu2
      IF (ias /= 1) THEN
        suma = czero
        DO  k = 1, kmax
          m = l1 + k
          suma = suma + p(k) * alfa(m)
          IF (ap(k) < atol) EXIT
        END DO
        asum = asum + suma * pp
        IF (pp < tol) ias = 1
      END IF
      IF (ibs /= 1) THEN
        sumb = czero
        DO  k = 1, kmax
          m = l2 + k
          sumb = sumb + p(k) * beta(m)
          IF (ap(k) < atol) EXIT
        END DO
        bsum = bsum + sumb * pp
        IF (pp < btol) ibs = 1
      END IF
      IF (ias == 1 .AND. ibs == 1) EXIT
      l1 = l1 + 30
      l2 = l2 + 30
    END DO
  END IF
 
  asum = asum + cone
  pp = rfnu * rfn13
  bsum = bsum * pp
END IF
 
100 RETURN
!-----------------------------------------------------------------------
!     ABS(W2) > 0.25_dp
!-----------------------------------------------------------------------
110 w = sqrt(w2)
wr = real(w, kind=dp)
wi = aimag(w)
IF (wr < 0.0_dp) wr = 0.0_dp
IF (wi < 0.0_dp) wi = 0.0_dp
w = cmplx(wr, wi, kind=dp)
za = (cone+w) / zb
zc = log(za)
zcr = real(zc, kind=dp)
zci = aimag(zc)
IF (zci < 0.0_dp) zci = 0.0_dp
IF (zci > hpi) zci = hpi
IF (zcr < 0.0_dp) zcr = 0.0_dp
zc = cmplx(zcr, zci, kind=dp)
zth = (zc-w) * 1.5_dp
cfnu = cmplx(fnu, 0.0_dp, kind=dp)
zeta1 = zc * cfnu
zeta2 = w * cfnu
azth = abs(zth)
zthr = real(zth, kind=dp)
zthi = aimag(zth)
ang = thpi
IF (zthr < 0.0_dp .OR. zthi >= 0.0_dp) THEN
  ang = hpi
  IF (zthr /= 0.0_dp) THEN
    ang = atan(zthi/zthr)
    IF (zthr < 0.0_dp) ang = ang + pi
  END IF
END IF
pp = azth ** ex2
ang = ang * ex2
zetar = pp * cos(ang)
zetai = pp * sin(ang)
IF (zetai < 0.0_dp) zetai = 0.0_dp
zeta = cmplx(zetar, zetai, kind=dp)
arg = zeta * fn23
rtzta = zth / zeta
za = rtzta / w
phi = sqrt(za+za) * rfn13
IF (ipmtr == 1) GO TO 100
tfn = cmplx(rfnu, 0.0_dp, kind=dp) / w
rzth = cmplx(rfnu, 0.0_dp, kind=dp) / zth
zc = rzth * ar(2)
t2 = cone / w2
up(2) = (t2*c(2) + c(3)) * tfn
bsum = up(2) + zc
asum = czero
IF (rfnu >= tol) THEN
  przth = rzth
  ptfn = tfn
  up(1) = cone
  pp = 1.0_dp
  bsumr = real(bsum, kind=dp)
  bsumi = aimag(bsum)
  btol = tol * (abs(bsumr) + abs(bsumi))
  ks = 0
  kp1 = 2
  l = 3
  ias = 0
  ibs = 0
  DO  lr = 2, 12, 2
    lrp1 = lr + 1
!-----------------------------------------------------------------------
!     COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN
!     NEXT SUMA AND SUMB
!-----------------------------------------------------------------------
    DO  k = lr, lrp1
      ks = ks + 1
      kp1 = kp1 + 1
      l = l + 1
      za = cmplx(c(l), 0.0_dp, kind=dp)
      DO  j = 2, kp1
        l = l + 1
        za = za * t2 + c(l)
      END DO
      ptfn = ptfn * tfn
      up(kp1) = ptfn * za
      cr(ks) = przth * br(ks+1)
      przth = przth * rzth
      dr(ks) = przth * ar(ks+2)
    END DO
    pp = pp * rfnu2
    IF (ias /= 1) THEN
      suma = up(lrp1)
      ju = lrp1
      DO  jr = 1, lr
        ju = ju - 1
        suma = suma + cr(jr) * up(ju)
      END DO
      asum = asum + suma
      asumr = real(asum, kind=dp)
      asumi = aimag(asum)
      test = abs(asumr) + abs(asumi)
      IF (pp < tol .AND. test < tol) ias = 1
    END IF
    IF (ibs /= 1) THEN
      sumb = up(lr+2) + up(lrp1) * zc
      ju = lrp1
      DO  jr = 1, lr
        ju = ju - 1
        sumb = sumb + dr(jr) * up(ju)
      END DO
      bsum = bsum + sumb
      bsumr = real(bsum, kind=dp)
      bsumi = aimag(bsum)
      test = abs(bsumr) + abs(bsumi)
      IF (pp < btol .AND. test < tol) ibs = 1
    END IF
    IF (ias == 1 .AND. ibs == 1) GO TO 170
  END DO
END IF
 
170 asum = asum + cone
bsum = -bsum * rfn13 / rtzta
GO TO 100
END SUBROUTINE cunhj
 
 
 
SUBROUTINE cseri(z, fnu, kode, n, y, nz, tol, elim, alim)
!***BEGIN PROLOGUE  CSERI
!***REFER TO  CBESI,CBESK
 
!     CSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0.0 BY
!     MEANS OF THE POWER SERIES FOR LARGE ABS(Z) IN THE
!     REGION ABS(Z) <= 2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
!     NZ > 0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
!     DUE TO UNDERFLOW. NZ < 0 MEANS UNDERFLOW OCCURRED, BUT THE
!     CONDITION ABS(Z) <= 2*SQRT(FNU+1) WAS VIOLATED AND THE
!     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
 
!***ROUTINES CALLED  CUCHK,GAMLN,R1MACH
!***END PROLOGUE  CSERI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: ak1, ck, coef, crsc, cz, hz, rz, s1, s2, w(2)
real(dp)     :: aa, acz, ak, arm, ascle, atol, az, dfnu,  &
                 fnup, rak1, rs, rtr1, s, ss, x
INTEGER       :: i, ib, iflag, il, k, l, m, nn, nw
real(dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
 
nz = 0
az = abs(z)
IF (az /= 0.0_dp) THEN
  x = real(z, kind=dp)
  arm = 1.0d+3 * tiny(0.0_dp)
  rtr1 = sqrt(arm)
  crsc = cmplx(1.0_dp, 0.0_dp, kind=dp)
  iflag = 0
  IF (az >= arm) THEN
    hz = z * 0.5_dp
    cz = czero
    IF (az > rtr1) cz = hz * hz
    acz = abs(cz)
    nn = n
    ck = log(hz)
 
    10 dfnu = fnu + (nn-1)
    fnup = dfnu + 1.0_dp
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
    ak1 = ck * dfnu
    ak = gamln(fnup)
    ak1 = ak1 - ak
    IF (kode == 2) ak1 = ak1 - x
    rak1 = real(ak1, kind=dp)
    IF (rak1 > -elim) GO TO 30
 
    20 nz = nz + 1
    y(nn) = czero
    IF (acz > dfnu) GO TO 120
    nn = nn - 1
    IF (nn == 0) RETURN
    GO TO 10
 
    30 IF (rak1 <= -alim) THEN
      iflag = 1
      ss = 1.0_dp / tol
      crsc = cmplx(tol, 0.0_dp, kind=dp)
      ascle = arm * ss
    END IF
    ak = aimag(ak1)
    aa = exp(rak1)
    IF (iflag == 1) aa = aa * ss
    coef = aa * cmplx(cos(ak), sin(ak), kind=dp)
    atol = tol * acz / fnup
    il = min(2,nn)
    DO  i = 1, il
      dfnu = fnu + (nn-i)
      fnup = dfnu + 1.0_dp
      s1 = cone
      IF (acz >= tol*fnup) THEN
        ak1 = cone
        ak = fnup + 2.0_dp
        s = fnup
        aa = 2.0_dp
 
        40 rs = 1.0_dp / s
        ak1 = ak1 * cz * rs
        s1 = s1 + ak1
        s = s + ak
        ak = ak + 2.0_dp
        aa = aa * acz * rs
        IF (aa > atol) GO TO 40
      END IF
      m = nn - i + 1
      s2 = s1 * coef
      w(i) = s2
      IF (iflag /= 0) THEN
        CALL cuchk(s2, nw, ascle, tol)
        IF (nw /= 0) GO TO 20
      END IF
      y(m) = s2 * crsc
      IF (i /= il) coef = coef * dfnu / hz
    END DO
    IF (nn <= 2) RETURN
    k = nn - 2
    ak = k
    rz = (cone+cone) / z
    IF (iflag == 1) GO TO 80
    ib = 3
 
    60 DO  i = ib, nn
      y(k) = cmplx(ak+fnu, 0.0_dp, kind=dp) * rz * y(k+1) + y(k+2)
      ak = ak - 1.0_dp
      k = k - 1
    END DO
    RETURN
!-----------------------------------------------------------------------
!     RECUR BACKWARD WITH SCALED VALUES
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
!     UNDERFLOW LIMIT = ASCLE = TINY(0.0_dp)*CSCL*1.0E+3
!-----------------------------------------------------------------------
    80 s1 = w(1)
    s2 = w(2)
    DO  l = 3, nn
      ck = s2
      s2 = s1 + cmplx(ak+fnu, 0.0_dp, kind=dp) * rz * s2
      s1 = ck
      ck = s2 * crsc
      y(k) = ck
      ak = ak - 1.0_dp
      k = k - 1
      IF (abs(ck) > ascle) GO TO 100
    END DO
    RETURN
 
    100 ib = l + 1
    IF (ib > nn) RETURN
    GO TO 60
  END IF
  nz = n
  IF (fnu == 0.0_dp) nz = nz - 1
END IF
y(1) = czero
IF (fnu == 0.0_dp) y(1) = cone
IF (n == 1) RETURN
y(2:n) = czero
RETURN
!-----------------------------------------------------------------------
!     RETURN WITH NZ < 0 IF ABS(Z*Z/4) > FNU+N-NZ-1 COMPLETE
!     THE CALCULATION IN CBINU WITH N=N-ABS(NZ)
!-----------------------------------------------------------------------
120 nz = -nz
RETURN
END SUBROUTINE cseri
 
 
 
SUBROUTINE casyi(z, fnu, kode, n, y, nz, rl, tol, elim, alim)
!***BEGIN PROLOGUE  CASYI
!***REFER TO  CBESI,CBESK
 
!     CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0.0 BY
!     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE
!     REGION ABS(Z) > MAX(RL,FNU*FNU/2).  NZ=0 IS A NORMAL RETURN.
!     NZ < 0 INDICATES AN OVERFLOW ON KODE=1.
 
!***ROUTINES CALLED  R1MACH
!***END PROLOGUE  CASYI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: rl
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: ak1, ck, cs1, cs2, cz, dk, ez, p1, rz, s2
real(dp)     :: aa, acz, aez, ak, arg, arm, atol, az, bb, bk, dfnu,  &
                 dnu2, fdn, rtr1, s, sgn, sqk, x, yy
INTEGER       :: i, ib, il, inu, j, jl, k, koded, m, nn
 
real(dp), PARAMETER     :: pi = 3.14159265358979324_dp, rtpi = 0.159154943091895336_dp
! rtpi = reciprocal of 2.pi
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
 
nz = 0
az = abs(z)
x = real(z, kind=dp)
arm = 1.0d+3 * tiny(0.0_dp)
rtr1 = sqrt(arm)
il = min(2,n)
dfnu = fnu + (n-il)
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
ak1 = rtpi / z
ak1 = sqrt(ak1)
cz = z
IF (kode == 2) cz = z - x
acz = real(cz, kind=dp)
IF (abs(acz) <= elim) THEN
  dnu2 = dfnu + dfnu
  koded = 1
  IF (.NOT.(abs(acz) > alim .AND. n > 2)) THEN
    koded = 0
    ak1 = ak1 * exp(cz)
  END IF
  fdn = 0.0_dp
  IF (dnu2 > rtr1) fdn = dnu2 * dnu2
  ez = z * 8.0_dp
!-----------------------------------------------------------------------
!     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.0_dp * az
  s = tol / aez
  jl = int(rl + rl + 2.0_dp)
  yy = aimag(z)
  p1 = czero
  IF (yy /= 0.0_dp) THEN
!-----------------------------------------------------------------------
!     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
!     SIGNIFICANCE WHEN FNU OR N IS LARGE
!-----------------------------------------------------------------------
    inu = int(fnu)
    arg = (fnu - inu) * pi
    inu = inu + n - il
    ak = -sin(arg)
    bk = cos(arg)
    IF (yy < 0.0_dp) bk = -bk
    p1 = cmplx(ak, bk, kind=dp)
    IF (mod(inu,2) == 1) p1 = -p1
  END IF
  DO  k = 1, il
    sqk = fdn - 1.0_dp
    atol = s * abs(sqk)
    sgn = 1.0_dp
    cs1 = cone
    cs2 = cone
    ck = cone
    ak = 0.0_dp
    aa = 1.0_dp
    bb = aez
    dk = ez
    DO  j = 1, jl
      ck = ck * sqk / dk
      cs2 = cs2 + ck
      sgn = -sgn
      cs1 = cs1 + ck * sgn
      dk = dk + ez
      aa = aa * abs(sqk) / bb
      bb = bb + aez
      ak = ak + 8.0_dp
      sqk = sqk - ak
      IF (aa <= atol) GO TO 20
    END DO
    GO TO 60
 
    20 s2 = cs1
    IF (x+x < elim) s2 = s2 + p1 * cs2 * exp(-z-z)
    fdn = fdn + 8.0_dp * dfnu + 4.0_dp
    p1 = -p1
    m = n - il + k
    y(m) = s2 * ak1
  END DO
  IF (n <= 2) RETURN
  nn = n
  k = nn - 2
  ak = k
  rz = (cone+cone) / z
  ib = 3
  DO  i = ib, nn
    y(k) = cmplx(ak+fnu, 0.0_dp, kind=dp) * rz * y(k+1) + y(k+2)
    ak = ak - 1.0_dp
    k = k - 1
  END DO
  IF (koded == 0) RETURN
  ck = exp(cz)
  y(1:nn) = y(1:nn) * ck
  RETURN
END IF
nz = -1
RETURN
 
60 nz = -2
RETURN
END SUBROUTINE casyi
 
 
 
SUBROUTINE cbunk(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
!***BEGIN PROLOGUE  CBUNK
!***REFER TO  CBESK,CBESH
 
!     CBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU > FNUL.
!     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR K(FNU,Z)
!     IN CUNK1 AND THE EXPANSION FOR H(2,FNU,Z) IN CUNK2
 
!***ROUTINES CALLED  CUNK1,CUNK2
!***END PROLOGUE  CBUNK
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: mr
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN OUT)  :: tol
real(dp), INTENT(IN OUT)  :: elim
real(dp), INTENT(IN OUT)  :: alim
 
real(dp)  :: ax, ay, xx, yy
 
nz = 0
xx = real(z, kind=dp)
yy = aimag(z)
ax = abs(xx) * 1.7321_dp
ay = abs(yy)
IF (ay <= ax) THEN
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR K(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!-----------------------------------------------------------------------
  CALL cunk1(z, fnu, kode, mr, n, y, nz, 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
!-----------------------------------------------------------------------
  CALL cunk2(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
END IF
RETURN
END SUBROUTINE cbunk
 
 
 
SUBROUTINE cunk1(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
!***BEGIN PROLOGUE  CUNK1
!***REFER TO  CBESK
 
!     CUNK1 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  CS1S2,CUCHK,CUNIK,R1MACH
!***END PROLOGUE  CUNK1
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: mr
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cfn, ck, crsc, cs, cscl, csgn, cspn, csr(3), css(3),  &
                 cwrk(16,3), cy(2), c1, c2, phi(2), rz, sum(2), s1, s2, &
                 zeta1(2), zeta2(2), zr, phid, zeta1d, zeta2d, sumd
real(dp)     :: ang, aphi, asc, ascle, bry(3), cpn, c2i, c2m, c2r,  &
                 fmr, fn, fnf, rs1, sgn, spn, x
INTEGER       :: i, ib, iflag, ifn, il, init(2), inu, iuf, k, kdflg, kflag, &
                 kk, m, nw, j, ipard, initd, ic
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
real(dp), PARAMETER     :: pi = 3.14159265358979324_dp
 
kdflg = 1
nz = 0
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
!     THE UNDERFLOW LIMIT
!-----------------------------------------------------------------------
cscl = 1.0_dp/tol
crsc = tol
css(1) = cscl
css(2) = cone
css(3) = crsc
csr(1) = crsc
csr(2) = cone
csr(3) = cscl
bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
bry(2) = 1.0_dp / bry(1)
bry(3) = huge(0.0_dp)
x = real(z, kind=dp)
zr = z
IF (x < 0.0_dp) zr = -z
j = 2
DO  i = 1, n
!-----------------------------------------------------------------------
!     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!-----------------------------------------------------------------------
  j = 3 - j
  fn = fnu + (i-1)
  init(j) = 0
  CALL cunik(zr, fn, 2, 0, tol, init(j), phi(j), zeta1(j), zeta2(j), sum(j), &
             cwrk(1:,j))
  IF (kode /= 1) THEN
    cfn = fn
    s1 = zeta1(j) - cfn * (cfn/(zr + zeta2(j)))
  ELSE
    s1 = zeta1(j) - zeta2(j)
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) <= elim) THEN
    IF (kdflg == 1) kflag = 2
    IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
      aphi = abs(phi(j))
      rs1 = rs1 + log(aphi)
      IF (abs(rs1) > elim) GO TO 10
      IF (kdflg == 1) kflag = 1
      IF (rs1 >= 0.0_dp) THEN
        IF (kdflg == 1) kflag = 3
      END IF
    END IF
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
    s2 = phi(j) * sum(j)
    c2r = real(s1, kind=dp)
    c2i = aimag(s1)
    c2m = exp(c2r) * real(css(kflag), kind=dp)
    s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
    s2 = s2 * s1
    IF (kflag == 1) THEN
      CALL cuchk(s2, nw, bry(1), tol)
      IF (nw /= 0) GO TO 10
    END IF
    cy(kdflg) = s2
    y(i) = s2 * csr(kflag)
    IF (kdflg == 2) GO TO 30
    kdflg = 2
    cycle
  END IF
 
  10 IF (rs1 > 0.0_dp) GO TO 150
!-----------------------------------------------------------------------
!     FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
  IF (x < 0.0_dp) GO TO 150
  kdflg = 1
  y(i) = czero
  nz = nz + 1
  IF (i /= 1) THEN
    IF (y(i-1) /= czero) THEN
      y(i-1) = czero
      nz = nz + 1
    END IF
  END IF
END DO
i = n
 
30 rz = 2.0_dp / zr
ck = fn * rz
ib = i + 1
IF (n >= ib) THEN
!-----------------------------------------------------------------------
!     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
  CALL cunik(zr, fn, 2, ipard, tol, initd, phid, zeta1d, zeta2d, sumd,  &
             cwrk(1:,3))
  IF (kode /= 1) THEN
    cfn = fn
    s1 = zeta1d - cfn * (cfn/(zr + zeta2d))
  ELSE
    s1 = zeta1d - zeta2d
  END IF
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) <= elim) THEN
    IF (abs(rs1) < alim) GO TO 50
!-----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!-----------------------------------------------------------------------
    aphi = abs(phid)
    rs1 = rs1 + log(aphi)
    IF (abs(rs1) < elim) GO TO 50
  END IF
  IF (rs1 > 0.0_dp) GO TO 150
!-----------------------------------------------------------------------
!     FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
  IF (x < 0.0_dp) GO TO 150
  nz = n
  y(1:n) = czero
  RETURN
!-----------------------------------------------------------------------
!     RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
!-----------------------------------------------------------------------
  50 s1 = cy(1)
  s2 = cy(2)
  c1 = csr(kflag)
  ascle = bry(kflag)
  DO  i = ib, n
    c2 = s2
    s2 = ck * s2 + s1
    s1 = c2
    ck = ck + rz
    c2 = s2 * c1
    y(i) = c2
    IF (kflag < 3) THEN
      c2r = real(c2, kind=dp)
      c2i = aimag(c2)
      c2r = abs(c2r)
      c2i = abs(c2i)
      c2m = max(c2r,c2i)
      IF (c2m > ascle) THEN
        kflag = kflag + 1
        ascle = bry(kflag)
        s1 = s1 * c1
        s2 = c2
        s1 = s1 * css(kflag)
        s2 = s2 * css(kflag)
        c1 = csr(kflag)
      END IF
    END IF
  END DO
END IF
IF (mr == 0) RETURN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z) < 0.0_dp
!-----------------------------------------------------------------------
nz = 0
fmr = mr
sgn = -sign(pi, fmr)
!-----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
!-----------------------------------------------------------------------
csgn = cmplx(0.0_dp, sgn, kind=dp)
inu = int(fnu)
fnf = fnu - inu
ifn = inu + n - 1
ang = fnf * sgn
cpn = cos(ang)
spn = sin(ang)
cspn = cmplx(cpn, spn, kind=dp)
IF (mod(ifn,2) == 1) cspn = -cspn
asc = bry(1)
kk = n
iuf = 0
kdflg = 1
ib = ib - 1
ic = ib - 1
DO  k = 1, n
  fn = fnu + (kk-1)
!-----------------------------------------------------------------------
!     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
!     FUNCTION ABOVE
!-----------------------------------------------------------------------
  m = 3
  IF (n > 2) GO TO 80
 
  70 initd = init(j)
  phid = phi(j)
  zeta1d = zeta1(j)
  zeta2d = zeta2(j)
  sumd = sum(j)
  m = j
  j = 3 - j
  GO TO 90
 
  80 IF (.NOT.(kk == n .AND. ib < n)) THEN
    IF (kk == ib .OR. kk == ic) GO TO 70
    initd = 0
  END IF
 
  90 CALL cunik(zr, fn, 1, 0, tol, initd, phid, zeta1d, zeta2d, sumd,  &
                cwrk(1:,m))
  IF (kode /= 1) THEN
    cfn = fn
    s1 = -zeta1d + cfn * (cfn/(zr + zeta2d))
  ELSE
    s1 = -zeta1d + zeta2d
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) > elim) GO TO 110
  IF (kdflg == 1) iflag = 2
  IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
    aphi = abs(phid)
    rs1 = rs1 + log(aphi)
    IF (abs(rs1) > elim) GO TO 110
    IF (kdflg == 1) iflag = 1
    IF (rs1 >= 0.0_dp) THEN
      IF (kdflg == 1) iflag = 3
    END IF
  END IF
  s2 = csgn * phid * sumd
  c2r = real(s1, kind=dp)
  c2i = aimag(s1)
  c2m = exp(c2r) * real(css(iflag), kind=dp)
  s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
  s2 = s2 * s1
  IF (iflag == 1) THEN
    CALL cuchk(s2, nw, bry(1), tol)
    IF (nw /= 0) s2 = 0.0_dp
  END IF
 
  100 cy(kdflg) = s2
  c2 = s2
  s2 = s2 * csr(iflag)
!-----------------------------------------------------------------------
!     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!-----------------------------------------------------------------------
  s1 = y(kk)
  IF (kode /= 1) THEN
    CALL cs1s2(zr, s1, s2, nw, asc, alim, iuf)
    nz = nz + nw
  END IF
  y(kk) = s1 * cspn + s2
  kk = kk - 1
  cspn = -cspn
  IF (c2 == czero) THEN
    kdflg = 1
    cycle
  END IF
  IF (kdflg == 2) GO TO 130
  kdflg = 2
  cycle
 
  110 IF (rs1 > 0.0_dp) GO TO 150
  s2 = czero
  GO TO 100
END DO
k = n
 
130 il = n - k
IF (il == 0) RETURN
!-----------------------------------------------------------------------
!     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(1)
s2 = cy(2)
cs = csr(iflag)
ascle = bry(iflag)
fn = inu + il
DO  i = 1, il
  c2 = s2
  s2 = s1 + (fn + fnf) * rz * s2
  s1 = c2
  fn = fn - 1.0_dp
  c2 = s2 * cs
  ck = c2
  c1 = y(kk)
  IF (kode /= 1) THEN
    CALL cs1s2(zr, c1, c2, nw, asc, alim, iuf)
    nz = nz + nw
  END IF
  y(kk) = c1 * cspn + c2
  kk = kk - 1
  cspn = -cspn
  IF (iflag < 3) THEN
    c2r = real(ck, kind=dp)
    c2i = aimag(ck)
    c2r = abs(c2r)
    c2i = abs(c2i)
    c2m = max(c2r, c2i)
    IF (c2m > ascle) THEN
      iflag = iflag + 1
      ascle = bry(iflag)
      s1 = s1 * cs
      s2 = ck
      s1 = s1 * css(iflag)
      s2 = s2 * css(iflag)
      cs = csr(iflag)
    END IF
  END IF
END DO
RETURN
 
150 nz = -1
RETURN
END SUBROUTINE cunk1
 
 
 
SUBROUTINE cunk2(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
!***BEGIN PROLOGUE  CUNK2
!***REFER TO  CBESK
 
!  CUNK2 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 INDICATES THE DIRECTION OF ROTATION FOR
!  ANALYTIC CONTINUATION.
!  NZ=-1 MEANS AN OVERFLOW WILL OCCUR
 
!***ROUTINES CALLED  CAIRY,CS1S2,CUCHK,CUNHJ,R1MACH
!***END PROLOGUE  CUNK2
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: mr
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: ai, arg(2), asum(2), bsum(2), cfn, ck, cs, csgn, cspn,  &
                 csr(3), css(3), cy(2), c1, c2, dai, phi(2), rz, s1, s2, &
                 zb, zeta1(2), zeta2(2), zn, zr, phid, argd, zeta1d, zeta2d, &
                 asumd, bsumd
real(dp)     :: aarg, ang, aphi, asc, ascle, bry(3), car, cpn, c2i,  &
                 c2m, c2r, crsc, cscl, fmr, fn, fnf, rs1, sar, sgn, spn, x, yy
INTEGER       :: i, ib, iflag, ifn, il, in, inu, iuf, k, kdflg, kflag, kk,  &
                 nai, ndai, nw, idum, j, ipard, ic
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp),  &
                            ci = (0.0_dp,1.0_dp),   &
                            cr1 = (1.0_dp, 1.73205080756887729_dp),  &
                            cr2 = (-0.5_dp, -8.66025403784438647d-01)
real(dp), PARAMETER     :: hpi = 1.57079632679489662_dp,  &
                            pi = 3.14159265358979324_dp,  &
                            aic = 1.26551212348464539_dp
COMPLEX (dp), PARAMETER  :: cip(4) = (/ (1.0_dp,0.0_dp), (0.0_dp,-1.0_dp),  &
                                        (-1.0_dp,0.0_dp), (0.0_dp,1.0_dp) /)
 
kdflg = 1
nz = 0
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
!     THE UNDERFLOW LIMIT
!-----------------------------------------------------------------------
cscl = 1.0_dp/tol
crsc = tol
css(1) = cscl
css(2) = cone
css(3) = crsc
csr(1) = crsc
csr(2) = cone
csr(3) = cscl
bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
bry(2) = 1.0_dp / bry(1)
bry(3) = huge(0.0_dp)
x = real(z, kind=dp)
zr = z
IF (x < 0.0_dp) zr = -z
yy = aimag(zr)
zn = -zr * ci
zb = zr
inu = int(fnu)
fnf = fnu - inu
ang = -hpi * fnf
car = cos(ang)
sar = sin(ang)
cpn = -hpi * car
spn = -hpi * sar
c2 = cmplx(-spn, cpn, kind=dp)
kk = mod(inu,4) + 1
cs = cr1 * c2 * cip(kk)
IF (yy <= 0.0_dp) THEN
  zn = conjg(-zn)
  zb = conjg(zb)
END IF
!-----------------------------------------------------------------------
!     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
DO  i = 1, n
!-----------------------------------------------------------------------
!     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!-----------------------------------------------------------------------
  j = 3 - j
  fn = fnu + (i-1)
  CALL cunhj(zn, fn, 0, tol, phi(j), arg(j), zeta1(j), zeta2(j), asum(j), &
             bsum(j))
  IF (kode /= 1) THEN
    cfn = cmplx(fn, 0.0_dp, kind=dp)
    s1 = zeta1(j) - cfn * (cfn/(zb + zeta2(j)))
  ELSE
    s1 = zeta1(j) - zeta2(j)
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) <= elim) THEN
    IF (kdflg == 1) kflag = 2
    IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
      aphi = abs(phi(j))
      aarg = abs(arg(j))
      rs1 = rs1 + log(aphi) - 0.25_dp * log(aarg) - aic
      IF (abs(rs1) > elim) GO TO 10
      IF (kdflg == 1) kflag = 1
      IF (rs1 >= 0.0_dp) THEN
        IF (kdflg == 1) kflag = 3
      END IF
    END IF
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
    c2 = arg(j) * cr2
    CALL cairy(c2, 0, 2, ai, nai, idum)
    CALL cairy(c2, 1, 2, dai, ndai, idum)
    s2 = cs * phi(j) * (ai*asum(j) + cr2*dai*bsum(j))
    c2r = real(s1, kind=dp)
    c2i = aimag(s1)
    c2m = exp(c2r) * real(css(kflag), kind=dp)
    s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
    s2 = s2 * s1
    IF (kflag == 1) THEN
      CALL cuchk(s2, nw, bry(1), tol)
      IF (nw /= 0) GO TO 10
    END IF
    IF (yy <= 0.0_dp) s2 = conjg(s2)
    cy(kdflg) = s2
    y(i) = s2 * csr(kflag)
    cs = -ci * cs
    IF (kdflg == 2) GO TO 30
    kdflg = 2
    cycle
  END IF
 
  10 IF (rs1 > 0.0_dp) GO TO 150
!-----------------------------------------------------------------------
!     FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
  IF (x < 0.0_dp) GO TO 150
  kdflg = 1
  y(i) = czero
  cs = -ci * cs
  nz = nz + 1
  IF (i /= 1) THEN
    IF (y(i-1) /= czero) THEN
      y(i-1) = czero
      nz = nz + 1
    END IF
  END IF
END DO
i = n
 
30 rz = 2.0_dp / zr
ck = fn * rz
ib = i + 1
IF (n >= ib) THEN
!-----------------------------------------------------------------------
!     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
!     ON UNDERFLOW
!-----------------------------------------------------------------------
  fn = fnu + (n-1)
  ipard = 1
  IF (mr /= 0) ipard = 0
  CALL cunhj(zn, fn, ipard, tol, phid, argd, zeta1d, zeta2d, asumd, bsumd)
  IF (kode /= 1) THEN
    cfn = fn
    s1 = zeta1d - cfn * (cfn/(zb + zeta2d))
  ELSE
    s1 = zeta1d - zeta2d
  END IF
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) <= elim) THEN
    IF (abs(rs1) < alim) GO TO 50
!-----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!-----------------------------------------------------------------------
    aphi = abs(phid)
    aarg = abs(argd)
    rs1 = rs1 + log(aphi) - 0.25_dp * log(aarg) - aic
    IF (abs(rs1) < elim) GO TO 50
  END IF
  IF (rs1 > 0.0_dp) GO TO 150
!-----------------------------------------------------------------------
!     FOR X < 0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
  IF (x < 0.0_dp) GO TO 150
  nz = n
  y(1:n) = czero
  RETURN
!-----------------------------------------------------------------------
!     SCALED FORWARD RECURRENCE FOR REMAINDER OF THE SEQUENCE
!-----------------------------------------------------------------------
  50 s1 = cy(1)
  s2 = cy(2)
  c1 = csr(kflag)
  ascle = bry(kflag)
  DO  i = ib, n
    c2 = s2
    s2 = ck * s2 + s1
    s1 = c2
    ck = ck + rz
    c2 = s2 * c1
    y(i) = c2
    IF (kflag < 3) THEN
      c2r = real(c2, kind=dp)
      c2i = aimag(c2)
      c2r = abs(c2r)
      c2i = abs(c2i)
      c2m = max(c2r,c2i)
      IF (c2m > ascle) THEN
        kflag = kflag + 1
        ascle = bry(kflag)
        s1 = s1 * c1
        s2 = c2
        s1 = s1 * css(kflag)
        s2 = s2 * css(kflag)
        c1 = csr(kflag)
      END IF
    END IF
  END DO
END IF
IF (mr == 0) RETURN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z) < 0.0_dp
!-----------------------------------------------------------------------
nz = 0
fmr = mr
sgn = -sign(pi, fmr)
!-----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
!-----------------------------------------------------------------------
csgn = cmplx(0.0_dp, sgn, kind=dp)
IF (yy <= 0.0_dp) csgn = conjg(csgn)
ifn = inu + n - 1
ang = fnf * sgn
cpn = cos(ang)
spn = sin(ang)
cspn = cmplx(cpn, spn, kind=dp)
IF (mod(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 = cmplx(car, -sar, kind=dp) * csgn
in = mod(ifn,4) + 1
c2 = cip(in)
cs = cs * conjg(c2)
asc = bry(1)
kk = n
kdflg = 1
ib = ib - 1
ic = ib - 1
iuf = 0
DO  k = 1, n
!-----------------------------------------------------------------------
!     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
!     FUNCTION ABOVE
!-----------------------------------------------------------------------
  fn = fnu + (kk-1)
  IF (n > 2) GO TO 80
 
  70 phid = phi(j)
  argd = arg(j)
  zeta1d = zeta1(j)
  zeta2d = zeta2(j)
  asumd = asum(j)
  bsumd = bsum(j)
  j = 3 - j
  GO TO 90
 
  80 IF (.NOT.(kk == n .AND. ib < n)) THEN
    IF (kk == ib .OR. kk == ic) GO TO 70
    CALL cunhj(zn, fn, 0, tol, phid, argd, zeta1d, zeta2d, asumd, bsumd)
  END IF
 
  90 IF (kode /= 1) THEN
    cfn = fn
    s1 = -zeta1d + cfn * (cfn/(zb + zeta2d))
  ELSE
    s1 = -zeta1d + zeta2d
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) > elim) GO TO 110
  IF (kdflg == 1) iflag = 2
  IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
    aphi = abs(phid)
    aarg = abs(argd)
    rs1 = rs1 + log(aphi) - 0.25_dp * log(aarg) - aic
    IF (abs(rs1) > elim) GO TO 110
    IF (kdflg == 1) iflag = 1
    IF (rs1 >= 0.0_dp) THEN
      IF (kdflg == 1) iflag = 3
    END IF
  END IF
  CALL cairy(argd, 0, 2, ai, nai, idum)
  CALL cairy(argd, 1, 2, dai, ndai, idum)
  s2 = cs * phid * (ai*asumd + dai*bsumd)
  c2r = real(s1, kind=dp)
  c2i = aimag(s1)
  c2m = exp(c2r) * real(css(iflag), kind=dp)
  s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
  s2 = s2 * s1
  IF (iflag == 1) THEN
    CALL cuchk(s2, nw, bry(1), tol)
    IF (nw /= 0) s2 = 0.0_dp
  END IF
 
  100 IF (yy <= 0.0_dp) s2 = conjg(s2)
  cy(kdflg) = s2
  c2 = s2
  s2 = s2 * csr(iflag)
!-----------------------------------------------------------------------
!     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!-----------------------------------------------------------------------
  s1 = y(kk)
  IF (kode /= 1) THEN
    CALL cs1s2(zr, s1, s2, nw, asc, alim, iuf)
    nz = nz + nw
  END IF
  y(kk) = s1 * cspn + s2
  kk = kk - 1
  cspn = -cspn
  cs = -cs * ci
  IF (c2 == czero) THEN
    kdflg = 1
    cycle
  END IF
  IF (kdflg == 2) GO TO 130
  kdflg = 2
  cycle
 
  110 IF (rs1 > 0.0_dp) GO TO 150
  s2 = czero
  GO TO 100
END DO
k = n
 
130 il = n - k
IF (il == 0) RETURN
!-----------------------------------------------------------------------
!     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(1)
s2 = cy(2)
cs = csr(iflag)
ascle = bry(iflag)
fn = inu + il
DO  i = 1, il
  c2 = s2
  s2 = s1 + cmplx(fn+fnf, 0.0_dp, kind=dp) * rz * s2
  s1 = c2
  fn = fn - 1.0_dp
  c2 = s2 * cs
  ck = c2
  c1 = y(kk)
  IF (kode /= 1) THEN
    CALL cs1s2(zr, c1, c2, nw, asc, alim, iuf)
    nz = nz + nw
  END IF
  y(kk) = c1 * cspn + c2
  kk = kk - 1
  cspn = -cspn
  IF (iflag < 3) THEN
    c2r = real(ck, kind=dp)
    c2i = aimag(ck)
    c2r = abs(c2r)
    c2i = abs(c2i)
    c2m = max(c2r, c2i)
    IF (c2m > ascle) THEN
      iflag = iflag + 1
      ascle = bry(iflag)
      s1 = s1 * cs
      s2 = ck
      s1 = s1 * css(iflag)
      s2 = s2 * css(iflag)
      cs = csr(iflag)
    END IF
  END IF
END DO
RETURN
 
150 nz = -1
RETURN
END SUBROUTINE cunk2
 
 
 
SUBROUTINE cbuni(z, fnu, kode, n, y, nz, nui, nlast, fnul, tol, elim, alim)
!***BEGIN PROLOGUE  CBUNI
!***REFER TO  CBESI,CBESK
 
!   CBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE ABS(Z) > FNUL AND
!   FNU+N-1 < 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  CUNI1,CUNI2,R1MACH
!***END PROLOGUE  CBUNI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(IN)        :: nui
INTEGER, INTENT(OUT)       :: nlast
real(dp), INTENT(IN)      :: fnul
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cscl, cscr, cy(2), rz, st, s1, s2
real(dp)     :: ax, ay, dfnu, fnui, gnu, xx, yy, ascle, bry(3), str, sti, &
                 stm
INTEGER       :: i, iflag, iform, k, nl, nw
 
nz = 0
xx = real(z, kind=dp)
yy = aimag(z)
ax = abs(xx) * 1.73205080756887_dp
ay = abs(yy)
iform = 1
IF (ay > ax) iform = 2
IF (nui == 0) GO TO 40
fnui = nui
dfnu = fnu + (n-1)
gnu = dfnu + fnui
IF (iform /= 2) THEN
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!-----------------------------------------------------------------------
  CALL cuni1(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
ELSE
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(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
!-----------------------------------------------------------------------
  CALL cuni2(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
END IF
IF (nw >= 0) THEN
  IF (nw /= 0) GO TO 50
  ay = abs(cy(1))
!----------------------------------------------------------------------
!     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
!----------------------------------------------------------------------
  bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
  bry(2) = 1.0_dp / bry(1)
  bry(3) = bry(2)
  iflag = 2
  ascle = bry(2)
  ax = 1.0_dp
  cscl = ax
  IF (ay <= bry(1)) THEN
    iflag = 1
    ascle = bry(1)
    ax = 1.0_dp / tol
    cscl = ax
  ELSE
    IF (ay >= bry(2)) THEN
      iflag = 3
      ascle = bry(3)
      ax = tol
      cscl = ax
    END IF
  END IF
  ay = 1.0_dp / ax
  cscr = ay
  s1 = cy(2) * cscl
  s2 = cy(1) * cscl
  rz = 2.0_dp / z
  DO  i = 1, nui
    st = s2
    s2 = cmplx(dfnu+fnui, 0.0_dp, kind=dp) * rz * s2 + s1
    s1 = st
    fnui = fnui - 1.0_dp
    IF (iflag < 3) THEN
      st = s2 * cscr
      str = real(st, kind=dp)
      sti = aimag(st)
      str = abs(str)
      sti = abs(sti)
      stm = max(str,sti)
      IF (stm > ascle) THEN
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1 * cscr
        s2 = st
        ax = ax * tol
        ay = 1.0_dp / ax
        cscl = ax
        cscr = ay
        s1 = s1 * cscl
        s2 = s2 * cscl
      END IF
    END IF
  END DO
  y(n) = s2 * cscr
  IF (n == 1) RETURN
  nl = n - 1
  fnui = nl
  k = nl
  DO  i = 1, nl
    st = s2
    s2 = cmplx(fnu+fnui, 0.0_dp, kind=dp) * rz * s2 + s1
    s1 = st
    st = s2 * cscr
    y(k) = st
    fnui = fnui - 1.0_dp
    k = k - 1
    IF (iflag < 3) THEN
      str = real(st, kind=dp)
      sti = aimag(st)
      str = abs(str)
      sti = abs(sti)
      stm = max(str,sti)
      IF (stm > ascle) THEN
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1 * cscr
        s2 = st
        ax = ax * tol
        ay = 1.0_dp / ax
        cscl = ax
        cscr = ay
        s1 = s1 * cscl
        s2 = s2 * cscl
      END IF
    END IF
  END DO
  RETURN
END IF
 
30 nz = -1
IF (nw == -2) nz = -2
RETURN
 
40 IF (iform /= 2) THEN
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3 <= ARG(Z) <= PI/3
!-----------------------------------------------------------------------
  CALL cuni1(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
ELSE
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(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
!-----------------------------------------------------------------------
  CALL cuni2(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
END IF
IF (nw < 0) GO TO 30
nz = nw
RETURN
 
50 nlast = n
RETURN
END SUBROUTINE cbuni
 
 
 
SUBROUTINE cuni1(z, fnu, kode, n, y, nz, nlast, fnul, tol, elim, alim)
!***BEGIN PROLOGUE  CUNI1
!***REFER TO  CBESI,CBESK
 
!     CUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
!     EXPANSION FOR I(FNU,Z) IN -PI/3 <= ARG Z <= 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 < FNUL.
!     Y(I)=CZERO FOR I=NLAST+1,N
 
!***ROUTINES CALLED  CUCHK,CUNIK,CUOIK,R1MACH
!***END PROLOGUE  CUNI1
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: nlast
real(dp), INTENT(IN)      :: fnul
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cfn, crsc, cscl, csr(3), css(3), c1, c2, cwrk(16), phi,  &
                 rz, sum, s1, s2, zeta1, zeta2, cy(2)
real(dp)     :: aphi, ascle, bry(3), c2i, c2m, c2r, fn, rs1, yy
INTEGER       :: i, iflag, init, k, m, nd, nn, nuf, nw
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
 
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 = cmplx(1.0_dp/tol, 0.0_dp, kind=dp)
crsc = tol
css(1) = cscl
css(2) = cone
css(3) = crsc
csr(1) = crsc
csr(2) = cone
csr(3) = cscl
bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
!-----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!-----------------------------------------------------------------------
fn = max(fnu, 1.0_dp)
init = 0
CALL cunik(z, fn, 1, 1, tol, init, phi, zeta1, zeta2, sum, cwrk)
IF (kode /= 1) THEN
  cfn = fn
  s1 = -zeta1 + cfn * (cfn/(z + zeta2))
ELSE
  s1 = -zeta1 + zeta2
END IF
rs1 = real(s1, kind=dp)
IF (abs(rs1) > elim) GO TO 70
 
10 nn = min(2,nd)
DO  i = 1, nn
  fn = fnu + (nd-i)
  init = 0
  CALL cunik(z, fn, 1, 0, tol, init, phi, zeta1, zeta2, sum, cwrk)
  IF (kode /= 1) THEN
    cfn = fn
    yy = aimag(z)
    s1 = -zeta1 + cfn * (cfn/(z+zeta2)) + cmplx(0.0_dp, yy, kind=dp)
  ELSE
    s1 = -zeta1 + zeta2
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) > elim) GO TO 50
  IF (i == 1) iflag = 2
  IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
    aphi = abs(phi)
    rs1 = rs1 + log(aphi)
    IF (abs(rs1) > elim) GO TO 50
    IF (i == 1) iflag = 1
    IF (rs1 >= 0.0_dp) THEN
      IF (i == 1) iflag = 3
    END IF
  END IF
!-----------------------------------------------------------------------
!     SCALE S1 IF ABS(S1) < ASCLE
!-----------------------------------------------------------------------
  s2 = phi * sum
  c2r = real(s1, kind=dp)
  c2i = aimag(s1)
  c2m = exp(c2r) * real(css(iflag), kind=dp)
  s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
  s2 = s2 * s1
  IF (iflag == 1) THEN
    CALL cuchk(s2, nw, bry(1), tol)
    IF (nw /= 0) GO TO 50
  END IF
  m = nd - i + 1
  cy(i) = s2
  y(m) = s2 * csr(iflag)
END DO
IF (nd > 2) THEN
  rz = 2.0_dp / z
  bry(2) = 1.0_dp / bry(1)
  bry(3) = huge(0.0_dp)
  s1 = cy(1)
  s2 = cy(2)
  c1 = csr(iflag)
  ascle = bry(iflag)
  k = nd - 2
  fn = k
  DO  i = 3, nd
    c2 = s2
    s2 = s1 + cmplx(fnu+fn, 0.0_dp, kind=dp) * rz * s2
    s1 = c2
    c2 = s2 * c1
    y(k) = c2
    k = k - 1
    fn = fn - 1.0_dp
    IF (iflag < 3) THEN
      c2r = real(c2, kind=dp)
      c2i = aimag(c2)
      c2r = abs(c2r)
      c2i = abs(c2i)
      c2m = max(c2r,c2i)
      IF (c2m > ascle) THEN
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1 * c1
        s2 = c2
        s1 = s1 * css(iflag)
        s2 = s2 * css(iflag)
        c1 = csr(iflag)
      END IF
    END IF
  END DO
END IF
 
40 RETURN
!-----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!-----------------------------------------------------------------------
50 IF (rs1 <= 0.0_dp) THEN
  y(nd) = czero
  nz = nz + 1
  nd = nd - 1
  IF (nd == 0) GO TO 40
  CALL cuoik(z, fnu, kode, 1, nd, y, nuf, tol, elim, alim)
  IF (nuf >= 0) THEN
    nd = nd - nuf
    nz = nz + nuf
    IF (nd == 0) GO TO 40
    fn = fnu + (nd-1)
    IF (fn >= fnul) GO TO 10
    nlast = nd
    RETURN
  END IF
END IF
 
60 nz = -1
RETURN
 
70 IF (rs1 > 0.0_dp) GO TO 60
nz = n
y(1:n) = czero
RETURN
END SUBROUTINE cuni1
 
 
 
SUBROUTINE cuni2(z, fnu, kode, n, y, nz, nlast, fnul, tol, elim, alim)
!***BEGIN PROLOGUE  CUNI2
!***REFER TO  CBESI,CBESK
 
!     CUNI2 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 < FNUL.
!     Y(I) = CZERO FOR I=NLAST+1,N
 
!***ROUTINES CALLED  CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH
!***END PROLOGUE  CUNI2
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
INTEGER, INTENT(OUT)       :: nlast
real(dp), INTENT(IN)      :: fnul
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: ai, arg, asum, bsum, cfn, cid, crsc, cscl, csr(3), css(3), &
                 cy(2), c1, c2, dai, phi, rz, s1, s2, zb, zeta1, zeta2,  &
                 zn, zar
real(dp)     :: aarg, ang, aphi, ascle, ay, bry(3), car, c2i, c2m,  &
                 c2r, fn, rs1, sar, yy
INTEGER       :: i, iflag, in, inu, j, k, nai, nd, ndai, nn, nuf, nw, idum
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp), &
                            ci = (0.0_dp,1.0_dp)
COMPLEX (dp), PARAMETER  :: cip(4) = (/ (1.0_dp,0.0_dp), (0.0_dp,1.0_dp),  &
                                        (-1.0_dp,0.0_dp), (0.0_dp,-1.0_dp) /)
real(dp), PARAMETER     :: hpi = 1.57079632679489662_dp, aic = 1.265512123484645396_dp
 
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 = cmplx(1.0_dp/tol, 0.0_dp, kind=dp)
crsc = tol
css(1) = cscl
css(2) = cone
css(3) = crsc
csr(1) = crsc
csr(2) = cone
csr(3) = cscl
bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
yy = aimag(z)
!-----------------------------------------------------------------------
!     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)
car = cos(ang)
sar = sin(ang)
c2 = cmplx(car, sar, kind=dp)
zar = c2
in = inu + n - 1
in = mod(in,4)
c2 = c2 * cip(in+1)
IF (yy <= 0.0_dp) THEN
  zn = conjg(-zn)
  zb = conjg(zb)
  cid = -cid
  c2 = conjg(c2)
END IF
!-----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!-----------------------------------------------------------------------
fn = max(fnu,1.0_dp)
CALL cunhj(zn, fn, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
IF (kode /= 1) THEN
  cfn = fnu
  s1 = -zeta1 + cfn * (cfn/(zb + zeta2))
ELSE
  s1 = -zeta1 + zeta2
END IF
rs1 = real(s1, kind=dp)
IF (abs(rs1) > elim) GO TO 70
 
10 nn = min(2,nd)
DO  i = 1, nn
  fn = fnu + (nd-i)
  CALL cunhj(zn, fn, 0, tol, phi, arg, zeta1, zeta2, asum, bsum)
  IF (kode /= 1) THEN
    cfn = fn
    ay = abs(yy)
    s1 = -zeta1 + cfn * (cfn/(zb+zeta2)) + cmplx(0.0_dp, ay, kind=dp)
  ELSE
    s1 = -zeta1 + zeta2
  END IF
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
  rs1 = real(s1, kind=dp)
  IF (abs(rs1) > elim) GO TO 50
  IF (i == 1) iflag = 2
  IF (abs(rs1) >= alim) THEN
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
    aphi = abs(phi)
    aarg = abs(arg)
    rs1 = rs1 + log(aphi) - 0.25_dp * log(aarg) - aic
    IF (abs(rs1) > elim) GO TO 50
    IF (i == 1) iflag = 1
    IF (rs1 >= 0.0_dp) THEN
      IF (i == 1) iflag = 3
    END IF
  END IF
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
  CALL cairy(arg, 0, 2, ai, nai, idum)
  CALL cairy(arg, 1, 2, dai, ndai, idum)
  s2 = phi * (ai*asum + dai*bsum)
  c2r = real(s1, kind=dp)
  c2i = aimag(s1)
  c2m = exp(c2r) * real(css(iflag), kind=dp)
  s1 = c2m * cmplx(cos(c2i), sin(c2i), kind=dp)
  s2 = s2 * s1
  IF (iflag == 1) THEN
    CALL cuchk(s2, nw, bry(1), tol)
    IF (nw /= 0) GO TO 50
  END IF
  IF (yy <= 0.0_dp) s2 = conjg(s2)
  j = nd - i + 1
  s2 = s2 * c2
  cy(i) = s2
  y(j) = s2 * csr(iflag)
  c2 = c2 * cid
END DO
IF (nd > 2) THEN
  rz = 2.0_dp / z
  bry(2) = 1.0_dp / bry(1)
  bry(3) = huge(0.0_dp)
  s1 = cy(1)
  s2 = cy(2)
  c1 = csr(iflag)
  ascle = bry(iflag)
  k = nd - 2
  fn = k
  DO  i = 3, nd
    c2 = s2
    s2 = s1 + (fnu + fn) * rz * s2
    s1 = c2
    c2 = s2 * c1
    y(k) = c2
    k = k - 1
    fn = fn - 1.0_dp
    IF (iflag < 3) THEN
      c2r = real(c2, kind=dp)
      c2i = aimag(c2)
      c2r = abs(c2r)
      c2i = abs(c2i)
      c2m = max(c2r,c2i)
      IF (c2m > ascle) THEN
        iflag = iflag + 1
        ascle = bry(iflag)
        s1 = s1 * c1
        s2 = c2
        s1 = s1 * css(iflag)
        s2 = s2 * css(iflag)
        c1 = csr(iflag)
      END IF
    END IF
  END DO
END IF
 
40 RETURN
 
50 IF (rs1 <= 0.0_dp) THEN
!-----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!-----------------------------------------------------------------------
  y(nd) = czero
  nz = nz + 1
  nd = nd - 1
  IF (nd == 0) GO TO 40
  CALL cuoik(z, fnu, kode, 1, nd, y, nuf, tol, elim, alim)
  IF (nuf >= 0) THEN
    nd = nd - nuf
    nz = nz + nuf
    IF (nd == 0) GO TO 40
    fn = fnu + (nd-1)
    IF (fn >= fnul) THEN
!      FN = AIMAG(CID)
!      J = NUF + 1
!      K = MOD(J,4) + 1
!      S1 = CIP(K)
!      IF (FN < 0.0_dp) S1 = CONJG(S1)
!      C2 = C2*S1
      in = inu + nd - 1
      in = mod(in,4) + 1
      c2 = zar * cip(in)
      IF (yy <= 0.0_dp) c2 = conjg(c2)
      GO TO 10
    END IF
    nlast = nd
    RETURN
  END IF
END IF
 
60 nz = -1
RETURN
 
70 IF (rs1 > 0.0_dp) GO TO 60
nz = n
y(1:n) = czero
RETURN
END SUBROUTINE cuni2
 
 
 
SUBROUTINE cs1s2(zr, s1, s2, nz, ascle, alim, iuf)
!***BEGIN PROLOGUE  CS1S2
!***REFER TO  CBESK,CAIRY
 
!     CS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE ADDITION OF
!     THE I AND K FUNCTIONS IN THE ANALYTIC CONTINUATION 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  (NONE)
!***END PROLOGUE  CS1S2
 
COMPLEX (dp), INTENT(IN)      :: zr
COMPLEX (dp), INTENT(IN OUT)  :: s1
COMPLEX (dp), INTENT(IN OUT)  :: s2
INTEGER, INTENT(OUT)          :: nz
real(dp), INTENT(IN)         :: ascle
real(dp), INTENT(IN)         :: alim
INTEGER, INTENT(IN OUT)       :: iuf
 
COMPLEX (dp)  :: c1, s1d
real(dp)     :: aa, aln, as1, as2, xx
 
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp)
 
nz = 0
as1 = abs(s1)
as2 = abs(s2)
aa = real(s1, kind=dp)
aln = aimag(s1)
IF (aa /= 0.0_dp .OR. aln /= 0.0_dp) THEN
  IF (as1 /= 0.0_dp) THEN
    xx = real(zr, kind=dp)
    aln = -xx - xx + log(as1)
    s1d = s1
    s1 = czero
    as1 = 0.0_dp
    IF (aln >= -alim) THEN
      c1 = log(s1d) - zr - zr
      s1 = exp(c1)
      as1 = abs(s1)
      iuf = iuf + 1
    END IF
  END IF
END IF
aa = max(as1,as2)
IF (aa > ascle) RETURN
s1 = czero
s2 = czero
nz = 1
iuf = 0
RETURN
END SUBROUTINE cs1s2
 
 
 
SUBROUTINE cshch(z, csh, cch)
!***BEGIN PROLOGUE  CSHCH
!***REFER TO  CBESK,CBESH
 
!     CSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+I*Y)
!     AND CCH=COSH(X+I*Y), WHERE I**2=-1.
 
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  CSHCH
 
COMPLEX (dp), INTENT(IN OUT)  :: z
COMPLEX (dp), INTENT(OUT)     :: csh
COMPLEX (dp), INTENT(OUT)     :: cch
 
real(dp)  :: cchi, cchr, ch, cn, cshi, cshr, sh, sn, x, y
 
x = real(z, kind=dp)
y = aimag(z)
sh = sinh(x)
ch = cosh(x)
sn = sin(y)
cn = cos(y)
cshr = sh * cn
cshi = ch * sn
csh = cmplx(cshr, cshi, kind=dp)
cchr = ch * cn
cchi = sh * sn
cch = cmplx(cchr, cchi, kind=dp)
RETURN
END SUBROUTINE cshch
 
 
 
SUBROUTINE crati(z, fnu, n, cy, tol)
!***BEGIN PROLOGUE  CRATI
!***REFER TO  CBESI,CBESK,CBESH
 
!   CRATI 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  (NONE)
!***END PROLOGUE  CRATI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
real(dp), INTENT(IN)      :: tol
 
COMPLEX (dp)  :: cdfnu, pt, p1, p2, rz, t1
real(dp)     :: ak, amagz, ap1, ap2, arg, az, dfnu, fdnu, flam, fnup,  &
                 rap1, rho, test, test1
INTEGER       :: i, id, idnu, inu, itime, k, kk, magz
 
real(dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp)
 
az = abs(z)
inu = int(fnu)
idnu = inu + n - 1
fdnu = idnu
magz = int(az)
amagz = magz + 1
fnup = max(amagz, fdnu)
id = idnu - magz - 1
itime = 1
k = 1
rz = (cone+cone) / z
t1 = fnup * rz
p2 = -t1
p1 = cone
t1 = t1 + rz
IF (id > 0) id = 0
ap2 = abs(p2)
ap1 = 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 = sqrt(arg)
test = test1
rap1 = 1.0_dp / ap1
p1 = p1 * rap1
p2 = p2 * rap1
ap2 = ap2 * rap1
 
10 k = k + 1
ap1 = ap2
pt = p2
p2 = p1 - t1 * p2
p1 = pt
t1 = t1 + rz
ap2 = abs(p2)
IF (ap1 <= test) GO TO 10
IF (itime /= 2) THEN
  ak = abs(t1) * 0.5_dp
  flam = ak + sqrt(ak*ak - 1.0_dp)
  rho = min(ap2/ap1, flam)
  test = test1 * sqrt(rho/(rho*rho - 1.0_dp))
  itime = 2
  GO TO 10
END IF
kk = k + 1 - id
ak = kk
dfnu = fnu + (n-1)
cdfnu = dfnu
t1 = ak
p1 = 1.0_dp/ap2
p2 = czero
DO  i = 1, kk
  pt = p1
  p1 = rz * (cdfnu+t1) * p1 + p2
  p2 = pt
  t1 = t1 - cone
END DO
IF (real(p1, kind=dp) == 0.0_dp .AND. aimag(p1) == 0.0_dp) THEN
  p1 = cmplx(tol, tol, kind=dp)
END IF
cy(n) = p2 / p1
IF (n == 1) RETURN
k = n - 1
ak = k
t1 = ak
cdfnu = fnu * rz
DO  i = 2, n
  pt = cdfnu + t1 * rz + cy(k+1)
  IF (real(pt, kind=dp) == 0.0_dp .AND. aimag(pt) == 0.0_dp) THEN
    pt = cmplx(tol, tol, kind=dp)
  END IF
  cy(k) = cone / pt
  t1 = t1 - cone
  k = k - 1
END DO
RETURN
END SUBROUTINE crati
 
 
 
SUBROUTINE cbknu(z, fnu, kode, n, y, nz, tol, elim, alim)
!***BEGIN PROLOGUE  CBKNU
!***REFER TO  CBESI,CBESK,CAIRY,CBESH
 
!     CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE
 
!***ROUTINES CALLED  CKSCL,CSHCH,GAMLN,I1MACH,R1MACH,CUCHK
!***END PROLOGUE  CBKNU
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cch, ck, coef, crsc, cs, cscl, csh, csr(3), css(3), cz,  &
                 f, fmu, p, pt, p1, p2, q, rz, smu, st, s1, s2, zd, celm, &
                 cy(2)
real(dp)     :: aa, ak, ascle, a1, a2, bb, bk, bry(3), caz, dnu, dnu2,  &
                 etest, fc, fhs, fk, fks, g1, g2, p2i, p2m, p2r, rk, s,  &
                 tm, t1, t2, xx, yy, helim, elm, xd, yd, alas, as
INTEGER       :: i, iflag, inu, k, kflag, kk, koded, nw, j, ic, inub
 
INTEGER, PARAMETER       :: kmax = 30
real(dp), PARAMETER     :: r1 = 2.0_dp
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp), cone = (1.0_dp,0.0_dp), &
                            ctwo = (2.0_dp,0.0_dp)
 
real(dp), PARAMETER  :: pi = 3.14159265358979324_dp,  &
    rthpi = 1.25331413731550025_dp, spi = 1.90985931710274403_dp,  &
    hpi = 1.57079632679489662_dp, fpi = 1.89769999331517738_dp,  &
    tth = 6.66666666666666666d-01
 
real(dp), PARAMETER  :: cc(8) = (/ 5.77215664901532861d-01,  &
 -4.20026350340952355d-02, -4.21977345555443367d-02, 7.21894324666309954d-03, &
 -2.15241674114950973d-04, -2.01348547807882387d-05, 1.13302723198169588d-06, &
  6.11609510448141582d-09 /)
 
xx = real(z, kind=dp)
yy = aimag(z)
caz = abs(z)
cscl = 1.0_dp/tol
crsc = tol
css(1) = cscl
css(2) = cone
css(3) = crsc
csr(1) = crsc
csr(2) = cone
csr(3) = cscl
bry(1) = 1.0e+3 * tiny(0.0_dp) / tol
bry(2) = 1.0_dp / bry(1)
bry(3) = huge(0.0_dp)
nz = 0
iflag = 0
koded = kode
rz = ctwo / z
inu = int(fnu + 0.5_dp)
dnu = fnu - inu
IF (abs(dnu) /= 0.5_dp) THEN
  dnu2 = 0.0_dp
  IF (abs(dnu) > tol) dnu2 = dnu * dnu
  IF (caz <= r1) THEN
!-----------------------------------------------------------------------
!     SERIES FOR ABS(Z) <= R1
!-----------------------------------------------------------------------
    fc = 1.0_dp
    smu = log(rz)
    fmu = smu * dnu
    CALL cshch(fmu, csh, cch)
    IF (dnu /= 0.0_dp) THEN
      fc = dnu * pi
      fc = fc / sin(fc)
      smu = csh / dnu
    END IF
    a2 = 1.0_dp + dnu
!-----------------------------------------------------------------------
!     GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
!-----------------------------------------------------------------------
    t2 = exp(-gamln(a2))
    t1 = 1.0_dp / (t2*fc)
    IF (abs(dnu) <= 0.1_dp) THEN
!-----------------------------------------------------------------------
!     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
!-----------------------------------------------------------------------
      ak = 1.0_dp
      s = cc(1)
      DO  k = 2, 8
        ak = ak * dnu2
        tm = cc(k) * ak
        s = s + tm
        IF (abs(tm) < tol) EXIT
      END DO
      g1 = -s
    ELSE
      g1 = (t1-t2) / (dnu+dnu)
    END IF
    g2 = 0.5_dp * (t1+t2) * fc
    g1 = g1 * fc
    f = g1 * cch + smu * g2
    pt = exp(fmu)
    p = cmplx(0.5_dp/t2, 0.0_dp, kind=dp) * pt
    q = cmplx(0.5_dp/t1, 0.0_dp, kind=dp) / pt
    s1 = f
    s2 = p
    ak = 1.0_dp
    a1 = 1.0_dp
    ck = cone
    bk = 1.0_dp - dnu2
    IF (inu <= 0 .AND. n <= 1) THEN
!-----------------------------------------------------------------------
!     GENERATE K(FNU,Z), 0.0D0  <=  FNU  <  0.5D0 AND N=1
!-----------------------------------------------------------------------
      IF (caz >= tol) THEN
        cz = z * z * 0.25_dp
        t1 = 0.25_dp * caz * caz
 
        30 f = (f*ak + p + q) / bk
        p = p / (ak-dnu)
        q = q / (ak+dnu)
        rk = 1.0_dp / ak
        ck = ck * cz * rk
        s1 = s1 + ck * f
        a1 = a1 * t1 * rk
        bk = bk + ak + ak + 1.0_dp
        ak = ak + 1.0_dp
        IF (a1 > tol) GO TO 30
      END IF
      y(1) = s1
      IF (koded == 1) RETURN
      y(1) = s1 * exp(z)
      RETURN
    END IF
!-----------------------------------------------------------------------
!     GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
!-----------------------------------------------------------------------
    IF (caz >= tol) THEN
      cz = z * z * 0.25_dp
      t1 = 0.25_dp * caz * caz
 
      40 f = (f*ak + p + q) / bk
      p = p / (ak-dnu)
      q = q / (ak+dnu)
      rk = 1.0_dp / ak
      ck = ck * cz * rk
      s1 = s1 + ck * f
      s2 = s2 + ck * (p - f*ak)
      a1 = a1 * t1 * rk
      bk = bk + ak + ak + 1.0_dp
      ak = ak + 1.0_dp
      IF (a1 > tol) GO TO 40
    END IF
    kflag = 2
    bk = real(smu, kind=dp)
    a1 = fnu + 1.0_dp
    ak = a1 * abs(bk)
    IF (ak > alim) kflag = 3
    p2 = s2 * css(kflag)
    s2 = p2 * rz
    s1 = s1 * css(kflag)
    IF (koded == 1) GO TO 100
    f = exp(z)
    s1 = s1 * f
    s2 = s2 * f
    GO TO 100
  END IF
END IF
!-----------------------------------------------------------------------
!     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 / sqrt(z)
kflag = 2
IF (koded /= 2) THEN
  IF (xx > alim) GO TO 200
!     BLANK LINE
  a1 = exp(-xx) * real(css(kflag), kind=dp)
  pt = a1 * cmplx(cos(yy), -sin(yy), kind=dp)
  coef = coef * pt
END IF
 
50 IF (abs(dnu) == 0.5_dp) GO TO 210
!-----------------------------------------------------------------------
!     MILLER ALGORITHM FOR ABS(Z) > R1
!-----------------------------------------------------------------------
ak = cos(pi*dnu)
ak = abs(ak)
IF (ak == 0.0_dp) GO TO 210
fhs = abs(0.25_dp - dnu2)
IF (fhs == 0.0_dp) GO TO 210
!-----------------------------------------------------------------------
!     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 = (digits(0.0_dp) - 1) * log10( real( radix(0.0_dp), kind=dp) ) * 3.321928094_dp
t1 = max(t1,12.0_dp)
t1 = min(t1,60.0_dp)
t2 = tth * t1 - 6.0_dp
IF (xx == 0.0_dp) THEN
  t1 = hpi
ELSE
  t1 = atan(yy/xx)
  t1 = abs(t1)
END IF
IF (t2 <= caz) THEN
!-----------------------------------------------------------------------
!     FORWARD RECURRENCE LOOP WHEN ABS(Z) >= R2
!-----------------------------------------------------------------------
  etest = ak / (pi*caz*tol)
  fk = 1.0_dp
  IF (etest < 1.0_dp) GO TO 80
  fks = 2.0_dp
  rk = caz + caz + 2.0_dp
  a1 = 0.0_dp
  a2 = 1.0_dp
  DO  i = 1, kmax
    ak = fhs / fks
    bk = rk / (fk+1.0_dp)
    tm = a2
    a2 = bk * a2 - ak * a1
    a1 = tm
    rk = rk + 2.0_dp
    fks = fks + fk + fk + 2.0_dp
    fhs = fhs + fk + fk
    fk = fk + 1.0_dp
    tm = abs(a2) * fk
    IF (etest < tm) GO TO 70
  END DO
  GO TO 220
 
  70 fk = fk + spi * t1 * sqrt(t2/caz)
  fhs = abs(0.25_dp-dnu2)
ELSE
!-----------------------------------------------------------------------
!     COMPUTE BACKWARD INDEX K FOR ABS(Z) < R2
!-----------------------------------------------------------------------
  a2 = sqrt(caz)
  ak = fpi * ak / (tol*sqrt(a2))
  aa = 3.0_dp * t1 / (1.0_dp+caz)
  bb = 14.7_dp * t1 / (28.0_dp+caz)
  ak = (log(ak) + caz*cos(aa)/(1.0_dp + 0.008_dp*caz)) / cos(bb)
  fk = 0.12125_dp * ak * ak / caz + 1.5_dp
END IF
 
80 k = int(fk)
!-----------------------------------------------------------------------
!     BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
!-----------------------------------------------------------------------
fk = k
fks = fk * fk
p1 = czero
p2 = tol
cs = p2
DO  i = 1, k
  a1 = fks - fk
  a2 = (fks+fk) / (a1+fhs)
  rk = 2.0_dp / (fk + 1.0_dp)
  t1 = (fk+xx) * rk
  t2 = yy * rk
  pt = p2
  p2 = (p2*cmplx(t1, t2, kind=dp) - p1) * a2
  p1 = pt
  cs = cs + p2
  fks = a1 - fk + 1.0_dp
  fk = fk - 1.0_dp
END DO
!-----------------------------------------------------------------------
!     COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER SCALING
!-----------------------------------------------------------------------
tm = abs(cs)
pt = cmplx(1.0_dp/tm, 0.0_dp, kind=dp)
s1 = pt * p2
cs = conjg(cs) * pt
s1 = coef * s1 * cs
IF (inu <= 0 .AND. n <= 1) THEN
  zd = z
  IF (iflag == 1) GO TO 190
  GO TO 130
END IF
!-----------------------------------------------------------------------
!     COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING
!-----------------------------------------------------------------------
tm = abs(p2)
pt = cmplx(1.0_dp/tm, 0.0_dp, kind=dp)
p1 = pt * p1
p2 = conjg(p2) * pt
pt = p1 * p2
s2 = s1 * (cone + (cmplx(dnu+0.5_dp, 0.0_dp, kind=dp) - pt)/z)
!-----------------------------------------------------------------------
!     FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
!     SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
!-----------------------------------------------------------------------
100 ck = cmplx(dnu+1.0_dp, 0.0_dp, kind=dp) * rz
IF (n == 1) inu = inu - 1
IF (inu <= 0) THEN
  IF (n == 1) s1 = s2
  zd = z
  IF (iflag == 1) GO TO 190
  GO TO 130
END IF
inub = 1
IF (iflag == 1) GO TO 160
 
110 p1 = csr(kflag)
ascle = bry(kflag)
DO  i = inub, inu
  st = s2
  s2 = ck * s2 + s1
  s1 = st
  ck = ck + rz
  IF (kflag < 3) THEN
    p2 = s2 * p1
    p2r = real(p2, kind=dp)
    p2i = aimag(p2)
    p2r = abs(p2r)
    p2i = abs(p2i)
    p2m = max(p2r,p2i)
    IF (p2m > ascle) THEN
      kflag = kflag + 1
      ascle = bry(kflag)
      s1 = s1 * p1
      s2 = p2
      s1 = s1 * css(kflag)
      s2 = s2 * css(kflag)
      p1 = csr(kflag)
    END IF
  END IF
END DO
IF (n == 1) s1 = s2
 
130 y(1) = s1 * csr(kflag)
IF (n == 1) RETURN
y(2) = s2 * csr(kflag)
IF (n == 2) RETURN
kk = 2
 
140 kk = kk + 1
IF (kk > n) RETURN
p1 = csr(kflag)
ascle = bry(kflag)
DO  i = kk, n
  p2 = s2
  s2 = ck * s2 + s1
  s1 = p2
  ck = ck + rz
  p2 = s2 * p1
  y(i) = p2
  IF (kflag < 3) THEN
    p2r = real(p2, kind=dp)
    p2i = aimag(p2)
    p2r = abs(p2r)
    p2i = abs(p2i)
    p2m = max(p2r,p2i)
    IF (p2m > ascle) THEN
      kflag = kflag + 1
      ascle = bry(kflag)
      s1 = s1 * p1
      s2 = p2
      s1 = s1 * css(kflag)
      s2 = s2 * css(kflag)
      p1 = csr(kflag)
    END IF
  END IF
END DO
RETURN
!-----------------------------------------------------------------------
!     IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
!-----------------------------------------------------------------------
160 helim = 0.5_dp * elim
elm = exp(-elim)
celm = elm
ascle = bry(1)
zd = z
xd = xx
yd = yy
ic = -1
j = 2
DO  i = 1, inu
  st = s2
  s2 = ck * s2 + s1
  s1 = st
  ck = ck + rz
  as = abs(s2)
  alas = log(as)
  p2r = -xd + alas
  IF (p2r >= -elim) THEN
    p2 = -zd + log(s2)
    p2r = real(p2, kind=dp)
    p2i = aimag(p2)
    p2m = exp(p2r) / tol
    p1 = p2m * cmplx(cos(p2i), sin(p2i), kind=dp)
    CALL cuchk(p1, nw, ascle, tol)
    IF (nw == 0) THEN
      j = 3 - j
      cy(j) = p1
      IF (ic == i-1) GO TO 180
      ic = i
      cycle
    END IF
  END IF
  IF (alas >= helim) THEN
    xd = xd - elim
    s1 = s1 * celm
    s2 = s2 * celm
    zd = cmplx(xd, yd, kind=dp)
  END IF
END DO
IF (n == 1) s1 = s2
GO TO 190
 
180 kflag = 1
inub = i + 1
s2 = cy(j)
j = 3 - j
s1 = cy(j)
IF (inub <= inu) GO TO 110
IF (n == 1) s1 = s2
GO TO 130
 
190 y(1) = s1
IF (n /= 1) THEN
  y(2) = s2
END IF
ascle = bry(1)
CALL ckscl(zd, fnu, n, y, nz, rz, ascle, tol, elim)
inu = n - nz
IF (inu <= 0) RETURN
kk = nz + 1
s1 = y(kk)
y(kk) = s1 * csr(1)
IF (inu == 1) RETURN
kk = nz + 2
s2 = y(kk)
y(kk) = s2 * csr(1)
IF (inu == 2) RETURN
t2 = fnu + (kk-1)
ck = t2 * rz
kflag = 1
GO TO 140
!-----------------------------------------------------------------------
!     SCALE BY EXP(Z), IFLAG = 1 CASES
!-----------------------------------------------------------------------
200 koded = 2
iflag = 1
kflag = 2
GO TO 50
!-----------------------------------------------------------------------
!     FNU=HALF ODD INTEGER CASE, DNU=-0.5
!-----------------------------------------------------------------------
210 s1 = coef
s2 = coef
GO TO 100
 
220 nz = -2
RETURN
END SUBROUTINE cbknu
 
 
 
SUBROUTINE ckscl(zr, fnu, n, y, nz, rz, ascle, tol, elim)
!***BEGIN PROLOGUE  CKSCL
!***REFER TO  CBKNU,CUNK1,CUNK2
 
!     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  CUCHK
!***END PROLOGUE  CKSCL
 
COMPLEX (dp), INTENT(IN)   :: zr
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
COMPLEX (dp), INTENT(IN)   :: rz
real(dp), INTENT(IN OUT)  :: ascle
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
 
COMPLEX (dp)  :: ck, cs, cy(2), s1, s2, zd, celm
real(dp)     :: aa, acs, as, csi, csr, fn, xx, zri, elm, alas, helim
INTEGER       :: i, ic, kk, nn, nw
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp,0.0_dp)
 
nz = 0
ic = 0
xx = real(zr, kind=dp)
nn = min(2,n)
DO  i = 1, nn
  s1 = y(i)
  cy(i) = s1
  as = abs(s1)
  acs = -xx + log(as)
  nz = nz + 1
  y(i) = czero
  IF (acs >= -elim) THEN
    cs = -zr + log(s1)
    csr = real(cs, kind=dp)
    csi = aimag(cs)
    aa = exp(csr) / tol
    cs = aa * cmplx(cos(csi), sin(csi), kind=dp)
    CALL cuchk(cs, nw, ascle, tol)
    IF (nw == 0) THEN
      y(i) = cs
      nz = nz - 1
      ic = i
    END IF
  END IF
END DO
IF (n == 1) RETURN
IF (ic <= 1) THEN
  y(1) = czero
  nz = 2
END IF
IF (n == 2) RETURN
IF (nz == 0) RETURN
fn = fnu + 1.0_dp
ck = fn * rz
s1 = cy(1)
s2 = cy(2)
helim = 0.5_dp * elim
elm = exp(-elim)
celm = elm
zri = aimag(zr)
zd = zr
 
!     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
!     S2 GETS LARGER THAN EXP(ELIM/2)
 
DO  i = 3, n
  kk = i
  cs = s2
  s2 = ck * s2 + s1
  s1 = cs
  ck = ck + rz
  as = abs(s2)
  alas = log(as)
  acs = -xx + alas
  nz = nz + 1
  y(i) = czero
  IF (acs >= -elim) THEN
    cs = -zd + log(s2)
    csr = real(cs, kind=dp)
    csi = aimag(cs)
    aa = exp(csr) / tol
    cs = aa * cmplx(cos(csi), sin(csi), kind=dp)
    CALL cuchk(cs, nw, ascle, tol)
    IF (nw == 0) THEN
      y(i) = cs
      nz = nz - 1
      IF (ic == kk-1) GO TO 30
      ic = kk
      cycle
    END IF
  END IF
  IF (alas >= helim) THEN
    xx = xx - elim
    s1 = s1 * celm
    s2 = s2 * celm
    zd = cmplx(xx, zri, kind=dp)
  END IF
END DO
nz = n
IF (ic == n) nz = n - 1
GO TO 40
 
30 nz = kk - 2
 
40 y(1:nz) = czero
RETURN
END SUBROUTINE ckscl
 
 
 
SUBROUTINE cacon(z, fnu, kode, mr, n, y, nz, rl, fnul, tol, elim, alim)
!***BEGIN PROLOGUE  CACON
!***REFER TO  CBESK,CBESH
 
!     CACON APPLIES THE ANALYTIC CONTINUATION FORMULA
 
!         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
!                 MP=PI*MR*CMPLX(0.0,1.0)
 
!     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
!     HALF Z PLANE
 
!***ROUTINES CALLED  CBINU,CBKNU,CS1S2,R1MACH
!***END PROLOGUE  CACON
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: mr
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN OUT)  :: rl
real(dp), INTENT(IN OUT)  :: fnul
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN OUT)  :: elim
real(dp), INTENT(IN OUT)  :: alim
 
COMPLEX (dp)  :: ck, cs, cscl, cscr, csgn, cspn, css(3), csr(3), c1, c2,  &
                 rz, sc1, sc2, st, s1, s2, zn, cy(2)
real(dp)     :: arg, ascle, as2, bscle, bry(3), cpn, c1i, c1m, c1r,  &
                 fmr, sgn, spn, yy
INTEGER       :: i, inu, iuf, kflag, nn, nw
real(dp), PARAMETER     :: pi = 3.14159265358979324_dp
COMPLEX (dp), PARAMETER  :: cone = (1.0_dp, 0.0_dp)
 
nz = 0
zn = -z
nn = n
CALL cbinu(zn, fnu, kode, nn, y, nw, rl, fnul, tol, elim, alim)
IF (nw >= 0) THEN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!-----------------------------------------------------------------------
  nn = min(2, n)
  CALL cbknu(zn, fnu, kode, nn, cy, nw, tol, elim, alim)
  IF (nw == 0) THEN
    s1 = cy(1)
    fmr = mr
    sgn = -sign(pi, fmr)
    csgn = cmplx(0.0_dp, sgn, kind=dp)
    IF (kode /= 1) THEN
      yy = -aimag(zn)
      cpn = cos(yy)
      spn = sin(yy)
      csgn = csgn * cmplx(cpn, spn, kind=dp)
    END IF
!-----------------------------------------------------------------------
!     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
    inu = int(fnu)
    arg = (fnu - inu) * sgn
    cpn = cos(arg)
    spn = sin(arg)
    cspn = cmplx(cpn, spn, kind=dp)
    IF (mod(inu, 2) == 1) cspn = -cspn
    iuf = 0
    c1 = s1
    c2 = y(1)
    ascle = 1.0e+3 * tiny(0.0_dp) / tol
    IF (kode /= 1) THEN
      CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
      nz = nz + nw
      sc1 = c1
    END IF
    y(1) = cspn * c1 + csgn * c2
    IF (n == 1) RETURN
    cspn = -cspn
    s2 = cy(2)
    c1 = s2
    c2 = y(2)
    IF (kode /= 1) THEN
      CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
      nz = nz + nw
      sc2 = c1
    END IF
    y(2) = cspn * c1 + csgn * c2
    IF (n == 2) RETURN
    cspn = -cspn
    rz = 2.0_dp / zn
    ck = cmplx(fnu+1.0_dp, 0.0_dp, kind=dp) * rz
!-----------------------------------------------------------------------
!     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
!-----------------------------------------------------------------------
    cscl = 1.0_dp/tol
    cscr = tol
    css(1) = cscl
    css(2) = cone
    css(3) = cscr
    csr(1) = cscr
    csr(2) = cone
    csr(3) = cscl
    bry(1) = ascle
    bry(2) = 1.0_dp / ascle
    bry(3) = huge(0.0_dp)
    as2 = abs(s2)
    kflag = 2
    IF (as2 <= bry(1)) THEN
      kflag = 1
    ELSE
      IF (as2 >= bry(2)) THEN
        kflag = 3
      END IF
    END IF
    bscle = bry(kflag)
    s1 = s1 * css(kflag)
    s2 = s2 * css(kflag)
    cs = csr(kflag)
    DO  i = 3, n
      st = s2
      s2 = ck * s2 + s1
      s1 = st
      c1 = s2 * cs
      st = c1
      c2 = y(i)
      IF (kode /= 1) THEN
        IF (iuf >= 0) THEN
          CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
          nz = nz + nw
          sc1 = sc2
          sc2 = c1
          IF (iuf == 3) THEN
            iuf = -4
            s1 = sc1 * css(kflag)
            s2 = sc2 * css(kflag)
            st = sc2
          END IF
        END IF
      END IF
      y(i) = cspn * c1 + csgn * c2
      ck = ck + rz
      cspn = -cspn
      IF (kflag < 3) THEN
        c1r = real(c1, kind=dp)
        c1i = aimag(c1)
        c1r = abs(c1r)
        c1i = abs(c1i)
        c1m = max(c1r, c1i)
        IF (c1m > bscle) THEN
          kflag = kflag + 1
          bscle = bry(kflag)
          s1 = s1 * cs
          s2 = st
          s1 = s1 * css(kflag)
          s2 = s2 * css(kflag)
          cs = csr(kflag)
        END IF
      END IF
    END DO
    RETURN
  END IF
END IF
nz = -1
IF (nw == -2) nz = -2
RETURN
END SUBROUTINE cacon
 
 
 
SUBROUTINE cbinu(z, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
!***BEGIN PROLOGUE  CBINU
!***REFER TO  CBESH,CBESI,CBESJ,CBESK,CAIRY,CBIRY
 
!     CBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
 
!***ROUTINES CALLED  CASYI,CBUNI,CMLRI,CSERI,CUOIK,CWRSK
!***END PROLOGUE  CBINU
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: cy(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: rl
real(dp), INTENT(IN)      :: fnul
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: cw(2)
real(dp)     :: az, dfnu
INTEGER       :: inw, nlast, nn, nui, nw
COMPLEX (dp), PARAMETER  :: czero = (0.0_dp, 0.0_dp)
 
nz = 0
az = abs(z)
nn = n
cy = czero
dfnu = fnu + (n-1)
IF (az > 2.0_dp) THEN
  IF (az*az*0.25_dp > dfnu+1.0_dp) GO TO 10
END IF
!-----------------------------------------------------------------------
!     POWER SERIES
!-----------------------------------------------------------------------
CALL cseri(z, fnu, kode, nn, cy, nw, tol, elim, alim)
inw = abs(nw)
nz = nz + inw
nn = nn - inw
IF (nn == 0) RETURN
IF (nw >= 0) GO TO 80
dfnu = fnu + (nn-1)
 
10 IF (az >= rl) THEN
  IF (dfnu > 1.0_dp) THEN
    IF (az+az < dfnu*dfnu) GO TO 20
  END IF
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR LARGE Z
!-----------------------------------------------------------------------
  CALL casyi(z, fnu, kode, nn, cy, nw, rl, tol, elim, alim)
  IF (nw < 0) GO TO 90
  GO TO 80
END IF
IF (dfnu <= 1.0_dp) GO TO 40
!-----------------------------------------------------------------------
!     OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
!-----------------------------------------------------------------------
20 CALL cuoik(z, fnu, kode, 1, nn, cy, nw, tol, elim, alim)
IF (nw < 0) GO TO 90
nz = nz + nw
nn = nn - nw
IF (nn == 0) RETURN
dfnu = fnu + (nn-1)
IF (dfnu > fnul) GO TO 70
IF (az > fnul) GO TO 70
 
30 IF (az > rl) GO TO 50
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE SERIES
!-----------------------------------------------------------------------
40 CALL cmlri(z, fnu, kode, nn, cy, nw, tol)
IF (nw < 0) GO TO 90
GO TO 80
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
!-----------------------------------------------------------------------
!     OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
!-----------------------------------------------------------------------
50 CALL cuoik(z, fnu, kode, 2, 2, cw, nw, tol, elim, alim)
IF (nw < 0) THEN
  nz = nn
  cy(1:nn) = czero
  RETURN
END IF
IF (nw > 0) GO TO 90
CALL cwrsk(z, fnu, kode, nn, cy, nw, cw, tol, elim, alim)
IF (nw < 0) GO TO 90
GO TO 80
!-----------------------------------------------------------------------
!     INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
!-----------------------------------------------------------------------
70 nui = int(fnul - dfnu + 1.0_dp)
nui = max(nui, 0)
CALL cbuni(z, fnu, kode, nn, cy, nw, nui, nlast, fnul, tol, elim, alim)
IF (nw < 0) GO TO 90
nz = nz + nw
IF (nlast /= 0) THEN
  nn = nlast
  GO TO 30
END IF
 
80 RETURN
 
90 nz = -1
IF (nw == -2) nz = -2
RETURN
END SUBROUTINE cbinu
 
 
 
FUNCTION gamln(z) RESULT(fn_val)
 
! N.B. Argument IERR has been removed.
 
!***BEGIN PROLOGUE  GAMLN
!***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
 
!   GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR Z > 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 <= ZMIN.
!   THE FUNCTION WAS MADE AS PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE
!   NUMBER OF BASE 10 DIGITS IN A WORD,
!   RLN = MAX(-LOG10(EPSILON(0.0_dp)), 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      - REAL ARGUMENT, Z > 0.0_dp
 
!      OUTPUT
!        GAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z
!        IERR   - ERROR FLAG
!                 IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
!                 IERR=1, Z <= 0.0_dp,    NO COMPUTATION
 
!***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 BY D. E. AMOS, SAND83-0083, MAY 1983.
!***ROUTINES CALLED  I1MACH,R1MACH
!***END PROLOGUE  GAMLN
 
real(dp), INTENT(IN)  :: z
real(dp)              :: fn_val
 
INTEGER    :: i, i1m, k, mz, nz
real(dp)  :: fln, fz, rln, s, tlg, trm, tst, t1, wdtol,  &
              zdmy, zinc, zm, zmin, zp, zsq
 
!           LNGAMMA(N), N=1,100
real(dp), PARAMETER  :: gln(100) = (/ 0.00000000000000000_dp,  &
  0.00000000000000000_dp, 6.93147180559945309d-01, 1.79175946922805500_dp,   &
  3.17805383034794562_dp, 4.78749174278204599_dp, 6.57925121201010100_dp,    &
  8.52516136106541430_dp, 1.06046029027452502d+01, 1.28018274800814696d+01,  &
  1.51044125730755153d+01, 1.75023078458738858d+01, 1.99872144956618861d+01, &
  2.25521638531234229d+01, 2.51912211827386815d+01, 2.78992713838408916d+01, &
  3.06718601060806728d+01, 3.35050734501368889d+01, 3.63954452080330536d+01, &
  3.93398841871994940d+01, 4.23356164607534850d+01, 4.53801388984769080d+01, &
  4.84711813518352239d+01, 5.16066755677643736d+01, 5.47847293981123192d+01, &
  5.80036052229805199d+01, 6.12617017610020020d+01, 6.45575386270063311d+01, &
  6.78897431371815350d+01, 7.12570389671680090d+01, 7.46582363488301644d+01, &
  7.80922235533153106d+01, 8.15579594561150372d+01, 8.50544670175815174d+01, &
  8.85808275421976788d+01, 9.21361756036870925d+01, 9.57196945421432025d+01, &
  9.93306124547874269d+01, 1.02968198614513813d+02, 1.06631760260643459d+02, &
  1.10320639714757395d+02, 1.14034211781461703d+02, 1.17771881399745072d+02, &
  1.21533081515438634d+02, 1.25317271149356895d+02, 1.29123933639127215d+02, &
  1.32952575035616310d+02, 1.36802722637326368d+02, 1.40673923648234259d+02, &
  1.44565743946344886d+02, 1.48477766951773032d+02, 1.52409592584497358d+02, &
  1.56360836303078785d+02, 1.60331128216630907d+02, 1.64320112263195181d+02, &
  1.68327445448427652d+02, 1.72352797139162802d+02, 1.76395848406997352d+02, &
  1.80456291417543771d+02, 1.84533828861449491d+02, 1.88628173423671591d+02, &
  1.92739047287844902d+02, 1.96866181672889994d+02, 2.01009316399281527d+02, &
  2.05168199482641199d+02, 2.09342586752536836d+02, 2.13532241494563261d+02, &
  2.17736934113954227d+02, 2.21956441819130334d+02, 2.26190548323727593d+02, &
  2.30439043565776952d+02, 2.34701723442818268d+02, 2.38978389561834323d+02, &
  2.43268849002982714d+02, 2.47572914096186884d+02, 2.51890402209723194d+02, &
  2.56221135550009525d+02, 2.60564940971863209d+02, 2.64921649798552801d+02, &
  2.69291097651019823d+02, 2.73673124285693704d+02, 2.78067573440366143d+02, &
  2.82474292687630396d+02, 2.86893133295426994d+02, 2.91323950094270308d+02, &
  2.95766601350760624d+02, 3.00220948647014132d+02, 3.04686856765668715d+02, &
  3.09164193580146922d+02, 3.13652829949879062d+02, 3.18152639620209327d+02, &
  3.22663499126726177d+02, 3.27185287703775217d+02, 3.31717887196928473d+02, &
  3.36261181979198477d+02, 3.40815058870799018d+02, 3.45379407062266854d+02, &
  3.49954118040770237d+02, 3.54539085519440809d+02, 3.59134205369575399d+02 /)
 
!             COEFFICIENTS OF ASYMPTOTIC EXPANSION
real(dp), PARAMETER  :: cf(22) = (/  &
    8.33333333333333333d-02, -2.77777777777777778d-03,  &
    7.93650793650793651d-04, -5.95238095238095238d-04,  &
    8.41750841750841751d-04, -1.91752691752691753d-03,  &
    6.41025641025641026d-03, -2.95506535947712418d-02,  &
    1.79644372368830573d-01, -1.39243221690590112_dp,   &
    1.34028640441683920d+01, -1.56848284626002017d+02,  &
    2.19310333333333333d+03, -3.61087712537249894d+04,  &
    6.91472268851313067d+05, -1.52382215394074162d+07,  &
    3.82900751391414141d+08, -1.08822660357843911d+10,  &
    3.47320283765002252d+11, -1.23696021422692745d+13,  &
    4.88788064793079335d+14, -2.13203339609193739d+16 /)
 
!             LN(2*PI)
real(dp), PARAMETER  :: con = 1.83787706640934548_dp
 
!***FIRST EXECUTABLE STATEMENT  GAMLN
IF (z > 0.0_dp) THEN
  IF (z <= 101.0_dp) THEN
    nz = int(z)
    fz = z - nz
    IF (fz <= 0.0_dp) THEN
      IF (nz <= 100) THEN
        fn_val = gln(nz)
        RETURN
      END IF
    END IF
  END IF
  wdtol = epsilon(0.0_dp)
  wdtol = max(wdtol, 0.5d-18)
  i1m = digits(0.0_dp)
  rln = log10( real( radix(0.0_dp), kind=dp) ) * i1m
  fln = min(rln,20.0_dp)
  fln = max(fln,3.0_dp)
  fln = fln - 3.0_dp
  zm = 1.8000_dp + 0.3875_dp * fln
  mz = int(zm + 1.0_dp)
  zmin = mz
  zdmy = z
  zinc = 0.0_dp
  IF (z < zmin) THEN
    zinc = zmin - nz
    zdmy = z + zinc
  END IF
  zp = 1.0_dp / zdmy
  t1 = cf(1) * zp
  s = t1
  IF (zp >= wdtol) THEN
    zsq = zp * zp
    tst = t1 * wdtol
    DO  k = 2, 22
      zp = zp * zsq
      trm = cf(k) * zp
      IF (abs(trm) < tst) EXIT
      s = s + trm
    END DO
  END IF
 
  IF (zinc == 0.0_dp) THEN
    tlg = log(z)
    fn_val = z * (tlg-1.0_dp) + 0.5_dp * (con-tlg) + s
    RETURN
  END IF
  zp = 1.0_dp
  nz = int(zinc)
  DO  i = 1, nz
    zp = zp * (z + (i-1))
  END DO
  tlg = log(zdmy)
  fn_val = zdmy * (tlg-1.0_dp) - log(zp) + 0.5_dp * (con-tlg) + s
  RETURN
END IF
 
WRITE(*, *) '** ERROR: Zero or -ve argument for function GAMLN **'
fn_val = 0.0
RETURN
END FUNCTION gamln
 
 
 
SUBROUTINE cuchk(y, nz, ascle, tol)
!***BEGIN PROLOGUE  CUCHK
!***REFER TO CSERI,CUOIK,CUNK1,CUNK2,CUNI1,CUNI2,CKSCL
 
!   Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN
!   EXP(-ALIM) = ASCLE = 1.0E+3*TINY(0.0_dp)/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  CUCHK
 
COMPLEX (dp), INTENT(IN)  :: y
INTEGER, INTENT(OUT)      :: nz
real(dp), INTENT(IN)     :: ascle
real(dp), INTENT(IN)     :: tol
 
real(dp)  :: ss, st, yr, yi
 
nz = 0
yr = real(y, kind=dp)
yi = aimag(y)
yr = abs(yr)
yi = abs(yi)
st = min(yr, yi)
IF (st > ascle) RETURN
ss = max(yr, yi)
st = st / tol
IF (ss < st) nz = 1
RETURN
END SUBROUTINE cuchk
 
 
 
SUBROUTINE cacai(z, fnu, kode, mr, n, y, nz, rl, tol, elim, alim)
!***BEGIN PROLOGUE  CACAI
!***REFER TO  CAIRY
 
!  CACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
 
!      K(FNU,ZN*EXP(MP)) = K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
!              MP = PI*MR*CMPLX(0.0,1.0)
 
!  TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT HALF Z PLANE
!  FOR USE WITH CAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
!  CACAI IS THE SAME AS CACON WITH THE PARTS FOR LARGER ORDERS AND
!  RECURRENCE REMOVED.  A RECURSIVE CALL TO CACON CAN RESULT IF CACON
!  IS CALLED FROM CAIRY.
 
!***ROUTINES CALLED  CASYI,CBKNU,CMLRI,CSERI,CS1S2,R1MACH
!***END PROLOGUE  CACAI
 
COMPLEX (dp), INTENT(IN)   :: z
real(dp), INTENT(IN)      :: fnu
INTEGER, INTENT(IN)        :: kode
INTEGER, INTENT(IN OUT)    :: mr
INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(OUT)  :: y(n)
INTEGER, INTENT(OUT)       :: nz
real(dp), INTENT(IN)      :: rl
real(dp), INTENT(IN)      :: tol
real(dp), INTENT(IN)      :: elim
real(dp), INTENT(IN)      :: alim
 
COMPLEX (dp)  :: csgn, cspn, c1, c2, zn, cy(2)
real(dp)     :: arg, ascle, az, cpn, dfnu, fmr, sgn, spn, yy
INTEGER       :: inu, iuf, nn, nw
real(dp), PARAMETER  :: pi = 3.14159265358979324_dp
 
nz = 0
zn = -z
az = abs(z)
nn = n
dfnu = fnu + (n-1)
IF (az > 2.0_dp) THEN
  IF (az*az*0.25_dp > dfnu+1.0_dp) GO TO 10
END IF
!-----------------------------------------------------------------------
!     POWER SERIES FOR THE I FUNCTION
!-----------------------------------------------------------------------
CALL cseri(zn, fnu, kode, nn, y, nw, tol, elim, alim)
GO TO 20
 
10 IF (az >= rl) THEN
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
!-----------------------------------------------------------------------
  CALL casyi(zn, fnu, kode, nn, y, nw, rl, tol, elim, alim)
  IF (nw < 0) GO TO 30
ELSE
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
!-----------------------------------------------------------------------
  CALL cmlri(zn, fnu, kode, nn, y, nw, tol)
  IF (nw < 0) GO TO 30
END IF
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!-----------------------------------------------------------------------
20 CALL cbknu(zn, fnu, kode, 1, cy, nw, tol, elim, alim)
IF (nw == 0) THEN
  fmr = mr
  sgn = -sign(pi, fmr)
  csgn = cmplx(0.0_dp, sgn, kind=dp)
  IF (kode /= 1) THEN
    yy = -aimag(zn)
    cpn = cos(yy)
    spn = sin(yy)
    csgn = csgn * cmplx(cpn, spn, kind=dp)
  END IF
!-----------------------------------------------------------------------
!     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
  inu = int(fnu)
  arg = (fnu - inu) * sgn
  cpn = cos(arg)
  spn = sin(arg)
  cspn = cmplx(cpn, spn, kind=dp)
  IF (mod(inu,2) == 1) cspn = -cspn
  c1 = cy(1)
  c2 = y(1)
  IF (kode /= 1) THEN
    iuf = 0
    ascle = 1.0e+3 * tiny(0.0_dp) / tol
    CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
    nz = nz + nw
  END IF
  y(1) = cspn * c1 + csgn * c2
  RETURN
END IF
 
30 nz = -1
IF (nw == -2) nz = -2
RETURN
END SUBROUTINE cacai
 
END Module complex_bessel