bessely.f90

module bessely
    implicit none
    private
    public zbesy
contains
!*****************************************************************
!* EVALUATE A Y-BESSEL FUNCTION OF COMPLEX ARGUMENT (SECOND KIND)*
!* ------------------------------------------------------------- *
!* SAMPLE RUN:                                                   *
!* (Evaluate Y0 to Y4 for argument Z=(1.0,2.0) ).                *
!*                                                               *
!* zr(0) =   1.367419                                            *
!* zi(0) =   1.521507                                            *
!* zr(1) =  -1.089470                                            *
!* zi(1) =   1.314951                                            *
!* zr(2) =  -0.751245                                            *
!* zi(2) =  -0.123950                                            *
!* zr(3) =   0.290153                                            *
!* zi(3) =  -0.212119                                            *
!* zr(4) =   0.590344                                            *
!* zi(4) =  -0.826960                                            *
!* NZ =          0                                               *
!* Error code:           0                                       *
!*                                                               *
!* ------------------------------------------------------------- *
!* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77    *
!*       [BIBLI 18].                                             *
!*                                                               *
!*                        F90 Release 1.0 By J-P Moreau, Paris   *
!*                                 (www.jpmoreau.fr)             *
!*****************************************************************
!PROGRAM TEST_ZBESY

!  real*8 zr, zi
!  real*8 cyr(10), cyi(10), cwr(10), cwi(10)

!  n=5
!  zr=1.d0; zi=2.d0

!  call ZBESY(zr,zi,0.d0,1,n,cyr,cyi,nz,cwr,cwi,ierr)

!  print *,' '
!  do i=1, n
!    write(*,10) i-1, cyr(i)
!    write(*,11) i-1, cyi(i)
!  end do
!  print *,' NZ=', NZ
!  print *,' Error code:', ierr
!  print *,' '
!  stop

!10 format('  zr(',I1,') = ',F10.6)
!11 format('  zi(',I1,') = ',F10.6)

!END


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


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

      DATA hpi /1.57079632679489662d0/

!***FIRST EXECUTABLE STATEMENT  ZBESH
      ierr = 0
      nz=0
      IF (zr.EQ.0.0d0 .AND. zi.EQ.0.0d0) ierr=1
      IF (fnu.LT.0.0d0) ierr=1
      IF (m.LT.1 .OR. m.GT.2) ierr=1
      IF (kode.LT.1 .OR. kode.GT.2) ierr=1
      IF (n.LT.1) ierr=1
      IF (ierr.NE.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).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
!     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
!     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
!     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
!     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
!     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
!-----------------------------------------------------------------------
      tol = dmax1(d1mach(4),1.0d-18)
      k1 = i1mach(15)
      k2 = i1mach(16)
      r1m5 = d1mach(5)
      k = min0(iabs(k1),iabs(k2))
      elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
      k1 = i1mach(14) - 1
      aa = r1m5*dble(float(k1))
      dig = dmin1(aa,18.0d0)
      aa = aa*2.303d0
      alim = elim + dmax1(-aa,-41.45d0)
      fnul = 10.0d0 + 6.0d0*(dig-3.0d0)
      rl = 1.2d0*dig + 3.0d0
      fn = fnu + dble(float(nn-1))
      mm = 3 - m - m
      fmm = dble(float(mm))
      znr = fmm*zi
      zni = -fmm*zr
!-----------------------------------------------------------------------
!     TEST FOR PROPER RANGE
!-----------------------------------------------------------------------
      az = zabs(zr,zi)
      aa = 0.5d0/tol
      bb=dble(float(i1mach(9)))*0.5d0
      aa = dmin1(aa,bb)
      IF (az.GT.aa) GO TO 260
      IF (fn.GT.aa) GO TO 260
      aa = dsqrt(aa)
      IF (az.GT.aa) ierr=3
      IF (fn.GT.aa) ierr=3
!-----------------------------------------------------------------------
!     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
!-----------------------------------------------------------------------
      ufl = dexp(-elim)
      IF (az.LT.ufl) GO TO 230
      IF (fnu.GT.fnul) GO TO 90
      IF (fn.LE.1.0d0) GO TO 70
      IF (fn.GT.2.0d0) GO TO 60
      IF (az.GT.tol) GO TO 70
      arg = 0.5d0*az
      aln = -fn*dlog(arg)
      IF (aln.GT.elim) GO TO 230
      GO TO 70
   60 CONTINUE
      CALL zuoik(znr, zni, fnu, kode, 2, nn, cyr, cyi, nuf, tol, elim, alim)
      IF (nuf.LT.0) GO TO 230
      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.EQ.0) GO TO 140
   70 CONTINUE
      IF ((znr.LT.0.0d0) .OR. (znr.EQ.0.0d0 .AND. zni.LT.0.0d0 .AND. &
      m.EQ.2)) GO TO 80
!-----------------------------------------------------------------------
!     RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
!     YN.GE.0. .OR. M=1)
!-----------------------------------------------------------------------
      CALL zbknu(znr, zni, fnu, kode, nn, cyr, cyi, nz, tol, elim, alim)
      GO TO 110
!-----------------------------------------------------------------------
!     LEFT HALF PLANE COMPUTATION
!-----------------------------------------------------------------------
   80 CONTINUE
      mr = -mm
      CALL zacon(znr, zni, fnu, kode, mr, nn, cyr, cyi, nw, rl, fnul, &
      tol, elim, alim)
      IF (nw.LT.0) GO TO 240
      nz=nw
      GO TO 110
   90 CONTINUE
!-----------------------------------------------------------------------
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
!-----------------------------------------------------------------------
      mr = 0
      IF ((znr.GE.0.0d0) .AND. (znr.NE.0.0d0 .OR. zni.GE.0.0d0 .OR. &
      m.NE.2)) GO TO 100
      mr = -mm
      IF (znr.NE.0.0d0 .OR. zni.GE.0.0d0) GO TO 100
      znr = -znr
      zni = -zni
  100 CONTINUE
      CALL zbunk(znr, zni, fnu, kode, mr, nn, cyr, cyi, nw, tol, elim, alim)
      IF (nw.LT.0) GO TO 240
      nz = nz + nw
  110 CONTINUE
!-----------------------------------------------------------------------
!     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
!-----------------------------------------------------------------------
      sgn = dsign(hpi,-fmm)
!-----------------------------------------------------------------------
!     CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
      inu = int(sngl(fnu))
      inuh = inu/2
      ir = inu - 2*inuh
      arg = (fnu-dble(float(inu-ir)))*sgn
      rhpi = 1.0d0/sgn
      zni = rhpi*dcos(arg)
      znr = -rhpi*dsin(arg)
      IF (mod(inuh,2).EQ.0) GO TO 120
      znr = -znr
      zni = -zni
  120 CONTINUE
      zti = -fmm
      DO 130 i=1,nn
        str = cyr(i)*znr - cyi(i)*zni
        cyi(i) = cyr(i)*zni + cyi(i)*znr
        cyr(i) = str
        str = -zni*zti
        zni = znr*zti
        znr = str
  130 CONTINUE
      RETURN
  140 CONTINUE
      IF (znr.LT.0.0d0) GO TO 230
      RETURN
  230 CONTINUE
      nz=0
      ierr=2
      RETURN
  240 CONTINUE
      IF(nw.EQ.(-1)) GO TO 230
      nz=0
      ierr=5
      RETURN
  260 CONTINUE
      nz=0
      ierr=4
      RETURN
END

SUBROUTINE zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZUOIK
!***REFER TO  ZBESI,ZBESK,ZBESH
!
!     ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
!     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
!     (IN LOGARITHMI! FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
!     WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
!     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
!     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
!     MULTIPLIERS (IN LOGARITHMI! FORM) IS MADE BASED ON ELIM. HERE
!     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
!     EXP(-ELIM)/TOL
!
!     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
!          =2 MEANS THE K SEQUENCE IS TESTED
!     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
!         =-1 MEANS AN OVERFLOW WOULD OCCUR
!     IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
!             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
!     IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
!     IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
!             ANOTHER ROUTINE
!
!***ROUTINES CALLED  ZUCHK,ZUNHJ,ZUNIK,D1MACH,ZABS,ZLOG
!***END PROLOGUE  ZUOIK
!     COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
!    *ZR
      DOUBLE PRECISION aarg, aic, alim, aphi, argi, argr, asumi, asumr, &
      ascle, ax, ay, bsumi, bsumr, cwrki, cwrkr, czi, czr, elim, fnn,  &
      fnu, gnn, gnu, phii, phir, rcz, str, sti, sumi, sumr, tol, yi,   &
      yr, zbi, zbr, zeroi, zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi,  &
      zni, znr, zr, zri, zrr
      INTEGER i, idum, iform, ikflg, init, kode, n, nn, nuf, nw
      dimension yr(1), yi(1), cwrkr(16), cwrki(16)
      DATA zeror,zeroi / 0.0d0, 0.0d0 /
      DATA aic / 1.265512123484645396d+00 /
      nuf = 0
      nn = n
      zrr = zr
      zri = zi
      IF (zr.GE.0.0d0) GO TO 10
      zrr = -zr
      zri = -zi
   10 CONTINUE
      zbr = zrr
      zbi = zri
      ax = dabs(zr)*1.7321d0
      ay = dabs(zi)
      iform = 1
      IF (ay.GT.ax) iform = 2
      gnu = dmax1(fnu,1.0d0)
      IF (ikflg.EQ.1) GO TO 20
      fnn = dble(float(nn))
      gnn = fnu + fnn - 1.0d0
      gnu = dmax1(gnn,fnn)
   20 CONTINUE
!-----------------------------------------------------------------------
!     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.EQ.2) GO TO 30
      init = 0
      CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii, &
      zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      GO TO 50
   30 CONTINUE
      znr = zri
      zni = -zrr
      IF (zi.GT.0.0d0) GO TO 40
      znr = -znr
   40 CONTINUE
      CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r, &
      zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      aarg = zabs(argr,argi)
   50 CONTINUE
      IF (kode.EQ.1) GO TO 60
      czr = czr - zbr
      czi = czi - zbi
   60 CONTINUE
      IF (ikflg.EQ.1) GO TO 70
      czr = -czr
      czi = -czi
   70 CONTINUE
      aphi = zabs(phir,phii)
      rcz = czr
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
      IF (rcz.GT.elim) GO TO 210
      IF (rcz.LT.alim) GO TO 80
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.elim) GO TO 210
      GO TO 130
   80 CONTINUE
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
      IF (rcz.LT.(-elim)) GO TO 90
      IF (rcz.GT.(-alim)) GO TO 130
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.(-elim)) GO TO 110
   90 CONTINUE
      DO 100 i=1,nn
        yr(i) = zeror
        yi(i) = zeroi
  100 CONTINUE
      nuf = nn
      RETURN
  110 CONTINUE
      ascle = 1.0d+3*d1mach(1)/tol
      CALL zlog(phir, phii, str, sti, idum)
      czr = czr + str
      czi = czi + sti
      IF (iform.EQ.1) GO TO 120
      CALL zlog(argr, argi, str, sti, idum)
      czr = czr - 0.25d0*str - aic
      czi = czi - 0.25d0*sti
  120 CONTINUE
      ax = dexp(rcz)/tol
      ay = czi
      czr = ax*dcos(ay)
      czi = ax*dsin(ay)
      CALL zuchk(czr, czi, nw, ascle, tol)
      IF (nw.NE.0) GO TO 90
  130 CONTINUE
      IF (ikflg.EQ.2) RETURN
      IF (n.EQ.1) RETURN
!-----------------------------------------------------------------------
!     SET UNDERFLOWS ON I SEQUENCE
!-----------------------------------------------------------------------
  140 CONTINUE
      gnu = fnu + dble(float(nn-1))
      IF (iform.EQ.2) GO TO 150
      init = 0
      CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii, &
      zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      GO TO 160
  150 CONTINUE
      CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r, &
      zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      czr = -zeta1r + zeta2r
      czi = -zeta1i + zeta2i
      aarg = zabs(argr,argi)
  160 CONTINUE
      IF (kode.EQ.1) GO TO 170
      czr = czr - zbr
      czi = czi - zbi
  170 CONTINUE
      aphi = zabs(phir,phii)
      rcz = czr
      IF (rcz.LT.(-elim)) GO TO 180
      IF (rcz.GT.(-alim)) RETURN
      rcz = rcz + dlog(aphi)
      IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
      IF (rcz.GT.(-elim)) GO TO 190
  180 CONTINUE
      yr(nn) = zeror
      yi(nn) = zeroi
      nn = nn - 1
      nuf = nuf + 1
      IF (nn.EQ.0) RETURN
      GO TO 140
  190 CONTINUE
      ascle = 1.0d+3*d1mach(1)/tol
      CALL zlog(phir, phii, str, sti, idum)
      czr = czr + str
      czi = czi + sti
      IF (iform.EQ.1) GO TO 200
      CALL zlog(argr, argi, str, sti, idum)
      czr = czr - 0.25d0*str - aic
      czi = czi - 0.25d0*sti
  200 CONTINUE
      ax = dexp(rcz)/tol
      ay = czi
      czr = ax*dcos(ay)
      czi = ax*dsin(ay)
      CALL zuchk(czr, czi, nw, ascle, tol)
      IF (nw.NE.0) GO TO 180
      RETURN
  210 CONTINUE
      nuf = -1
      RETURN
      END

SUBROUTINE zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZBKNU
!***REFER TO  ZBESI,ZBESK,ZAIRY,ZBESH
!
!     ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
!
!***ROUTINES CALLED  DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,ZABS,ZDIV,
!                    ZEXP,ZLOG,ZMLT,ZSQRT
!***END PROLOGUE  ZBKNU
!
      DOUBLE PRECISION aa, ak, alim, ascle, a1, a2, bb, bk, bry, caz, &
      cbi, cbr, cc, cchi, cchr, cki, ckr, coefi, coefr, conei, coner, &
      crscr, csclr, cshi, cshr, csi, csr, csrr, cssr, ctwoi, ctwor,   &
      czeroi, czeror, czi, czr, dnu, dnu2, dpi, elim, etest, fc, fhs, &
      fi, fk, fks, fmui, fmur, fnu, fpi, fr, g1, g2, hpi, pi, pr, pti,&
      ptr, p1i, p1r, p2i, p2m, p2r, qi, qr, rak, rcaz, rthpi, rzi,    &
      rzr, r1, s, smui, smur, spi, sti, str, s1i, s1r, s2i, s2r, tm,  &
      tol, tth, t1, t2, yi, yr, zi, zr, elm, celmr, zdr, zdi, &
      as, alas, helim, cyr, cyi!, DGAMLN
      INTEGER i, iflag, inu, k, kflag, kk, kmax, kode, koded, n, nz,  &
      idum, j, ic, inub, nw
      dimension yr(n), yi(n), cc(8), cssr(3), csrr(3), bry(3), cyr(2),&
      cyi(2)
!     COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
!     COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK

      DATA kmax / 30 /
      DATA czeror,czeroi,coner,conei,ctwor,ctwoi,r1/ &
        0.0d0 , 0.0d0 , 1.0d0 , 0.0d0 , 2.0d0 , 0.0d0 , 2.0d0 /
      DATA dpi, rthpi, spi ,hpi, fpi, tth /  &
           3.14159265358979324d0,       1.25331413731550025d0, &
           1.90985931710274403d0,       1.57079632679489662d0, &
           1.89769999331517738d0,       6.66666666666666666d-01/
      DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), 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/

      caz = zabs(zr,zi)
      csclr = 1.0d0/tol
      crscr = tol
      cssr(1) = csclr
      cssr(2) = 1.0d0
      cssr(3) = crscr
      csrr(1) = crscr
      csrr(2) = 1.0d0
      csrr(3) = csclr
      bry(1) = 1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      nz = 0
      iflag = 0
      koded = kode
      rcaz = 1.0d0/caz
      str = zr*rcaz
      sti = -zi*rcaz
      rzr = (str+str)*rcaz
      rzi = (sti+sti)*rcaz
      inu = int(sngl(fnu+0.5d0))
      dnu = fnu - dble(float(inu))
      IF (dabs(dnu).EQ.0.5d0) GO TO 110
      dnu2 = 0.0d0
      IF (dabs(dnu).GT.tol) dnu2 = dnu*dnu
      IF (caz.GT.r1) GO TO 110
!-----------------------------------------------------------------------
!     SERIES FOR CABS(Z).LE.R1
!-----------------------------------------------------------------------
      fc = 1.0d0
      CALL zlog(rzr, rzi, smur, smui, idum)
      fmur = smur*dnu
      fmui = smui*dnu
      CALL zshch(fmur, fmui, cshr, cshi, cchr, cchi)
      IF (dnu.EQ.0.0d0) GO TO 10
      fc = dnu*dpi
      fc = fc/dsin(fc)
      smur = cshr/dnu
      smui = cshi/dnu
   10 CONTINUE
      a2 = 1.0d0 + dnu
!-----------------------------------------------------------------------
!     GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
!-----------------------------------------------------------------------
      t2 = dexp(-dgamln(a2,idum))
      t1 = 1.0d0/(t2*fc)
      IF (dabs(dnu).GT.0.1d0) GO TO 40
!-----------------------------------------------------------------------
!     SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
!-----------------------------------------------------------------------
      ak = 1.0d0
      s = cc(1)
      DO 20 k=2,8
        ak = ak*dnu2
        tm = cc(k)*ak
        s = s + tm
        IF (dabs(tm).LT.tol) GO TO 30
   20 CONTINUE
   30 g1 = -s
      GO TO 50
   40 CONTINUE
      g1 = (t1-t2)/(dnu+dnu)
   50 CONTINUE
      g2 = (t1+t2)*0.5d0
      fr = fc*(cchr*g1+smur*g2)
      fi = fc*(cchi*g1+smui*g2)
      CALL zexp(fmur, fmui, str, sti)
      pr = 0.5d0*str/t2
      pi = 0.5d0*sti/t2
      CALL zdiv(0.5d0, 0.0d0, str, sti, ptr, pti)
      qr = ptr/t1
      qi = pti/t1
      s1r = fr
      s1i = fi
      s2r = pr
      s2i = pi
      ak = 1.0d0
      a1 = 1.0d0
      ckr = coner
      cki = conei
      bk = 1.0d0 - dnu2
      IF (inu.GT.0 .OR. n.GT.1) GO TO 80
!-----------------------------------------------------------------------
!     GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
!-----------------------------------------------------------------------
      IF (caz.LT.tol) GO TO 70
      CALL zmlt(zr, zi, zr, zi, czr, czi)
      czr = 0.25d0*czr
      czi = 0.25d0*czi
      t1 = 0.25d0*caz*caz
   60 CONTINUE
      fr = (fr*ak+pr+qr)/bk
      fi = (fi*ak+pi+qi)/bk
      str = 1.0d0/(ak-dnu)
      pr = pr*str
      pi = pi*str
      str = 1.0d0/(ak+dnu)
      qr = qr*str
      qi = qi*str
      str = ckr*czr - cki*czi
      rak = 1.0d0/ak
      cki = (ckr*czi+cki*czr)*rak
      ckr = str*rak
      s1r = ckr*fr - cki*fi + s1r
      s1i = ckr*fi + cki*fr + s1i
      a1 = a1*t1*rak
      bk = bk + ak + ak + 1.0d0
      ak = ak + 1.0d0
      IF (a1.GT.tol) GO TO 60
   70 CONTINUE
      yr(1) = s1r
      yi(1) = s1i
      IF (koded.EQ.1) RETURN
      CALL zexp(zr, zi, str, sti)
      CALL zmlt(s1r, s1i, str, sti, yr(1), yi(1))
      RETURN
!-----------------------------------------------------------------------
!     GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
!-----------------------------------------------------------------------
   80 CONTINUE
      IF (caz.LT.tol) GO TO 100
      CALL zmlt(zr, zi, zr, zi, czr, czi)
      czr = 0.25d0*czr
      czi = 0.25d0*czi
      t1 = 0.25d0*caz*caz
   90 CONTINUE
      fr = (fr*ak+pr+qr)/bk
      fi = (fi*ak+pi+qi)/bk
      str = 1.0d0/(ak-dnu)
      pr = pr*str
      pi = pi*str
      str = 1.0d0/(ak+dnu)
      qr = qr*str
      qi = qi*str
      str = ckr*czr - cki*czi
      rak = 1.0d0/ak
      cki = (ckr*czi+cki*czr)*rak
      ckr = str*rak
      s1r = ckr*fr - cki*fi + s1r
      s1i = ckr*fi + cki*fr + s1i
      str = pr - fr*ak
      sti = pi - fi*ak
      s2r = ckr*str - cki*sti + s2r
      s2i = ckr*sti + cki*str + s2i
      a1 = a1*t1*rak
      bk = bk + ak + ak + 1.0d0
      ak = ak + 1.0d0
      IF (a1.GT.tol) GO TO 90
  100 CONTINUE
      kflag = 2
      a1 = fnu + 1.0d0
      ak = a1*dabs(smur)
      IF (ak.GT.alim) kflag = 3
      str = cssr(kflag)
      p2r = s2r*str
      p2i = s2i*str
      CALL zmlt(p2r, p2i, rzr, rzi, s2r, s2i)
      s1r = s1r*str
      s1i = s1i*str
      IF (koded.EQ.1) GO TO 210
      CALL zexp(zr, zi, fr, fi)
      CALL zmlt(s1r, s1i, fr, fi, s1r, s1i)
      CALL zmlt(s2r, s2i, fr, fi, s2r, s2i)
      GO TO 210
!-----------------------------------------------------------------------
!     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
!-----------------------------------------------------------------------
  110 CONTINUE
      CALL zsqrt(zr, zi, str, sti)
      CALL zdiv(rthpi, czeroi, str, sti, coefr, coefi)
      kflag = 2
      IF (koded.EQ.2) GO TO 120
      IF (zr.GT.alim) GO TO 290
!     BLANK LINE
      str = dexp(-zr)*cssr(kflag)
      sti = -str*dsin(zi)
      str = str*dcos(zi)
      CALL zmlt(coefr, coefi, str, sti, coefr, coefi)
  120 CONTINUE
      IF (dabs(dnu).EQ.0.5d0) GO TO 300
!-----------------------------------------------------------------------
!     MILLER ALGORITHM FOR CABS(Z).GT.R1
!-----------------------------------------------------------------------
      ak = dcos(dpi*dnu)
      ak = dabs(ak)
      IF (ak.EQ.czeror) GO TO 300
      fhs = dabs(0.25d0-dnu2)
      IF (fhs.EQ.czeror) GO TO 300
!-----------------------------------------------------------------------
!     COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
!     DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
!     12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
!     TOL WHERE B IS THE BASE OF THE ARITHMETIC.
!-----------------------------------------------------------------------
      t1 = dble(float(i1mach(14)-1))
      t1 = t1*d1mach(5)*3.321928094d0
      t1 = dmax1(t1,12.0d0)
      t1 = dmin1(t1,60.0d0)
      t2 = tth*t1 - 6.0d0
      IF (zr.NE.0.0d0) GO TO 130
      t1 = hpi
      GO TO 140
  130 CONTINUE
      t1 = datan(zi/zr)
      t1 = dabs(t1)
  140 CONTINUE
      IF (t2.GT.caz) GO TO 170
!-----------------------------------------------------------------------
!     FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
!-----------------------------------------------------------------------
      etest = ak/(dpi*caz*tol)
      fk = coner
      IF (etest.LT.coner) GO TO 180
      fks = ctwor
      ckr = caz + caz + ctwor
      p1r = czeror
      p2r = coner
      DO 150 i=1,kmax
        ak = fhs/fks
        cbr = ckr/(fk+coner)
        ptr = p2r
        p2r = cbr*p2r - p1r*ak
        p1r = ptr
        ckr = ckr + ctwor
        fks = fks + fk + fk + ctwor
        fhs = fhs + fk + fk
        fk = fk + coner
        str = dabs(p2r)*fk
        IF (etest.LT.str) GO TO 160
  150 CONTINUE
      GO TO 310
  160 CONTINUE
      fk = fk + spi*t1*dsqrt(t2/caz)
      fhs = dabs(0.25d0-dnu2)
      GO TO 180
  170 CONTINUE
!-----------------------------------------------------------------------
!     COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
!-----------------------------------------------------------------------
      a2 = dsqrt(caz)
      ak = fpi*ak/(tol*dsqrt(a2))
      aa = 3.0d0*t1/(1.0d0+caz)
      bb = 14.7d0*t1/(28.0d0+caz)
      ak = (dlog(ak)+caz*dcos(aa)/(1.0d0+0.008d0*caz))/dcos(bb)
      fk = 0.12125d0*ak*ak/caz + 1.5d0
  180 CONTINUE
!-----------------------------------------------------------------------
!     BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
!-----------------------------------------------------------------------
      k = int(sngl(fk))
      fk = dble(float(k))
      fks = fk*fk
      p1r = czeror
      p1i = czeroi
      p2r = tol
      p2i = czeroi
      csr = p2r
      csi = p2i
      DO 190 i=1,k
        a1 = fks - fk
        ak = (fks+fk)/(a1+fhs)
        rak = 2.0d0/(fk+coner)
        cbr = (fk+zr)*rak
        cbi = zi*rak
        ptr = p2r
        pti = p2i
        p2r = (ptr*cbr-pti*cbi-p1r)*ak
        p2i = (pti*cbr+ptr*cbi-p1i)*ak
        p1r = ptr
        p1i = pti
        csr = csr + p2r
        csi = csi + p2i
        fks = a1 - fk + coner
        fk = fk - coner
  190 CONTINUE
!-----------------------------------------------------------------------
!     COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
!     SCALING
!-----------------------------------------------------------------------
      tm = zabs(csr,csi)
      ptr = 1.0d0/tm
      s1r = p2r*ptr
      s1i = p2i*ptr
      csr = csr*ptr
      csi = -csi*ptr
      CALL zmlt(coefr, coefi, s1r, s1i, str, sti)
      CALL zmlt(str, sti, csr, csi, s1r, s1i)
      IF (inu.GT.0 .OR. n.GT.1) GO TO 200
      zdr = zr
      zdi = zi
      IF(iflag.EQ.1) GO TO 270
      GO TO 240
  200 CONTINUE
!-----------------------------------------------------------------------
!     COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
!-----------------------------------------------------------------------
      tm = zabs(p2r,p2i)
      ptr = 1.0d0/tm
      p1r = p1r*ptr
      p1i = p1i*ptr
      p2r = p2r*ptr
      p2i = -p2i*ptr
      CALL zmlt(p1r, p1i, p2r, p2i, ptr, pti)
      str = dnu + 0.5d0 - ptr
      sti = -pti
      CALL zdiv(str, sti, zr, zi, str, sti)
      str = str + 1.0d0
      CALL zmlt(str, sti, s1r, s1i, s2r, s2i)
!-----------------------------------------------------------------------
!     FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
!     SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
!-----------------------------------------------------------------------
  210 CONTINUE
      str = dnu + 1.0d0
      ckr = str*rzr
      cki = str*rzi
      IF (n.EQ.1) inu = inu - 1
      IF (inu.GT.0) GO TO 220
      IF (n.GT.1) GO TO 215
      s1r = s2r
      s1i = s2i
  215 CONTINUE
      zdr = zr
      zdi = zi
      IF(iflag.EQ.1) GO TO 270
      GO TO 240
  220 CONTINUE
      inub = 1
      IF(iflag.EQ.1) GO TO 261
  225 CONTINUE
      p1r = csrr(kflag)
      ascle = bry(kflag)
      DO 230 i=inub,inu
        str = s2r
        sti = s2i
        s2r = ckr*str - cki*sti + s1r
        s2i = ckr*sti + cki*str + s1i
        s1r = str
        s1i = sti
        ckr = ckr + rzr
        cki = cki + rzi
        IF (kflag.GE.3) GO TO 230
        p2r = s2r*p1r
        p2i = s2i*p1r
        str = dabs(p2r)
        sti = dabs(p2i)
        p2m = dmax1(str,sti)
        IF (p2m.LE.ascle) GO TO 230
        kflag = kflag + 1
        ascle = bry(kflag)
        s1r = s1r*p1r
        s1i = s1i*p1r
        s2r = p2r
        s2i = p2i
        str = cssr(kflag)
        s1r = s1r*str
        s1i = s1i*str
        s2r = s2r*str
        s2i = s2i*str
        p1r = csrr(kflag)
  230 CONTINUE
      IF (n.NE.1) GO TO 240
      s1r = s2r
      s1i = s2i
  240 CONTINUE
      str = csrr(kflag)
      yr(1) = s1r*str
      yi(1) = s1i*str
      IF (n.EQ.1) RETURN
      yr(2) = s2r*str
      yi(2) = s2i*str
      IF (n.EQ.2) RETURN
      kk = 2
  250 CONTINUE
      kk = kk + 1
      IF (kk.GT.n) RETURN
      p1r = csrr(kflag)
      ascle = bry(kflag)
      DO 260 i=kk,n
        p2r = s2r
        p2i = s2i
        s2r = ckr*p2r - cki*p2i + s1r
        s2i = cki*p2r + ckr*p2i + s1i
        s1r = p2r
        s1i = p2i
        ckr = ckr + rzr
        cki = cki + rzi
        p2r = s2r*p1r
        p2i = s2i*p1r
        yr(i) = p2r
        yi(i) = p2i
        IF (kflag.GE.3) GO TO 260
        str = dabs(p2r)
        sti = dabs(p2i)
        p2m = dmax1(str,sti)
        IF (p2m.LE.ascle) GO TO 260
        kflag = kflag + 1
        ascle = bry(kflag)
        s1r = s1r*p1r
        s1i = s1i*p1r
        s2r = p2r
        s2i = p2i
        str = cssr(kflag)
        s1r = s1r*str
        s1i = s1i*str
        s2r = s2r*str
        s2i = s2i*str
        p1r = csrr(kflag)
  260 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
!-----------------------------------------------------------------------
  261 CONTINUE
      helim = 0.5d0*elim
      elm = dexp(-elim)
      celmr = elm
      ascle = bry(1)
      zdr = zr
      zdi = zi
      ic = -1
      j = 2
      DO 262 i=1,inu
        str = s2r
        sti = s2i
        s2r = str*ckr-sti*cki+s1r
        s2i = sti*ckr+str*cki+s1i
        s1r = str
        s1i = sti
        ckr = ckr+rzr
        cki = cki+rzi
        as = zabs(s2r,s2i)
        alas = dlog(as)
        p2r = -zdr+alas
        IF(p2r.LT.(-elim)) GO TO 263
        CALL zlog(s2r,s2i,str,sti,idum)
        p2r = -zdr+str
        p2i = -zdi+sti
        p2m = dexp(p2r)/tol
        p1r = p2m*dcos(p2i)
        p1i = p2m*dsin(p2i)
        CALL zuchk(p1r,p1i,nw,ascle,tol)
        IF(nw.NE.0) GO TO 263
        j = 3 - j
        cyr(j) = p1r
        cyi(j) = p1i
        IF(ic.EQ.(i-1)) GO TO 264
        ic = i
        GO TO 262
  263   CONTINUE
        IF(alas.LT.helim) GO TO 262
        zdr = zdr-elim
        s1r = s1r*celmr
        s1i = s1i*celmr
        s2r = s2r*celmr
        s2i = s2i*celmr
  262 CONTINUE
      IF(n.NE.1) GO TO 270
      s1r = s2r
      s1i = s2i
      GO TO 270
  264 CONTINUE
      kflag = 1
      inub = i+1
      s2r = cyr(j)
      s2i = cyi(j)
      j = 3 - j
      s1r = cyr(j)
      s1i = cyi(j)
      IF(inub.LE.inu) GO TO 225
      IF(n.NE.1) GO TO 240
      s1r = s2r
      s1i = s2i
      GO TO 240
  270 CONTINUE
      yr(1) = s1r
      yi(1) = s1i
      IF(n.EQ.1) GO TO 280
      yr(2) = s2r
      yi(2) = s2i
  280 CONTINUE
      ascle = bry(1)
      CALL zkscl(zdr,zdi,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim)
      inu = n - nz
      IF (inu.LE.0) RETURN
      kk = nz + 1
      s1r = yr(kk)
      s1i = yi(kk)
      yr(kk) = s1r*csrr(1)
      yi(kk) = s1i*csrr(1)
      IF (inu.EQ.1) RETURN
      kk = nz + 2
      s2r = yr(kk)
      s2i = yi(kk)
      yr(kk) = s2r*csrr(1)
      yi(kk) = s2i*csrr(1)
      IF (inu.EQ.2) RETURN
      t2 = fnu + dble(float(kk-1))
      ckr = t2*rzr
      cki = t2*rzi
      kflag = 1
      GO TO 250
  290 CONTINUE
!-----------------------------------------------------------------------
!     SCALE BY DEXP(Z), IFLAG = 1 CASES
!-----------------------------------------------------------------------
      koded = 2
      iflag = 1
      kflag = 2
      GO TO 120
!-----------------------------------------------------------------------
!     FNU=HALF ODD INTEGER CASE, DNU=-0.5
!-----------------------------------------------------------------------
  300 CONTINUE
      s1r = coefr
      s1i = coefi
      s2r = coefr
      s2i = coefi
      GO TO 210

  310 CONTINUE
      nz=-2
      RETURN
END

SUBROUTINE zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZACON
!***REFER TO  ZBESK,ZBESH
!
!     ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
!
!         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
!                 MP=PI*MR*CMPLX(0.0,1.0)
!
!     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
!     HALF Z PLANE
!
!***ROUTINES CALLED  ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT
!***END PROLOGUE  ZACON
!     COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
!    *S1,S2,Y,Z,ZN
      DOUBLE PRECISION alim, arg, ascle, as2, azn, bry, bscle, cki,    &
      ckr, conei, coner, cpn, cscl, cscr, csgni, csgnr, cspni, cspnr,  &
      csr, csrr, cssr, cyi, cyr, c1i, c1m, c1r, c2i, c2r, elim, fmr,   &
      fn, fnu, fnul, pi, pti, ptr, razn, rl, rzi, rzr, sc1i, sc1r,     &
      sc2i, sc2r, sgn, spn, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, &
      yy, zeroi, zeror, zi, zni, znr, zr
      INTEGER i, inu, iuf, kflag, kode, mr, n, nn, nw, nz
      dimension yr(n), yi(n), cyr(2), cyi(2), cssr(3), csrr(3), bry(3)
      DATA pi / 3.14159265358979324d0 /
      DATA zeror,zeroi,coner,conei / 0.0d0,0.0d0,1.0d0,0.0d0 /
      nz = 0
      znr = -zr
      zni = -zi
      nn = n
      CALL zbinu(znr, zni, fnu, kode, nn, yr, yi, nw, rl, fnul, tol, elim, alim)
      IF (nw.LT.0) GO TO 90
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!-----------------------------------------------------------------------
      nn = min0(2,n)
      CALL zbknu(znr, zni, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
      IF (nw.NE.0) GO TO 90
      s1r = cyr(1)
      s1i = cyi(1)
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
      csgnr = zeror
      csgni = sgn
      IF (kode.EQ.1) GO TO 10
      yy = -zni
      cpn = dcos(yy)
      spn = dsin(yy)
      CALL zmlt(csgnr, csgni, cpn, spn, csgnr, csgni)
   10 CONTINUE
!-----------------------------------------------------------------------
!     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
      inu = int(sngl(fnu))
      arg = (fnu-dble(float(inu)))*sgn
      cpn = dcos(arg)
      spn = dsin(arg)
      cspnr = cpn
      cspni = spn
      IF (mod(inu,2).EQ.0) GO TO 20
      cspnr = -cspnr
      cspni = -cspni
   20 CONTINUE
      iuf = 0
      c1r = s1r
      c1i = s1i
      c2r = yr(1)
      c2i = yi(1)
      ascle = 1.0d+3*d1mach(1)/tol
      IF (kode.EQ.1) GO TO 30
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
      sc1r = c1r
      sc1i = c1i
   30 CONTINUE
      CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
      CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
      yr(1) = str + ptr
      yi(1) = sti + pti
      IF (n.EQ.1) RETURN
      cspnr = -cspnr
      cspni = -cspni
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = s2r
      c1i = s2i
      c2r = yr(2)
      c2i = yi(2)
      IF (kode.EQ.1) GO TO 40
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
      sc2r = c1r
      sc2i = c1i
   40 CONTINUE
      CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
      CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
      yr(2) = str + ptr
      yi(2) = sti + pti
      IF (n.EQ.2) RETURN
      cspnr = -cspnr
      cspni = -cspni
      azn = zabs(znr,zni)
      razn = 1.0d0/azn
      str = znr*razn
      sti = -zni*razn
      rzr = (str+str)*razn
      rzi = (sti+sti)*razn
      fn = fnu + 1.0d0
      ckr = fn*rzr
      cki = fn*rzi
!-----------------------------------------------------------------------
!     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
!-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      cscr = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = cscr
      csrr(1) = cscr
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = ascle
      bry(2) = 1.0d0/ascle
      bry(3) = d1mach(2)
      as2 = zabs(s2r,s2i)
      kflag = 2
      IF (as2.GT.bry(1)) GO TO 50
      kflag = 1
      GO TO 60
   50 CONTINUE
      IF (as2.LT.bry(2)) GO TO 60
      kflag = 3
   60 CONTINUE
      bscle = bry(kflag)
      s1r = s1r*cssr(kflag)
      s1i = s1i*cssr(kflag)
      s2r = s2r*cssr(kflag)
      s2i = s2i*cssr(kflag)
      csr = csrr(kflag)
      DO 80 i=3,n
        str = s2r
        sti = s2i
        s2r = ckr*str - cki*sti + s1r
        s2i = ckr*sti + cki*str + s1i
        s1r = str
        s1i = sti
        c1r = s2r*csr
        c1i = s2i*csr
        str = c1r
        sti = c1i
        c2r = yr(i)
        c2i = yi(i)
        IF (kode.EQ.1) GO TO 70
        IF (iuf.LT.0) GO TO 70
        CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
        nz = nz + nw
        sc1r = sc2r
        sc1i = sc2i
        sc2r = c1r
        sc2i = c1i
        IF (iuf.NE.3) GO TO 70
        iuf = -4
        s1r = sc1r*cssr(kflag)
        s1i = sc1i*cssr(kflag)
        s2r = sc2r*cssr(kflag)
        s2i = sc2i*cssr(kflag)
        str = sc2r
        sti = sc2i
   70   CONTINUE
        ptr = cspnr*c1r - cspni*c1i
        pti = cspnr*c1i + cspni*c1r
        yr(i) = ptr + csgnr*c2r - csgni*c2i
        yi(i) = pti + csgnr*c2i + csgni*c2r
        ckr = ckr + rzr
        cki = cki + rzi
        cspnr = -cspnr
        cspni = -cspni
        IF (kflag.GE.3) GO TO 80
        ptr = dabs(c1r)
        pti = dabs(c1i)
        c1m = dmax1(ptr,pti)
        IF (c1m.LE.bscle) GO TO 80
        kflag = kflag + 1
        bscle = bry(kflag)
        s1r = s1r*csr
        s1i = s1i*csr
        s2r = str
        s2i = sti
        s1r = s1r*cssr(kflag)
        s1i = s1i*cssr(kflag)
        s2r = s2r*cssr(kflag)
        s2i = s2i*cssr(kflag)
        csr = csrr(kflag)
   80 CONTINUE
      RETURN
   90 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
END

SUBROUTINE zbunk(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)
!***BEGIN PROLOGUE  ZBUNK
!***REFER TO  ZBESK,ZBESH
!
!     ZBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU.GT.FNUL.
!     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR K(FNU,Z)
!     IN ZUNK1 AND THE EXPANSION FOR H(2,FNU,Z) IN ZUNK2
!
!***ROUTINES CALLED  ZUNK1,ZUNK2
!***END PROLOGUE  ZBUNK
!     COMPLEX Y,Z
      DOUBLE PRECISION alim, ax, ay, elim, fnu, tol, yi, yr, zi, zr
      INTEGER kode, mr, n, nz
      dimension yr(1), yi(1)
      nz = 0
      ax = dabs(zr)*1.7321d0
      ay = dabs(zi)
      IF (ay.GT.ax) GO TO 10
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR K(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3.LE.ARG(Z).LE.PI/3
!-----------------------------------------------------------------------
      CALL zunk1(zr, zi, fnu, kode, mr, n, yr, yi, nz, tol, elim, alim)
      GO TO 20
   10 CONTINUE
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR H(2,FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!-----------------------------------------------------------------------
      CALL zunk2(zr, zi, fnu, kode, mr, n, yr, yi, nz, tol, elim, alim)
   20 CONTINUE
      RETURN
END

SUBROUTINE zunik(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR,  &
       PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
USE complex
!***BEGIN PROLOGUE  ZUNIK
!***REFER TO  ZBESI,ZBESK
!
!        ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC
!        EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2
!        RESPECTIVELY BY
!
!        W(FNU,ZR) = PHI*EXP(ZETA)*SUM
!
!        WHERE       ZETA=-ZETA1 + ZETA2       OR
!                          ZETA1 - ZETA2
!
!        THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE
!        SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG=
!        1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK
!        ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI,
!        ZETA1,ZETA2.
!
!***ROUTINES CALLED  ZDIV,ZLOG,ZSQRT
!***END PROLOGUE  ZUNIK
!     COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,PHI,S,SR,SUM,T,T2,ZETA1,
!    *ZETA2,ZN,ZR
      DOUBLE PRECISION ac, c, con, conei, coner, crfni, crfnr, cwrki,  &
      cwrkr, fnu, phii, phir, rfn, si, sr, sri, srr, sti, str, sumi,   &
      sumr, test, ti, tol, tr, t2i, t2r, zeroi, zeror, zeta1i, zeta1r, &
      zeta2i, zeta2r, zni, znr, zri, zrr
      INTEGER i, idum, ikflg, init, ipmtr, j, k, l
      dimension c(120), cwrkr(16), cwrki(16), con(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
      DATA con(1), con(2)  /  &
       3.98942280401432678d-01,  1.25331413731550025d+00 /
      DATA c(1), c(2), c(3), c(4), c(5), c(6), c(7), c(8), c(9), c(10), &
           c(11), c(12), c(13), c(14), c(15), c(16), c(17), c(18),      &
           c(19), c(20), c(21), c(22), c(23), c(24)/                    &
           1.00000000000000000d+00,    -2.08333333333333333d-01,        &
           1.25000000000000000d-01,     3.34201388888888889d-01,        &
          -4.01041666666666667d-01,     7.03125000000000000d-02,        &
          -1.02581259645061728d+00,     1.84646267361111111d+00,        &
          -8.91210937500000000d-01,     7.32421875000000000d-02,        &
           4.66958442342624743d+00,    -1.12070026162229938d+01,        &
           8.78912353515625000d+00,    -2.36408691406250000d+00,        &
           1.12152099609375000d-01,    -2.82120725582002449d+01,        &
           8.46362176746007346d+01,    -9.18182415432400174d+01,        &
           4.25349987453884549d+01,    -7.36879435947963170d+00,        &
           2.27108001708984375d-01,     2.12570130039217123d+02,        &
          -7.65252468141181642d+02,     1.05999045252799988d+03/
      DATA c(25), c(26), c(27), c(28), c(29), c(30), c(31), c(32),      &
           c(33), c(34), c(35), c(36), c(37), c(38), c(39), c(40),      &
           c(41), c(42), c(43), c(44), c(45), c(46), c(47), c(48)/      &
          -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.72772750258445740d+00,        &
           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.07404200127348304d+00,    -2.42919187900551333d+05,        &
           1.31176361466297720d+06,    -2.99801591853810675d+06/
      DATA c(49), c(50), c(51), c(52), c(53), c(54), c(55), c(56),      &
           c(57), c(58), c(59), c(60), c(61), c(62), c(63), c(64),      &
           c(65), c(66), c(67), c(68), c(69), c(70), c(71), c(72)/      &
           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/
      DATA c(73), c(74), c(75), c(76), c(77), c(78), c(79), c(80),      &
           c(81), c(82), c(83), c(84), c(85), c(86), c(87), c(88),      &
           c(89), c(90), c(91), c(92), c(93), c(94), c(95), c(96)/      &
          -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/
      DATA c(97), c(98), c(99), c(100), c(101), c(102), c(103), c(104), &
           c(105), c(106), c(107), c(108), c(109), c(110), c(111),      &
           c(112), c(113), c(114), c(115), c(116), c(117), c(118)/      &
           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/
      DATA c(119), c(120)/  &
          -2.91883881222208134d+07,     1.18838426256783253d+05/

      IF (init.NE.0) GO TO 40
!-----------------------------------------------------------------------
!     INITIALIZE ALL VARIABLES
!-----------------------------------------------------------------------
      rfn = 1.0d0/fnu
      tr = zrr*rfn
      ti = zri*rfn
      sr = coner + (tr*tr-ti*ti)
      si = conei + (tr*ti+ti*tr)
      CALL zsqrt(sr, si, srr, sri)
      str = coner + srr
      sti = conei + sri
      CALL zdiv(str, sti, tr, ti, znr, zni)
      CALL zlog(znr, zni, str, sti, idum)
      zeta1r = fnu*str
      zeta1i = fnu*sti
      zeta2r = fnu*srr
      zeta2i = fnu*sri
      CALL zdiv(coner, conei, srr, sri, tr, ti)
      srr = tr*rfn
      sri = ti*rfn
      CALL zsqrt(srr, sri, cwrkr(16), cwrki(16))
      phir = cwrkr(16)*con(ikflg)
      phii = cwrki(16)*con(ikflg)
      IF (ipmtr.NE.0) RETURN
      CALL zdiv(coner, conei, sr, si, t2r, t2i)
      cwrkr(1) = coner
      cwrki(1) = conei
      crfnr = coner
      crfni = conei
      ac = 1.0d0
      l = 1
      DO 20 k=2,15
        sr = zeror
        si = zeroi
        DO 10 j=1,k
          l = l + 1
          str = sr*t2r - si*t2i + c(l)
          si = sr*t2i + si*t2r
          sr = str
   10   CONTINUE
        str = crfnr*srr - crfni*sri
        crfni = crfnr*sri + crfni*srr
        crfnr = str
        cwrkr(k) = crfnr*sr - crfni*si
        cwrki(k) = crfnr*si + crfni*sr
        ac = ac*rfn
        test = dabs(cwrkr(k)) + dabs(cwrki(k))
        IF (ac.LT.tol .AND. test.LT.tol) GO TO 30
   20 CONTINUE
      k = 15
   30 CONTINUE
      init = k
   40 CONTINUE
      IF (ikflg.EQ.2) GO TO 60
!-----------------------------------------------------------------------
!     COMPUTE SUM FOR THE I FUNCTION
!-----------------------------------------------------------------------
      sr = zeror
      si = zeroi
      DO 50 i=1,init
        sr = sr + cwrkr(i)
        si = si + cwrki(i)
   50 CONTINUE
      sumr = sr
      sumi = si
      phir = cwrkr(16)*con(1)
      phii = cwrki(16)*con(1)
      RETURN
   60 CONTINUE
!-----------------------------------------------------------------------
!     COMPUTE SUM FOR THE K FUNCTION
!-----------------------------------------------------------------------
      sr = zeror
      si = zeroi
      tr = coner
      DO 70 i=1,init
        sr = sr + tr*cwrkr(i)
        si = si + tr*cwrki(i)
        tr = -tr
   70 CONTINUE
      sumr = sr
      sumi = si
      phir = cwrkr(16)*con(2)
      phii = cwrki(16)*con(2)
      RETURN
END

SUBROUTINE zunk1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZUNK1
!***REFER TO  ZBESK
!
!     ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
!     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
!     UNIFORM ASYMPTOTIC EXPANSION.
!     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
!     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
!
!***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS
!***END PROLOGUE  ZUNK1
!     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
!    *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
      DOUBLE PRECISION alim, ang, aphi, asc, ascle, bry, cki, ckr,     &
      conei, coner, crsc, cscl, csgni, cspni, cspnr, csr, csrr, cssr,  &
      cwrki, cwrkr, cyi, cyr, c1i, c1r, c2i, c2m, c2r, elim, fmr, fn,  &
      fnf, fnu, phidi, phidr, phii, phir, pi, rast, razr, rs1, rzi,    &
      rzr, sgn, sti, str, sumdi, sumdr, sumi, sumr, s1i, s1r, s2i,     &
      s2r, tol, yi, yr, zeroi, zeror, zeta1i, zeta1r, zeta2i, zeta2r,  &
      zet1di, zet1dr, zet2di, zet2dr, zi, zr, zri, zrr
      INTEGER i, ib, iflag, ifn, il, init, inu, iuf, k, kdflg, kflag,  &
      kk, kode, mr, n, nw, nz, initd, ic, ipard, j, m
      dimension bry(3), init(2), yr(1), yi(1), sumr(2), sumi(2),       &
      zeta1r(2), zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2),      &
      cwrkr(16,3), cwrki(16,3), cssr(3), csrr(3), phir(2), phii(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
      DATA pi / 3.14159265358979324d0 /

      kdflg = 1
      nz = 0
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
!     THE UNDERFLOW LIMIT
!-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      crsc = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = crsc
      csrr(1) = crsc
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = 1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      zrr = zr
      zri = zi
      IF (zr.GE.0.0d0) GO TO 10
      zrr = -zr
      zri = -zi
   10 CONTINUE
      j = 2
      DO 70 i=1,n
!-----------------------------------------------------------------------
!     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!-----------------------------------------------------------------------
        j = 3 - j
        fn = fnu + dble(float(i-1))
        init(j) = 0
        CALL zunik(zrr, zri, fn, 2, 0, tol, init(j), phir(j), phii(j), &
        zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), sumr(j), sumi(j),  &
        cwrkr(1,j), cwrki(1,j))
        IF (kode.EQ.1) GO TO 20
        str = zrr + zeta2r(j)
        sti = zri + zeta2i(j)
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = zeta1r(j) - str
        s1i = zeta1i(j) - sti
        GO TO 30
   20   CONTINUE
        s1r = zeta1r(j) - zeta2r(j)
        s1i = zeta1i(j) - zeta2i(j)
   30   CONTINUE
        rs1 = s1r
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        IF (dabs(rs1).GT.elim) GO TO 60
        IF (kdflg.EQ.1) kflag = 2
        IF (dabs(rs1).LT.alim) GO TO 40
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phir(j),phii(j))
        rs1 = rs1 + dlog(aphi)
        IF (dabs(rs1).GT.elim) GO TO 60
        IF (kdflg.EQ.1) kflag = 1
        IF (rs1.LT.0.0d0) GO TO 40
        IF (kdflg.EQ.1) kflag = 3
   40   CONTINUE
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
        s2r = phir(j)*sumr(j) - phii(j)*sumi(j)
        s2i = phir(j)*sumi(j) + phii(j)*sumr(j)
        str = dexp(s1r)*cssr(kflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s1r*s2i + s2r*s1i
        s2r = str
        IF (kflag.NE.1) GO TO 50
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 60
   50   CONTINUE
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        yr(i) = s2r*csrr(kflag)
        yi(i) = s2i*csrr(kflag)
        IF (kdflg.EQ.2) GO TO 75
        kdflg = 2
        GO TO 70
   60   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 300
!-----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
        IF (zr.LT.0.0d0) GO TO 300
        kdflg = 1
        yr(i)=zeror
        yi(i)=zeroi
        nz=nz+1
        IF (i.EQ.1) GO TO 70
        IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi)) GO TO 70
        yr(i-1)=zeror
        yi(i-1)=zeroi
        nz=nz+1
   70 CONTINUE
      i = n
   75 CONTINUE
      razr = 1.0d0/zabs(zrr,zri)
      str = zrr*razr
      sti = -zri*razr
      rzr = (str+str)*razr
      rzi = (sti+sti)*razr
      ckr = fn*rzr
      cki = fn*rzi
      ib = i + 1
      IF (n.LT.ib) GO TO 160
!-----------------------------------------------------------------------
!     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
!     ON UNDERFLOW.
!-----------------------------------------------------------------------
      fn = fnu + dble(float(n-1))
      ipard = 1
      IF (mr.NE.0) ipard = 0
      initd = 0
      CALL zunik(zrr, zri, fn, 2, ipard, tol, initd, phidr, phidi, &
      zet1dr, zet1di, zet2dr, zet2di, sumdr, sumdi, cwrkr(1,3),    &
      cwrki(1,3))
      IF (kode.EQ.1) GO TO 80
      str = zrr + zet2dr
      sti = zri + zet2di
      rast = fn/zabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = zet1dr - str
      s1i = zet1di - sti
      GO TO 90
   80 CONTINUE
      s1r = zet1dr - zet2dr
      s1i = zet1di - zet2di
   90 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 95
      IF (dabs(rs1).LT.alim) GO TO 100
!-----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!-----------------------------------------------------------------------
      aphi = zabs(phidr,phidi)
      rs1 = rs1+dlog(aphi)
      IF (dabs(rs1).LT.elim) GO TO 100
   95 CONTINUE
      IF (dabs(rs1).GT.0.0d0) GO TO 300
!-----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
      IF (zr.LT.0.0d0) GO TO 300
      nz = n
      DO 96 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
   96 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
!-----------------------------------------------------------------------
  100 CONTINUE
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(kflag)
      ascle = bry(kflag)
      DO 120 i=ib,n
        c2r = s2r
        c2i = s2i
        s2r = ckr*c2r - cki*c2i + s1r
        s2i = ckr*c2i + cki*c2r + s1i
        s1r = c2r
        s1i = c2i
        ckr = ckr + rzr
        cki = cki + rzi
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(i) = c2r
        yi(i) = c2i
        IF (kflag.GE.3) GO TO 120
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 120
        kflag = kflag + 1
        ascle = bry(kflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(kflag)
        s1i = s1i*cssr(kflag)
        s2r = s2r*cssr(kflag)
        s2i = s2i*cssr(kflag)
        c1r = csrr(kflag)
  120 CONTINUE
  160 CONTINUE
      IF (mr.EQ.0) RETURN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
!-----------------------------------------------------------------------
      nz = 0
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
!-----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
!-----------------------------------------------------------------------
      csgni = sgn
      inu = int(sngl(fnu))
      fnf = fnu - dble(float(inu))
      ifn = inu + n - 1
      ang = fnf*sgn
      cspnr = dcos(ang)
      cspni = dsin(ang)
      IF (mod(ifn,2).EQ.0) GO TO 170
      cspnr = -cspnr
      cspni = -cspni
  170 CONTINUE
      asc = bry(1)
      iuf = 0
      kk = n
      kdflg = 1
      ib = ib - 1
      ic = ib - 1
      DO 270 k=1,n
        fn = fnu + dble(float(kk-1))
!-----------------------------------------------------------------------
!     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
!     FUNCTION ABOVE
!-----------------------------------------------------------------------
        m=3
        IF (n.GT.2) GO TO 175
  172   CONTINUE
        initd = init(j)
        phidr = phir(j)
        phidi = phii(j)
        zet1dr = zeta1r(j)
        zet1di = zeta1i(j)
        zet2dr = zeta2r(j)
        zet2di = zeta2i(j)
        sumdr = sumr(j)
        sumdi = sumi(j)
        m = j
        j = 3 - j
        GO TO 180
  175   CONTINUE
        IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 180
        IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 172
        initd = 0
  180   CONTINUE
        CALL zunik(zrr, zri, fn, 1, 0, tol, initd, phidr, phidi, &
        zet1dr, zet1di, zet2dr, zet2di, sumdr, sumdi,            &
        cwrkr(1,m), cwrki(1,m))
        IF (kode.EQ.1) GO TO 200
        str = zrr + zet2dr
        sti = zri + zet2di
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zet1dr + str
        s1i = -zet1di + sti
        GO TO 210
  200   CONTINUE
        s1r = -zet1dr + zet2dr
        s1i = -zet1di + zet2di
  210   CONTINUE
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 260
        IF (kdflg.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 220
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phidr,phidi)
        rs1 = rs1 + dlog(aphi)
        IF (dabs(rs1).GT.elim) GO TO 260
        IF (kdflg.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 220
        IF (kdflg.EQ.1) iflag = 3
  220   CONTINUE
        str = phidr*sumdr - phidi*sumdi
        sti = phidr*sumdi + phidi*sumdr
        s2r = -csgni*sti
        s2i = csgni*str
        str = dexp(s1r)*cssr(iflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s2r*s1i + s2i*s1r
        s2r = str
        IF (iflag.NE.1) GO TO 230
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.EQ.0) GO TO 230
        s2r = zeror
        s2i = zeroi
  230   CONTINUE
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        c2r = s2r
        c2i = s2i
        s2r = s2r*csrr(iflag)
        s2i = s2i*csrr(iflag)
!-----------------------------------------------------------------------
!     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!-----------------------------------------------------------------------
        s1r = yr(kk)
        s1i = yi(kk)
        IF (kode.EQ.1) GO TO 250
        CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
        nz = nz + nw
  250   CONTINUE
        yr(kk) = s1r*cspnr - s1i*cspni + s2r
        yi(kk) = cspnr*s1i + cspni*s1r + s2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0) GO TO 255
        kdflg = 1
        GO TO 270
  255   CONTINUE
        IF (kdflg.EQ.2) GO TO 275
        kdflg = 2
        GO TO 270
  260   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 300
        s2r = zeror
        s2i = zeroi
        GO TO 230
  270 CONTINUE
      k = n
  275 CONTINUE
      il = n - k
      IF (il.EQ.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.
!-----------------------------------------------------------------------
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      csr = csrr(iflag)
      ascle = bry(iflag)
      fn = dble(float(inu+il))
      DO 290 i=1,il
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        fn = fn - 1.0d0
        c2r = s2r*csr
        c2i = s2i*csr
        ckr = c2r
        cki = c2i
        c1r = yr(kk)
        c1i = yi(kk)
        IF (kode.EQ.1) GO TO 280
        CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
        nz = nz + nw
  280   CONTINUE
        yr(kk) = c1r*cspnr - c1i*cspni + c2r
        yi(kk) = c1r*cspni + c1i*cspnr + c2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        IF (iflag.GE.3) GO TO 290
        c2r = dabs(ckr)
        c2i = dabs(cki)
        c2m = dmax1(c2r,c2i)
        IF (c2m.LE.ascle) GO TO 290
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*csr
        s1i = s1i*csr
        s2r = ckr
        s2i = cki
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        csr = csrr(iflag)
  290 CONTINUE
      RETURN
  300 CONTINUE
      nz = -1
      RETURN
END

SUBROUTINE zunk2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZUNK2
!***REFER TO  ZBESK
!
!     ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
!     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
!     UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
!     WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
!     -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
!     HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
!     ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
!     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
!
!***ROUTINES CALLED  ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS
!***END PROLOGUE  ZUNK2
!     COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
!    *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
!    *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
      DOUBLE PRECISION aarg, aic, aii, air, alim, ang, aphi, argdi,   &
      argdr, argi, argr, asc, ascle, asumdi, asumdr, asumi, asumr,    &
      bry, bsumdi, bsumdr, bsumi, bsumr, car, cipi, cipr, cki, ckr,   &
      conei, coner, crsc, cr1i, cr1r, cr2i, cr2r, cscl, csgni, csi,   &
      cspni, cspnr, csr, csrr, cssr, cyi, cyr, c1i, c1r, c2i, c2m,    &
      c2r, daii, dair, elim, fmr, fn, fnf, fnu, hpi, phidi, phidr,    &
      phii, phir, pi, pti, ptr, rast, razr, rs1, rzi, rzr, sar, sgn,  &
      sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, yy, zbi, zbr, zeroi, &
      zeror, zeta1i, zeta1r, zeta2i, zeta2r, zet1di, zet1dr, zet2di,  &
      zet2dr, zi, zni, znr, zr, zri, zrr
      INTEGER i, ib, iflag, ifn, il, in, inu, iuf, k, kdflg, kflag, kk, &
      kode, mr, n, nai, ndai, nw, nz, idum, j, ipard, ic
      dimension bry(3), yr(1), yi(1), asumr(2), asumi(2), bsumr(2),   &
      bsumi(2), phir(2), phii(2), argr(2), argi(2), zeta1r(2),        &
      zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2), cipr(4),       &
      cipi(4), cssr(3), csrr(3)
      DATA zeror,zeroi,coner,conei,cr1r,cr1i,cr2r,cr2i /              &
               0.0d0, 0.0d0, 1.0d0, 0.0d0,                            &
       1.0d0,1.73205080756887729d0 , -0.5d0,-8.66025403784438647d-01 /
      DATA hpi, pi, aic /                                             &
           1.57079632679489662d+00,     3.14159265358979324d+00,      &
           1.26551212348464539d+00/
      DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),   &
       cipi(4) /                                                      &
       1.0d0,0.0d0 ,  0.0d0,-1.0d0 ,  -1.0d0,0.0d0 ,  0.0d0,1.0d0 /

      kdflg = 1
      nz = 0
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
!     THE UNDERFLOW LIMIT
!-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      crsc = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = crsc
      csrr(1) = crsc
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = 1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      zrr = zr
      zri = zi
      IF (zr.GE.0.0d0) GO TO 10
      zrr = -zr
      zri = -zi
   10 CONTINUE
      yy = zri
      znr = zri
      zni = -zrr
      zbr = zrr
      zbi = zri
      inu = int(sngl(fnu))
      fnf = fnu - dble(float(inu))
      ang = -hpi*fnf
      car = dcos(ang)
      sar = dsin(ang)
      c2r = hpi*sar
      c2i = -hpi*car
      kk = mod(inu,4) + 1
      str = c2r*cipr(kk) - c2i*cipi(kk)
      sti = c2r*cipi(kk) + c2i*cipr(kk)
      csr = cr1r*str - cr1i*sti
      csi = cr1r*sti + cr1i*str
      IF (yy.GT.0.0d0) GO TO 20
      znr = -znr
      zbi = -zbi
   20 CONTINUE
!-----------------------------------------------------------------------
!     K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
!     QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
!     CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
!-----------------------------------------------------------------------
      j = 2
      DO 80 i=1,n
!-----------------------------------------------------------------------
!     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!-----------------------------------------------------------------------
        j = 3 - j
        fn = fnu + dble(float(i-1))
        CALL zunhj(znr, zni, fn, 0, tol, phir(j), phii(j), argr(j),     &
        argi(j), zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), asumr(j),  &
        asumi(j), bsumr(j), bsumi(j))
        IF (kode.EQ.1) GO TO 30
        str = zbr + zeta2r(j)
        sti = zbi + zeta2i(j)
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = zeta1r(j) - str
        s1i = zeta1i(j) - sti
        GO TO 40
   30   CONTINUE
        s1r = zeta1r(j) - zeta2r(j)
        s1i = zeta1i(j) - zeta2i(j)
   40   CONTINUE
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 70
        IF (kdflg.EQ.1) kflag = 2
        IF (dabs(rs1).LT.alim) GO TO 50
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phir(j),phii(j))
        aarg = zabs(argr(j),argi(j))
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 70
        IF (kdflg.EQ.1) kflag = 1
        IF (rs1.LT.0.0d0) GO TO 50
        IF (kdflg.EQ.1) kflag = 3
   50   CONTINUE
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
        c2r = argr(j)*cr2r - argi(j)*cr2i
        c2i = argr(j)*cr2i + argi(j)*cr2r
        CALL zairy(c2r, c2i, 0, 2, air, aii, nai, idum)
        CALL zairy(c2r, c2i, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumr(j) - daii*bsumi(j)
        sti = dair*bsumi(j) + daii*bsumr(j)
        ptr = str*cr2r - sti*cr2i
        pti = str*cr2i + sti*cr2r
        str = ptr + (air*asumr(j)-aii*asumi(j))
        sti = pti + (air*asumi(j)+aii*asumr(j))
        ptr = str*phir(j) - sti*phii(j)
        pti = str*phii(j) + sti*phir(j)
        s2r = ptr*csr - pti*csi
        s2i = ptr*csi + pti*csr
        str = dexp(s1r)*cssr(kflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s1r*s2i + s2r*s1i
        s2r = str
        IF (kflag.NE.1) GO TO 60
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 70
   60   CONTINUE
        IF (yy.LE.0.0d0) s2i = -s2i
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        yr(i) = s2r*csrr(kflag)
        yi(i) = s2i*csrr(kflag)
        str = csi
        csi = -csr
        csr = str
        IF (kdflg.EQ.2) GO TO 85
        kdflg = 2
        GO TO 80
   70   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 320
!-----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
        IF (zr.LT.0.0d0) GO TO 320
        kdflg = 1
        yr(i)=zeror
        yi(i)=zeroi
        nz=nz+1
        str = csi
        csi =-csr
        csr = str
        IF (i.EQ.1) GO TO 80
        IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi)) GO TO 80
        yr(i-1)=zeror
        yi(i-1)=zeroi
        nz=nz+1
   80 CONTINUE
      i = n
   85 CONTINUE
      razr = 1.0d0/zabs(zrr,zri)
      str = zrr*razr
      sti = -zri*razr
      rzr = (str+str)*razr
      rzi = (sti+sti)*razr
      ckr = fn*rzr
      cki = fn*rzi
      ib = i + 1
      IF (n.LT.ib) GO TO 180
!-----------------------------------------------------------------------
!     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
!     ON UNDERFLOW.
!-----------------------------------------------------------------------
      fn = fnu + dble(float(n-1))
      ipard = 1
      IF (mr.NE.0) ipard = 0
      CALL zunhj(znr, zni, fn, ipard, tol, phidr, phidi, argdr, argdi, &
      zet1dr, zet1di, zet2dr, zet2di, asumdr, asumdi, bsumdr, bsumdi)
      IF (kode.EQ.1) GO TO 90
      str = zbr + zet2dr
      sti = zbi + zet2di
      rast = fn/zabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = zet1dr - str
      s1i = zet1di - sti
      GO TO 100
   90 CONTINUE
      s1r = zet1dr - zet2dr
      s1i = zet1di - zet2di
  100 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 105
      IF (dabs(rs1).LT.alim) GO TO 120
!-----------------------------------------------------------------------
!     REFINE ESTIMATE AND TEST
!-----------------------------------------------------------------------
      aphi = zabs(phidr,phidi)
      rs1 = rs1+dlog(aphi)
      IF (dabs(rs1).LT.elim) GO TO 120
  105 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 320
!-----------------------------------------------------------------------
!     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!-----------------------------------------------------------------------
      IF (zr.LT.0.0d0) GO TO 320
      nz = n
      DO 106 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  106 CONTINUE
      RETURN
  120 CONTINUE
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(kflag)
      ascle = bry(kflag)
      DO 130 i=ib,n
        c2r = s2r
        c2i = s2i
        s2r = ckr*c2r - cki*c2i + s1r
        s2i = ckr*c2i + cki*c2r + s1i
        s1r = c2r
        s1i = c2i
        ckr = ckr + rzr
        cki = cki + rzi
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(i) = c2r
        yi(i) = c2i
        IF (kflag.GE.3) GO TO 130
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 130
        kflag = kflag + 1
        ascle = bry(kflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(kflag)
        s1i = s1i*cssr(kflag)
        s2r = s2r*cssr(kflag)
        s2i = s2i*cssr(kflag)
        c1r = csrr(kflag)
  130 CONTINUE
  180 CONTINUE
      IF (mr.EQ.0) RETURN
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
!-----------------------------------------------------------------------
      nz = 0
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
!-----------------------------------------------------------------------
!     CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
!-----------------------------------------------------------------------
      csgni = sgn
      IF (yy.LE.0.0d0) csgni = -csgni
      ifn = inu + n - 1
      ang = fnf*sgn
      cspnr = dcos(ang)
      cspni = dsin(ang)
      IF (mod(ifn,2).EQ.0) GO TO 190
      cspnr = -cspnr
      cspni = -cspni
  190 CONTINUE
!-----------------------------------------------------------------------
!     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.LE.0.0E0) ARE COMPUTED BY
!     CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
!-----------------------------------------------------------------------
      csr = sar*csgni
      csi = car*csgni
      in = mod(ifn,4) + 1
      c2r = cipr(in)
      c2i = cipi(in)
      str = csr*c2r + csi*c2i
      csi = -csr*c2i + csi*c2r
      csr = str
      asc = bry(1)
      iuf = 0
      kk = n
      kdflg = 1
      ib = ib - 1
      ic = ib - 1
      DO 290 k=1,n
        fn = fnu + dble(float(kk-1))
!-----------------------------------------------------------------------
!     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
!     FUNCTION ABOVE
!-----------------------------------------------------------------------
        IF (n.GT.2) GO TO 175
  172   CONTINUE
        phidr = phir(j)
        phidi = phii(j)
        argdr = argr(j)
        argdi = argi(j)
        zet1dr = zeta1r(j)
        zet1di = zeta1i(j)
        zet2dr = zeta2r(j)
        zet2di = zeta2i(j)
        asumdr = asumr(j)
        asumdi = asumi(j)
        bsumdr = bsumr(j)
        bsumdi = bsumi(j)
        j = 3 - j
        GO TO 210
  175   CONTINUE
        IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 210
        IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 172
        CALL zunhj(znr, zni, fn, 0, tol, phidr, phidi, argdr,  &
        argdi, zet1dr, zet1di, zet2dr, zet2di, asumdr,         &
        asumdi, bsumdr, bsumdi)
  210   CONTINUE
        IF (kode.EQ.1) GO TO 220
        str = zbr + zet2dr
        sti = zbi + zet2di
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zet1dr + str
        s1i = -zet1di + sti
        GO TO 230
  220   CONTINUE
        s1r = -zet1dr + zet2dr
        s1i = -zet1di + zet2di
  230   CONTINUE
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 280
        IF (kdflg.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 240
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phidr,phidi)
        aarg = zabs(argdr,argdi)
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 280
        IF (kdflg.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 240
        IF (kdflg.EQ.1) iflag = 3
  240   CONTINUE
        CALL zairy(argdr, argdi, 0, 2, air, aii, nai, idum)
        CALL zairy(argdr, argdi, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumdr - daii*bsumdi
        sti = dair*bsumdi + daii*bsumdr
        str = str + (air*asumdr-aii*asumdi)
        sti = sti + (air*asumdi+aii*asumdr)
        ptr = str*phidr - sti*phidi
        pti = str*phidi + sti*phidr
        s2r = ptr*csr - pti*csi
        s2i = ptr*csi + pti*csr
        str = dexp(s1r)*cssr(iflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s2r*s1i + s2i*s1r
        s2r = str
        IF (iflag.NE.1) GO TO 250
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.EQ.0) GO TO 250
        s2r = zeror
        s2i = zeroi
  250   CONTINUE
        IF (yy.LE.0.0d0) s2i = -s2i
        cyr(kdflg) = s2r
        cyi(kdflg) = s2i
        c2r = s2r
        c2i = s2i
        s2r = s2r*csrr(iflag)
        s2i = s2i*csrr(iflag)
!-----------------------------------------------------------------------
!     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!-----------------------------------------------------------------------
        s1r = yr(kk)
        s1i = yi(kk)
        IF (kode.EQ.1) GO TO 270
        CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
        nz = nz + nw
  270   CONTINUE
        yr(kk) = s1r*cspnr - s1i*cspni + s2r
        yi(kk) = s1r*cspni + s1i*cspnr + s2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        str = csi
        csi = -csr
        csr = str
        IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0) GO TO 255
        kdflg = 1
        GO TO 290
  255   CONTINUE
        IF (kdflg.EQ.2) GO TO 295
        kdflg = 2
        GO TO 290
  280   CONTINUE
        IF (rs1.GT.0.0d0) GO TO 320
        s2r = zeror
        s2i = zeroi
        GO TO 250
  290 CONTINUE
      k = n
  295 CONTINUE
      il = n - k
      IF (il.EQ.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.
!-----------------------------------------------------------------------
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      csr = csrr(iflag)
      ascle = bry(iflag)
      fn = dble(float(inu+il))
      DO 310 i=1,il
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        fn = fn - 1.0d0
        c2r = s2r*csr
        c2i = s2i*csr
        ckr = c2r
        cki = c2i
        c1r = yr(kk)
        c1i = yi(kk)
        IF (kode.EQ.1) GO TO 300
        CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
        nz = nz + nw
  300   CONTINUE
        yr(kk) = c1r*cspnr - c1i*cspni + c2r
        yi(kk) = c1r*cspni + c1i*cspnr + c2i
        kk = kk - 1
        cspnr = -cspnr
        cspni = -cspni
        IF (iflag.GE.3) GO TO 310
        c2r = dabs(ckr)
        c2i = dabs(cki)
        c2m = dmax1(c2r,c2i)
        IF (c2m.LE.ascle) GO TO 310
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*csr
        s1i = s1i*csr
        s2r = ckr
        s2i = cki
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        csr = csrr(iflag)
  310 CONTINUE
      RETURN
  320 CONTINUE
      nz = -1
      RETURN
END

SUBROUTINE zunhj(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI,  &
       ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
USE complex
!***BEGIN PROLOGUE  ZUNHJ
!***REFER TO  ZBESI,ZBESK
!
!     REFERENCES
!         HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A.
!         STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9.
!
!         ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC
!         PRESS, N.Y., 1974, PAGE 420
!
!     ABSTRACT
!         ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) =
!         J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU
!         BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION
!
!         C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) )
!
!         FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS
!         AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE.
!
!               (2/3)*FNU*ZETA**1.5 = ZETA1-ZETA2,
!
!         ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING
!         PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY.
!
!         MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND
!         MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR=
!         1 COMPUTES ALL EXCEPT ASUM AND BSUM.
!
!***ROUTINES CALLED  ZABS,ZDIV,ZLOG,ZSQRT
!***END PROLOGUE  ZUNHJ
!     COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN,
!    *RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1,
!    *ZETA2,ZTH
      DOUBLE PRECISION alfa, ang, ap, ar, argi, argr, asumi, asumr,      &
       atol, aw2, azth, beta, br, bsumi, bsumr, btol, c, conei, coner,   &
       cri, crr, dri, drr, ex1, ex2, fnu, fn13, fn23, gama, gpi, hpi,    &
       phii, phir, pi, pp, pr, przthi, przthr, ptfni, ptfnr, raw, raw2,  &
       razth, rfnu, rfnu2, rfn13, rtzti, rtztr, rzthi, rzthr, sti, str,  &
       sumai, sumar, sumbi, sumbr, test, tfni, tfnr, thpi, tol, tzai,    &
       tzar, t2i, t2r, upi, upr, wi, wr, w2i, w2r, zai, zar, zbi, zbr,   &
       zci, zcr, zeroi, zeror, zetai, zetar, zeta1i, zeta1r, zeta2i,     &
       zeta2r, zi, zr, zthi, zthr
      INTEGER ias, ibs, ipmtr, is, j, jr, ju, k, kmax, kp1, ks, l, lr,   &
       lrp1, l1, l2, m, idum
      dimension ar(14), br(14), c(105), alfa(180), beta(210), gama(30),  &
       ap(30), pr(30), pi(30), upr(14), upi(14), crr(14), cri(14),       &
       drr(14), dri(14)
      DATA ar(1), ar(2), ar(3), ar(4), ar(5), ar(6), ar(7), ar(8),  &
           ar(9), ar(10), ar(11), ar(12), ar(13), ar(14)/           &
           1.00000000000000000d+00,     1.04166666666666667d-01,    &
           8.35503472222222222d-02,     1.28226574556327160d-01,    &
           2.91849026464140464d-01,     8.81627267443757652d-01,    &
           3.32140828186276754d+00,     1.49957629868625547d+01,    &
           7.89230130115865181d+01,     4.74451538868264323d+02,    &
           3.20749009089066193d+03,     2.40865496408740049d+04,    &
           1.98923119169509794d+05,     1.79190200777534383d+06/
      DATA br(1), br(2), br(3), br(4), br(5), br(6), br(7), br(8),  &
           br(9), br(10), br(11), br(12), br(13), br(14)/           &
           1.00000000000000000d+00,    -1.45833333333333333d-01,    &
          -9.87413194444444444d-02,    -1.43312053915895062d-01,    &
          -3.17227202678413548d-01,    -9.42429147957120249d-01,    &
          -3.51120304082635426d+00,    -1.57272636203680451d+01,    &
          -8.22814390971859444d+01,    -4.92355370523670524d+02,    &
          -3.31621856854797251d+03,    -2.48276742452085896d+04,    &
          -2.04526587315129788d+05,    -1.83844491706820990d+06/
      DATA c(1), c(2), c(3), c(4), c(5), c(6), c(7), c(8), c(9), c(10), &
           c(11), c(12), c(13), c(14), c(15), c(16), c(17), c(18),      &
           c(19), c(20), c(21), c(22), c(23), c(24)/                    &
           1.00000000000000000d+00,    -2.08333333333333333d-01,        &
           1.25000000000000000d-01,     3.34201388888888889d-01,        &
          -4.01041666666666667d-01,     7.03125000000000000d-02,        &
          -1.02581259645061728d+00,     1.84646267361111111d+00,        &
          -8.91210937500000000d-01,     7.32421875000000000d-02,        &
           4.66958442342624743d+00,    -1.12070026162229938d+01,        &
           8.78912353515625000d+00,    -2.36408691406250000d+00,        &
           1.12152099609375000d-01,    -2.82120725582002449d+01,        &
           8.46362176746007346d+01,    -9.18182415432400174d+01,        &
           4.25349987453884549d+01,    -7.36879435947963170d+00,        &
           2.27108001708984375d-01,     2.12570130039217123d+02,        &
          -7.65252468141181642d+02,     1.05999045252799988d+03/
      DATA c(25), c(26), c(27), c(28), c(29), c(30), c(31), c(32),      &
           c(33), c(34), c(35), c(36), c(37), c(38), c(39), c(40),      &
           c(41), c(42), c(43), c(44), c(45), c(46), c(47), c(48)/      &
          -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.72772750258445740d+00,        &
           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.07404200127348304d+00,    -2.42919187900551333d+05,        &
           1.31176361466297720d+06,    -2.99801591853810675d+06/
      DATA c(49), c(50), c(51), c(52), c(53), c(54), c(55), c(56),      &
           c(57), c(58), c(59), c(60), c(61), c(62), c(63), c(64),      &
           c(65), c(66), c(67), c(68), c(69), c(70), c(71), c(72)/      &
           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/
      DATA c(73), c(74), c(75), c(76), c(77), c(78), c(79), c(80),      &
           c(81), c(82), c(83), c(84), c(85), c(86), c(87), c(88),      &
           c(89), c(90), c(91), c(92), c(93), c(94), c(95), c(96)/      &
          -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/
      DATA c(97), c(98), c(99), c(100), c(101), c(102), c(103), c(104), &
           c(105)/                                                      &
           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/
      DATA alfa(1), alfa(2), alfa(3), alfa(4), alfa(5), alfa(6),        &
           alfa(7), alfa(8), alfa(9), alfa(10), alfa(11), alfa(12),     &
           alfa(13), alfa(14), alfa(15), alfa(16), alfa(17), alfa(18),  &
           alfa(19), alfa(20), alfa(21), alfa(22)/                      &
          -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/
      DATA alfa(23), alfa(24), alfa(25), alfa(26), alfa(27), alfa(28),  &
           alfa(29), alfa(30), alfa(31), alfa(32), alfa(33), alfa(34),  &
           alfa(35), alfa(36), alfa(37), alfa(38), alfa(39), alfa(40),  &
           alfa(41), alfa(42), alfa(43), alfa(44)/                      &
           9.32120517249503256d-05,     8.85710852478711718d-05,        &
           8.42963105715700223d-05,     8.03497548407791151d-05,        &
           7.66981345359207388d-05,     7.33122157481777809d-05,        &
           7.01662625163141333d-05,     6.72375633790160292d-05,        &
           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/
      DATA alfa(45), alfa(46), alfa(47), alfa(48), alfa(49), alfa(50),  &
           alfa(51), alfa(52), alfa(53), alfa(54), alfa(55), alfa(56),  &
           alfa(57), alfa(58), alfa(59), alfa(60), alfa(61), alfa(62),  &
           alfa(63), alfa(64), alfa(65), alfa(66)/                      &
          -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,        &
          -3.54211971457743841d-04,    -1.56161263945159416d-04,        &
           3.04465503594936410d-05,     1.30198655773242693d-04,        &
           1.67471106699712269d-04,     1.70222587683592569d-04/
      DATA alfa(67), alfa(68), alfa(69), alfa(70), alfa(71), alfa(72),  &
           alfa(73), alfa(74), alfa(75), alfa(76), alfa(77), alfa(78),  &
           alfa(79), alfa(80), alfa(81), alfa(82), alfa(83), alfa(84),  &
           alfa(85), alfa(86), alfa(87), alfa(88)/                      &
           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/
      DATA alfa(89), alfa(90), alfa(91), alfa(92), alfa(93), alfa(94),  &
           alfa(95), alfa(96), alfa(97), alfa(98), alfa(99), alfa(100), &
           alfa(101), alfa(102), alfa(103), alfa(104), alfa(105),       &
           alfa(106), alfa(107), alfa(108), alfa(109), alfa(110)/       &
          -1.17519675674556414d-05,    -1.19987364870944141d-05,        &
           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/
      DATA alfa(111), alfa(112), alfa(113), alfa(114), alfa(115),       &
           alfa(116), alfa(117), alfa(118), alfa(119), alfa(120),       &
           alfa(121), alfa(122), alfa(123), alfa(124), alfa(125),       &
           alfa(126), alfa(127), alfa(128), alfa(129), alfa(130)/       &
           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,        &
          -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/
      DATA alfa(131), alfa(132), alfa(133), alfa(134), alfa(135),       &
           alfa(136), alfa(137), alfa(138), alfa(139), alfa(140),       &
           alfa(141), alfa(142), alfa(143), alfa(144), alfa(145),       &
           alfa(146), alfa(147), alfa(148), alfa(149), alfa(150)/       &
           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/
      DATA alfa(151), alfa(152), alfa(153), alfa(154), alfa(155),       &
           alfa(156), alfa(157), alfa(158), alfa(159), alfa(160),       &
           alfa(161), alfa(162), alfa(163), alfa(164), alfa(165),       &
           alfa(166), alfa(167), alfa(168), alfa(169), alfa(170)/       &
           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/
      DATA alfa(171), alfa(172), alfa(173), alfa(174), alfa(175),       &
           alfa(176), alfa(177), alfa(178), alfa(179), alfa(180)/       &
           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/
      DATA beta(1), beta(2), beta(3), beta(4), beta(5), beta(6),        &
           beta(7), beta(8), beta(9), beta(10), beta(11), beta(12),     &
           beta(13), beta(14), beta(15), beta(16), beta(17), beta(18),  &
           beta(19), beta(20), beta(21), beta(22)/                      &
           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/
      DATA beta(23), beta(24), beta(25), beta(26), beta(27), beta(28),  &
           beta(29), beta(30), beta(31), beta(32), beta(33), beta(34),  &
           beta(35), beta(36), beta(37), beta(38), beta(39), beta(40),  &
           beta(41), beta(42), beta(43), beta(44)/                      &
           9.41998224204237509d-05,     8.74130545753834437d-05,        &
           8.13466262162801467d-05,     7.59002269646219339d-05,        &
           7.09906300634153481d-05,     6.65482874842468183d-05,        &
           6.25146958969275078d-05,     5.88403394426251749d-05,        &
          -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/
      DATA beta(45), beta(46), beta(47), beta(48), beta(49), beta(50),  &
           beta(51), beta(52), beta(53), beta(54), beta(55), beta(56),  &
           beta(57), beta(58), beta(59), beta(60), beta(61), beta(62),  &
           beta(63), beta(64), beta(65), beta(66)/                      &
           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,        &
           5.52213076721292790d-04,     4.47932581552384646d-04,        &
           2.79520653992020589d-04,     1.52468156198446602d-04,        &
           6.93271105657043598d-05,     1.76258683069991397d-05/
      DATA beta(67), beta(68), beta(69), beta(70), beta(71), beta(72),  &
           beta(73), beta(74), beta(75), beta(76), beta(77), beta(78),  &
           beta(79), beta(80), beta(81), beta(82), beta(83), beta(84),  &
           beta(85), beta(86), beta(87), beta(88)/                      &
          -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/
      DATA beta(89), beta(90), beta(91), beta(92), beta(93), beta(94),  &
           beta(95), beta(96), beta(97), beta(98), beta(99), beta(100), &
           beta(101), beta(102), beta(103), beta(104), beta(105),       &
           beta(106), beta(107), beta(108), beta(109), beta(110)/       &
          -1.74532699844210224d-05,    -1.63997823854497997d-05,        &
          -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/
      DATA beta(111), beta(112), beta(113), beta(114), beta(115),       &
           beta(116), beta(117), beta(118), beta(119), beta(120),       &
           beta(121), beta(122), beta(123), beta(124), beta(125),       &
           beta(126), beta(127), beta(128), beta(129), beta(130)/       &
           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,        &
           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/
      DATA beta(131), beta(132), beta(133), beta(134), beta(135),       &
           beta(136), beta(137), beta(138), beta(139), beta(140),       &
           beta(141), beta(142), beta(143), beta(144), beta(145),       &
           beta(146), beta(147), beta(148), beta(149), beta(150)/       &
          -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/
      DATA beta(151), beta(152), beta(153), beta(154), beta(155),       &
           beta(156), beta(157), beta(158), beta(159), beta(160),       &
           beta(161), beta(162), beta(163), beta(164), beta(165),       &
           beta(166), beta(167), beta(168), beta(169), beta(170)/       &
          -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/
      DATA beta(171), beta(172), beta(173), beta(174), beta(175),       &
           beta(176), beta(177), beta(178), beta(179), beta(180),       &
           beta(181), beta(182), beta(183), beta(184), beta(185),       &
           beta(186), beta(187), beta(188), beta(189), beta(190)/       &
          -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,        &
           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/
      DATA beta(191), beta(192), beta(193), beta(194), beta(195),       &
           beta(196), beta(197), beta(198), beta(199), beta(200),       &
           beta(201), beta(202), beta(203), beta(204), beta(205),       &
           beta(206), beta(207), beta(208), beta(209), beta(210)/       &
          -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/
      DATA gama(1), gama(2), gama(3), gama(4), gama(5), gama(6),        &
           gama(7), gama(8), gama(9), gama(10), gama(11), gama(12),     &
           gama(13), gama(14), gama(15), gama(16), gama(17), gama(18),  &
           gama(19), gama(20), gama(21), gama(22)/                      &
           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/
      DATA gama(23), gama(24), gama(25), gama(26), gama(27), gama(28),  &
           gama(29), gama(30)/                                          &
           1.57865285987918445d-02,     1.50729501494095594d-02,        &
           1.44193250839954639d-02,     1.38184805735341786d-02,        &
           1.32643378994276568d-02,     1.27517121970498651d-02,        &
           1.22761545318762767d-02,     1.18338262398482403d-02/
      DATA ex1, ex2, hpi, gpi, thpi /                                   &
           3.33333333333333333d-01,     6.66666666666666667d-01,        &
           1.57079632679489662d+00,     3.14159265358979324d+00,        &
           4.71238898038468986d+00/
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /

      rfnu = 1.0d0/fnu
      zbr = zr*rfnu
      zbi = zi*rfnu
      rfnu2 = rfnu*rfnu
!-----------------------------------------------------------------------
!     COMPUTE IN THE FOURTH QUADRANT
!-----------------------------------------------------------------------
      fn13 = fnu**ex1
      fn23 = fn13*fn13
      rfn13 = 1.0d0/fn13
      w2r = coner - zbr*zbr + zbi*zbi
      w2i = conei - zbr*zbi - zbr*zbi
      aw2 = zabs(w2r,w2i)
      IF (aw2.GT.0.25d0) GO TO 130
!-----------------------------------------------------------------------
!     POWER SERIES FOR CABS(W2).LE.0.25D0
!-----------------------------------------------------------------------
      k = 1
      pr(1) = coner
      pi(1) = conei
      sumar = gama(1)
      sumai = zeroi
      ap(1) = 1.0d0
      IF (aw2.LT.tol) GO TO 20
      DO 10 k=2,30
        pr(k) = pr(k-1)*w2r - pi(k-1)*w2i
        pi(k) = pr(k-1)*w2i + pi(k-1)*w2r
        sumar = sumar + pr(k)*gama(k)
        sumai = sumai + pi(k)*gama(k)
        ap(k) = ap(k-1)*aw2
        IF (ap(k).LT.tol) GO TO 20
   10 CONTINUE
      k = 30
   20 CONTINUE
      kmax = k
      zetar = w2r*sumar - w2i*sumai
      zetai = w2r*sumai + w2i*sumar
      argr = zetar*fn23
      argi = zetai*fn23
      CALL zsqrt(sumar, sumai, zar, zai)
      CALL zsqrt(w2r, w2i, str, sti)
      zeta2r = str*fnu
      zeta2i = sti*fnu
      str = coner + ex2*(zetar*zar-zetai*zai)
      sti = conei + ex2*(zetar*zai+zetai*zar)
      zeta1r = str*zeta2r - sti*zeta2i
      zeta1i = str*zeta2i + sti*zeta2r
      zar = zar + zar
      zai = zai + zai
      CALL zsqrt(zar, zai, str, sti)
      phir = str*rfn13
      phii = sti*rfn13
      IF (ipmtr.EQ.1) GO TO 120
!-----------------------------------------------------------------------
!     SUM SERIES FOR ASUM AND BSUM
!-----------------------------------------------------------------------
      sumbr = zeror
      sumbi = zeroi
      DO 30 k=1,kmax
        sumbr = sumbr + pr(k)*beta(k)
        sumbi = sumbi + pi(k)*beta(k)
   30 CONTINUE
      asumr = zeror
      asumi = zeroi
      bsumr = sumbr
      bsumi = sumbi
      l1 = 0
      l2 = 30
      btol = tol*(dabs(bsumr)+dabs(bsumi))
      atol = tol
      pp = 1.0d0
      ias = 0
      ibs = 0
      IF (rfnu2.LT.tol) GO TO 110
      DO 100 is=2,7
        atol = atol/rfnu2
        pp = pp*rfnu2
        IF (ias.EQ.1) GO TO 60
        sumar = zeror
        sumai = zeroi
        DO 40 k=1,kmax
          m = l1 + k
          sumar = sumar + pr(k)*alfa(m)
          sumai = sumai + pi(k)*alfa(m)
          IF (ap(k).LT.atol) GO TO 50
   40   CONTINUE
   50   CONTINUE
        asumr = asumr + sumar*pp
        asumi = asumi + sumai*pp
        IF (pp.LT.tol) ias = 1
   60   CONTINUE
        IF (ibs.EQ.1) GO TO 90
        sumbr = zeror
        sumbi = zeroi
        DO 70 k=1,kmax
          m = l2 + k
          sumbr = sumbr + pr(k)*beta(m)
          sumbi = sumbi + pi(k)*beta(m)
          IF (ap(k).LT.atol) GO TO 80
   70   CONTINUE
   80   CONTINUE
        bsumr = bsumr + sumbr*pp
        bsumi = bsumi + sumbi*pp
        IF (pp.LT.btol) ibs = 1
   90   CONTINUE
        IF (ias.EQ.1 .AND. ibs.EQ.1) GO TO 110
        l1 = l1 + 30
        l2 = l2 + 30
  100 CONTINUE
  110 CONTINUE
      asumr = asumr + coner
      pp = rfnu*rfn13
      bsumr = bsumr*pp
      bsumi = bsumi*pp
  120 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     CABS(W2).GT.0.25D0
!-----------------------------------------------------------------------
  130 CONTINUE
      CALL zsqrt(w2r, w2i, wr, wi)
      IF (wr.LT.0.0d0) wr = 0.0d0
      IF (wi.LT.0.0d0) wi = 0.0d0
      str = coner + wr
      sti = wi
      CALL zdiv(str, sti, zbr, zbi, zar, zai)
      CALL zlog(zar, zai, zcr, zci, idum)
      IF (zci.LT.0.0d0) zci = 0.0d0
      IF (zci.GT.hpi) zci = hpi
      IF (zcr.LT.0.0d0) zcr = 0.0d0
      zthr = (zcr-wr)*1.5d0
      zthi = (zci-wi)*1.5d0
      zeta1r = zcr*fnu
      zeta1i = zci*fnu
      zeta2r = wr*fnu
      zeta2i = wi*fnu
      azth = zabs(zthr,zthi)
      ang = thpi
      IF (zthr.GE.0.0d0 .AND. zthi.LT.0.0d0) GO TO 140
      ang = hpi
      IF (zthr.EQ.0.0d0) GO TO 140
      ang = datan(zthi/zthr)
      IF (zthr.LT.0.0d0) ang = ang + gpi
  140 CONTINUE
      pp = azth**ex2
      ang = ang*ex2
      zetar = pp*dcos(ang)
      zetai = pp*dsin(ang)
      IF (zetai.LT.0.0d0) zetai = 0.0d0
      argr = zetar*fn23
      argi = zetai*fn23
      CALL zdiv(zthr, zthi, zetar, zetai, rtztr, rtzti)
      CALL zdiv(rtztr, rtzti, wr, wi, zar, zai)
      tzar = zar + zar
      tzai = zai + zai
      CALL zsqrt(tzar, tzai, str, sti)
      phir = str*rfn13
      phii = sti*rfn13
      IF (ipmtr.EQ.1) GO TO 120
      raw = 1.0d0/dsqrt(aw2)
      str = wr*raw
      sti = -wi*raw
      tfnr = str*rfnu*raw
      tfni = sti*rfnu*raw
      razth = 1.0d0/azth
      str = zthr*razth
      sti = -zthi*razth
      rzthr = str*razth*rfnu
      rzthi = sti*razth*rfnu
      zcr = rzthr*ar(2)
      zci = rzthi*ar(2)
      raw2 = 1.0d0/aw2
      str = w2r*raw2
      sti = -w2i*raw2
      t2r = str*raw2
      t2i = sti*raw2
      str = t2r*c(2) + c(3)
      sti = t2i*c(2)
      upr(2) = str*tfnr - sti*tfni
      upi(2) = str*tfni + sti*tfnr
      bsumr = upr(2) + zcr
      bsumi = upi(2) + zci
      asumr = zeror
      asumi = zeroi
      IF (rfnu.LT.tol) GO TO 220
      przthr = rzthr
      przthi = rzthi
      ptfnr = tfnr
      ptfni = tfni
      upr(1) = coner
      upi(1) = conei
      pp = 1.0d0
      btol = tol*(dabs(bsumr)+dabs(bsumi))
      ks = 0
      kp1 = 2
      l = 3
      ias = 0
      ibs = 0
      DO 210 lr=2,12,2
        lrp1 = lr + 1
!-----------------------------------------------------------------------
!     COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN
!     NEXT SUMA AND SUMB
!-----------------------------------------------------------------------
        DO 160 k=lr,lrp1
          ks = ks + 1
          kp1 = kp1 + 1
          l = l + 1
          zar = c(l)
          zai = zeroi
          DO 150 j=2,kp1
            l = l + 1
            str = zar*t2r - t2i*zai + c(l)
            zai = zar*t2i + zai*t2r
            zar = str
  150     CONTINUE
          str = ptfnr*tfnr - ptfni*tfni
          ptfni = ptfnr*tfni + ptfni*tfnr
          ptfnr = str
          upr(kp1) = ptfnr*zar - ptfni*zai
          upi(kp1) = ptfni*zar + ptfnr*zai
          crr(ks) = przthr*br(ks+1)
          cri(ks) = przthi*br(ks+1)
          str = przthr*rzthr - przthi*rzthi
          przthi = przthr*rzthi + przthi*rzthr
          przthr = str
          drr(ks) = przthr*ar(ks+2)
          dri(ks) = przthi*ar(ks+2)
  160   CONTINUE
        pp = pp*rfnu2
        IF (ias.EQ.1) GO TO 180
        sumar = upr(lrp1)
        sumai = upi(lrp1)
        ju = lrp1
        DO 170 jr=1,lr
          ju = ju - 1
          sumar = sumar + crr(jr)*upr(ju) - cri(jr)*upi(ju)
          sumai = sumai + crr(jr)*upi(ju) + cri(jr)*upr(ju)
  170   CONTINUE
        asumr = asumr + sumar
        asumi = asumi + sumai
        test = dabs(sumar) + dabs(sumai)
        IF (pp.LT.tol .AND. test.LT.tol) ias = 1
  180   CONTINUE
        IF (ibs.EQ.1) GO TO 200
        sumbr = upr(lr+2) + upr(lrp1)*zcr - upi(lrp1)*zci
        sumbi = upi(lr+2) + upr(lrp1)*zci + upi(lrp1)*zcr
        ju = lrp1
        DO 190 jr=1,lr
          ju = ju - 1
          sumbr = sumbr + drr(jr)*upr(ju) - dri(jr)*upi(ju)
          sumbi = sumbi + drr(jr)*upi(ju) + dri(jr)*upr(ju)
  190   CONTINUE
        bsumr = bsumr + sumbr
        bsumi = bsumi + sumbi
        test = dabs(sumbr) + dabs(sumbi)
        IF (pp.LT.btol .AND. test.LT.btol) ibs = 1
  200   CONTINUE
        IF (ias.EQ.1 .AND. ibs.EQ.1) GO TO 220
  210 CONTINUE
  220 CONTINUE
      asumr = asumr + coner
      str = -bsumr*rfn13
      sti = -bsumi*rfn13
      CALL zdiv(str, sti, rtztr, rtzti, bsumr, bsumi)
      GO TO 120
END

SUBROUTINE zuchk(YR, YI, NZ, ASCLE, TOL)
!***BEGIN PROLOGUE  ZUCHK
!***REFER TO ZSERI,ZUOIK,ZUNK1,ZUNK2,ZUNI1,ZUNI2,ZKSCL
!
!      Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN
!      EXP(-ALIM)=ASCLE=1.0E+3*D1MACH(1)/TOL. THE TEST IS MADE TO SEE
!      IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW
!      WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED
!      IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE
!      OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE
!      ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED.
!
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZUCHK
!
!     COMPLEX Y
      DOUBLE PRECISION ascle, ss, st, tol, wr, wi, yr, yi
      INTEGER nz
      nz = 0
      wr = dabs(yr)
      wi = dabs(yi)
      st = dmin1(wr,wi)
      IF (st.GT.ascle) RETURN
      ss = dmax1(wr,wi)
      st = st/tol
      IF (ss.LT.st) nz = 1
      RETURN
END

SUBROUTINE zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, ELIM, ALIM)
USE complex
!***BEGIN PROLOGUE  ZBINU
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY

!   ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE

!***ROUTINES CALLED  ZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
!***END PROLOGUE  ZBINU
      DOUBLE PRECISION alim, az, cwi, cwr, cyi, cyr, dfnu, elim, fnu, &
      fnul, rl, tol, zeroi, zeror, zi, zr
      INTEGER i, inw, kode, n, nlast, nn, nui, nw, nz
      dimension cyr(1), cyi(1), cwr(2), cwi(2)
      DATA zeror,zeroi / 0.0d0, 0.0d0 /

      nz = 0
      az = zabs(zr,zi)
      nn = n
      dfnu = fnu + dble(float(n-1))
      IF (az.LE.2.0d0) GO TO 10
      IF (az*az*0.25d0.GT.dfnu+1.0d0) GO TO 20
   10 CONTINUE
!-----------------------------------------------------------------------
!     POWER SERIES
!-----------------------------------------------------------------------
      CALL zseri(zr, zi, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
      inw = iabs(nw)
      nz = nz + inw
      nn = nn - inw
      IF (nn.EQ.0) RETURN
      IF (nw.GE.0) GO TO 120
      dfnu = fnu + dble(float(nn-1))
   20 CONTINUE
      IF (az.LT.rl) GO TO 40
      IF (dfnu.LE.1.0d0) GO TO 30
      IF (az+az.LT.dfnu*dfnu) GO TO 50
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR LARGE Z
!-----------------------------------------------------------------------
   30 CONTINUE
      CALL zasyi(zr, zi, fnu, kode, nn, cyr, cyi, nw, rl, tol, elim, alim)
      IF (nw.LT.0) GO TO 130
      GO TO 120
   40 CONTINUE
      IF (dfnu.LE.1.0d0) GO TO 70
   50 CONTINUE
!-----------------------------------------------------------------------
!     OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
!-----------------------------------------------------------------------
      CALL zuoik(zr, zi, fnu, kode, 1, nn, cyr, cyi, nw, tol, elim, alim)
      IF (nw.LT.0) GO TO 130
      nz = nz + nw
      nn = nn - nw
      IF (nn.EQ.0) RETURN
      dfnu = fnu+dble(float(nn-1))
      IF (dfnu.GT.fnul) GO TO 110
      IF (az.GT.fnul) GO TO 110
   60 CONTINUE
      IF (az.GT.rl) GO TO 80
   70 CONTINUE
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE SERIES
!-----------------------------------------------------------------------
      CALL zmlri(zr, zi, fnu, kode, nn, cyr, cyi, nw, tol)
      IF(nw.LT.0) GO TO 130
      GO TO 120
   80 CONTINUE
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
!-----------------------------------------------------------------------
      CALL zuoik(zr, zi, fnu, kode, 2, 2, cwr, cwi, nw, tol, elim, alim)
      IF (nw.GE.0) GO TO 100
      nz = nn
      DO 90 i=1,nn
        cyr(i) = zeror
        cyi(i) = zeroi
   90 CONTINUE
      RETURN
  100 CONTINUE
      IF (nw.GT.0) GO TO 130
      CALL zwrsk(zr, zi, fnu, kode, nn, cyr, cyi, nw, cwr, cwi, tol, elim, alim)
      IF (nw.LT.0) GO TO 130
      GO TO 120
  110 CONTINUE
!-----------------------------------------------------------------------
!     INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
!-----------------------------------------------------------------------
      nui = int(sngl(fnul-dfnu)) + 1
      nui = max0(nui,0)
      CALL zbuni(zr, zi, fnu, kode, nn, cyr, cyi, nw, nui, nlast, fnul, &
      tol, elim, alim)
      IF (nw.LT.0) GO TO 130
      nz = nz + nw
      IF (nlast.EQ.0) GO TO 120
      nn = nlast
      GO TO 60
  120 CONTINUE
      RETURN
  130 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
      END

SUBROUTINE zshch(ZR, ZI, CSHR, CSHI, CCHR, CCHI)
!***BEGIN PROLOGUE  ZSHCH
!***REFER TO  ZBESK,ZBESH
!
!     ZSHCH COMPUTES THE COMPLEX HYPERBOLI! FUNCTIONS CSH=SINH(X+I*Y)
!     AND CCH=COSH(X+I*Y), WHERE I**2=-1.
!
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZSHCH
!
  DOUBLE PRECISION cchi, cchr, ch, cn, cshi, cshr, sh, sn, zi, zr, dcosh, dsinh
      sh = dsinh(zr)
      ch = dcosh(zr)
      sn = dsin(zi)
      cn = dcos(zi)
      cshr = sh*cn
      cshi = ch*sn
      cchr = ch*cn
      cchi = sh*sn
  RETURN
END

DOUBLE PRECISION FUNCTION dgamln(Z,IERR)
USE utilit
!***BEGIN PROLOGUE  DGAMLN
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5F
!***KEYWORDS  GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION
!***DESCRIPTION
!
!               **** A DOUBLE PRECISION ROUTINE ****
!         DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR
!         Z.GT.0.  THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES
!         GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION
!         G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN.  THE FUNCTION WAS MADE AS
!         PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE
!         10 DIGITS IN A WORD, RLN=AMAX1(-ALOG10(R1MACH(4)),0.5E-18)
!         LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY.
!
!         SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100
!         VALUES IS USED FOR SPEED OF EXECUTION.
!
!     DESCRIPTION OF ARGUMENTS
!
!         INPUT      Z IS D0UBLE PRECISION
!           Z      - ARGUMENT, Z.GT.0.0D0
!
!         OUTPUT      DGAMLN IS DOUBLE PRECISION
!           DGAMLN  - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0
!           IERR    - ERROR FLAG
!                     IERR=0, NORMAL RETURN, COMPUTATION COMPLETED
!                     IERR=1, Z.LE.0.0D0,    NO COMPUTATION
!
!
!***REFERENCES  COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
!***ROUTINES CALLED  I1MACH,D1MACH
!***END PROLOGUE  DGAMLN
      DOUBLE PRECISION cf, con, fln, fz, gln, rln, s, tlg, trm, tst, &
      t1, wdtol, z, zdmy, zinc, zm, zmin, zp, zsq
      INTEGER i, ierr, i1m, k, mz, nz
      dimension cf(22), gln(100)
!       LNGAMMA(N), N=1,100
      DATA gln(1), gln(2), gln(3), gln(4), gln(5), gln(6), gln(7),       &
           gln(8), gln(9), gln(10), gln(11), gln(12), gln(13), gln(14),  &
           gln(15), gln(16), gln(17), gln(18), gln(19), gln(20),         &
           gln(21), gln(22)/                                             &
           0.00000000000000000d+00,     0.00000000000000000d+00,         &
           6.93147180559945309d-01,     1.79175946922805500d+00,         &
           3.17805383034794562d+00,     4.78749174278204599d+00,         &
           6.57925121201010100d+00,     8.52516136106541430d+00,         &
           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/
      DATA gln(23), gln(24), gln(25), gln(26), gln(27), gln(28),         &
           gln(29), gln(30), gln(31), gln(32), gln(33), gln(34),         &
           gln(35), gln(36), gln(37), gln(38), gln(39), gln(40),         &
           gln(41), gln(42), gln(43), gln(44)/                           &
           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/
      DATA gln(45), gln(46), gln(47), gln(48), gln(49), gln(50),         &
           gln(51), gln(52), gln(53), gln(54), gln(55), gln(56),         &
           gln(57), gln(58), gln(59), gln(60), gln(61), gln(62),         &
           gln(63), gln(64), gln(65), gln(66)/                           &
           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/
      DATA gln(67), gln(68), gln(69), gln(70), gln(71), gln(72),         &
           gln(73), gln(74), gln(75), gln(76), gln(77), gln(78),         &
           gln(79), gln(80), gln(81), gln(82), gln(83), gln(84),         &
           gln(85), gln(86), gln(87), gln(88)/                           &
           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/
      DATA gln(89), gln(90), gln(91), gln(92), gln(93), gln(94),         &
           gln(95), gln(96), gln(97), gln(98), gln(99), gln(100)/        &
           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
      DATA cf(1), cf(2), cf(3), cf(4), cf(5), cf(6), cf(7), cf(8),       &
           cf(9), cf(10), cf(11), cf(12), cf(13), cf(14), cf(15),        &
           cf(16), cf(17), cf(18), cf(19), cf(20), cf(21), 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.39243221690590112d+00,         &
           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)
      DATA con                    /     1.83787706640934548d+00/

!***FIRST EXECUTABLE STATEMENT  DGAMLN
      ierr=0
      IF (z.LE.0.0d0) GO TO 70
      IF (z.GT.101.0d0) GO TO 10
      nz = int(sngl(z))
      fz = z - float(nz)
      IF (fz.GT.0.0d0) GO TO 10
      IF (nz.GT.100) GO TO 10
      dgamln = gln(nz)
      RETURN
   10 CONTINUE
      wdtol = d1mach(4)
      wdtol = dmax1(wdtol,0.5d-18)
      i1m = i1mach(14)
      rln = d1mach(5)*float(i1m)
      fln = dmin1(rln,20.0d0)
      fln = dmax1(fln,3.0d0)
      fln = fln - 3.0d0
      zm = 1.8000d0 + 0.3875d0*fln
      mz = int(sngl(zm)) + 1
      zmin = float(mz)
      zdmy = z
      zinc = 0.0d0
      IF (z.GE.zmin) GO TO 20
      zinc = zmin - float(nz)
      zdmy = z + zinc
   20 CONTINUE
      zp = 1.0d0/zdmy
      t1 = cf(1)*zp
      s = t1
      IF (zp.LT.wdtol) GO TO 40
      zsq = zp*zp
      tst = t1*wdtol
      DO 30 k=2,22
        zp = zp*zsq
        trm = cf(k)*zp
        IF (dabs(trm).LT.tst) GO TO 40
        s = s + trm
   30 CONTINUE
   40 CONTINUE
      IF (zinc.NE.0.0d0) GO TO 50
      tlg = dlog(z)
      dgamln = z*(tlg-1.0d0) + 0.5d0*(con-tlg) + s
      RETURN
   50 CONTINUE
      zp = 1.0d0
      nz = int(sngl(zinc))
      DO 60 i=1,nz
        zp = zp*(z+float(i-1))
   60 CONTINUE
      tlg = dlog(zdmy)
      dgamln = zdmy*(tlg-1.0d0) - dlog(zp) + 0.5d0*(con-tlg) + s
      RETURN

   70 CONTINUE
      ierr=1
      RETURN
END

SUBROUTINE zkscl(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
USE complex
!***BEGIN PROLOGUE  ZKSCL
!***REFER TO  ZBESK
!
!     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
!     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
!     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
!
!***ROUTINES CALLED  ZUCHK,ZABS,ZLOG
!***END PROLOGUE  ZKSCL
!     COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
      DOUBLE PRECISION acs, as, ascle, cki, ckr, csi, csr, cyi, &
      cyr, elim, fn, fnu, rzi, rzr, str, s1i, s1r, s2i, s2r,    &
      tol, yi, yr, zeroi, zeror, zri, zrr, zdr, zdi, celmr,     &
      elm, helim, alas
      INTEGER i, ic, idum, kk, n, nn, nw, nz
      dimension yr(1), yi(1), cyr(2), cyi(2)
      DATA zeror,zeroi / 0.0d0 , 0.0d0 /

      nz = 0
      ic = 0
      nn = min0(2,n)
      DO 10 i=1,nn
        s1r = yr(i)
        s1i = yi(i)
        cyr(i) = s1r
        cyi(i) = s1i
        as = zabs(s1r,s1i)
        acs = -zrr + dlog(as)
        nz = nz + 1
        yr(i) = zeror
        yi(i) = zeroi
        IF (acs.LT.(-elim)) GO TO 10
        CALL zlog(s1r, s1i, csr, csi, idum)
        csr = csr - zrr
        csi = csi - zri
        str = dexp(csr)/tol
        csr = str*dcos(csi)
        csi = str*dsin(csi)
        CALL zuchk(csr, csi, nw, ascle, tol)
        IF (nw.NE.0) GO TO 10
        yr(i) = csr
        yi(i) = csi
        ic = i
        nz = nz - 1
   10 CONTINUE
      IF (n.EQ.1) RETURN
      IF (ic.GT.1) GO TO 20
      yr(1) = zeror
      yi(1) = zeroi
      nz = 2
   20 CONTINUE
      IF (n.EQ.2) RETURN
      IF (nz.EQ.0) RETURN
      fn = fnu + 1.0d0
      ckr = fn*rzr
      cki = fn*rzi
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      helim = 0.5d0*elim
      elm = dexp(-elim)
      celmr = elm
      zdr = zrr
      zdi = zri

!     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
!     S2 GETS LARGER THAN EXP(ELIM/2)

      DO 30 i=3,n
        kk = i
        csr = s2r
        csi = s2i
        s2r = ckr*csr - cki*csi + s1r
        s2i = cki*csr + ckr*csi + s1i
        s1r = csr
        s1i = csi
        ckr = ckr + rzr
        cki = cki + rzi
        as = zabs(s2r,s2i)
        alas = dlog(as)
        acs = -zdr + alas
        nz = nz + 1
        yr(i) = zeror
        yi(i) = zeroi
        IF (acs.LT.(-elim)) GO TO 25
        CALL zlog(s2r, s2i, csr, csi, idum)
        csr = csr - zdr
        csi = csi - zdi
        str = dexp(csr)/tol
        csr = str*dcos(csi)
        csi = str*dsin(csi)
        CALL zuchk(csr, csi, nw, ascle, tol)
        IF (nw.NE.0) GO TO 25
        yr(i) = csr
        yi(i) = csi
        nz = nz - 1
        IF (ic.EQ.kk-1) GO TO 40
        ic = kk
        GO TO 30
   25   CONTINUE
        IF(alas.LT.helim) GO TO 30
        zdr = zdr - elim
        s1r = s1r*celmr
        s1i = s1i*celmr
        s2r = s2r*celmr
        s2i = s2i*celmr
   30 CONTINUE
      nz = n
      IF(ic.EQ.n) nz=n-1
      GO TO 45
   40 CONTINUE
      nz = kk - 2
   45 CONTINUE
      DO 50 i=1,nz
        yr(i) = zeror
        yi(i) = zeroi
   50 CONTINUE
      RETURN
END

SUBROUTINE zacai(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZACAI
!***REFER TO  ZAIRY
!
!     ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
!
!         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
!                 MP=PI*MR*CMPLX(0.0,1.0)
!
!     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
!     HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
!     ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
!     RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
!     IS CALLED FROM ZAIRY.
!
!***ROUTINES CALLED  ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS
!***END PROLOGUE  ZACAI
!     COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
      DOUBLE PRECISION alim, arg, ascle, az, csgnr, csgni, cspnr,     &
      cspni, c1r, c1i, c2r, c2i, cyr, cyi, dfnu, elim, fmr, fnu, pi,  &
      rl, sgn, tol, yy, yr, yi, zr, zi, znr, zni
      INTEGER inu, iuf, kode, mr, n, nn, nw, nz
      dimension yr(1), yi(1), cyr(2), cyi(2)
      DATA pi / 3.14159265358979324d0 /
      nz = 0
      znr = -zr
      zni = -zi
      az = zabs(zr,zi)
      nn = n
      dfnu = fnu + dble(float(n-1))
      IF (az.LE.2.0d0) GO TO 10
      IF (az*az*0.25d0.GT.dfnu+1.0d0) GO TO 20
   10 CONTINUE
!-----------------------------------------------------------------------
!     POWER SERIES FOR THE I FUNCTION
!-----------------------------------------------------------------------
      CALL zseri(znr, zni, fnu, kode, nn, yr, yi, nw, tol, elim, alim)
      GO TO 40
   20 CONTINUE
      IF (az.LT.rl) GO TO 30
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
!-----------------------------------------------------------------------
      CALL zasyi(znr, zni, fnu, kode, nn, yr, yi, nw, rl, tol, elim, alim)
      IF (nw.LT.0) GO TO 80
      GO TO 40
   30 CONTINUE
!-----------------------------------------------------------------------
!     MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
!-----------------------------------------------------------------------
      CALL zmlri(znr, zni, fnu, kode, nn, yr, yi, nw, tol)
      IF(nw.LT.0) GO TO 80
   40 CONTINUE
!-----------------------------------------------------------------------
!     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!-----------------------------------------------------------------------
      CALL zbknu(znr, zni, fnu, kode, 1, cyr, cyi, nw, tol, elim, alim)
      IF (nw.NE.0) GO TO 80
      fmr = dble(float(mr))
      sgn = -dsign(pi,fmr)
      csgnr = 0.0d0
      csgni = sgn
      IF (kode.EQ.1) GO TO 50
      yy = -zni
      csgnr = -csgni*dsin(yy)
      csgni = csgni*dcos(yy)
   50 CONTINUE
!-----------------------------------------------------------------------
!     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
      inu = int(sngl(fnu))
      arg = (fnu-dble(float(inu)))*sgn
      cspnr = dcos(arg)
      cspni = dsin(arg)
      IF (mod(inu,2).EQ.0) GO TO 60
      cspnr = -cspnr
      cspni = -cspni
   60 CONTINUE
      c1r = cyr(1)
      c1i = cyi(1)
      c2r = yr(1)
      c2i = yi(1)
      IF (kode.EQ.1) GO TO 70
      iuf = 0
      ascle = 1.0d+3*d1mach(1)/tol
      CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
      nz = nz + nw
   70 CONTINUE
      yr(1) = cspnr*c1r - cspni*c1i + csgnr*c2r - csgni*c2i
      yi(1) = cspnr*c1i + cspni*c1r + csgnr*c2i + csgni*c2r
      RETURN
   80 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
END

SUBROUTINE zs1s2(ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF)
USE complex
!***BEGIN PROLOGUE  ZS1S2
!***REFER TO  ZBESK,ZAIRY
!
!     ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE
!     ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON-
!     TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION.
!     ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF
!     MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER
!     OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE
!     PRECISION ABOVE THE UNDERFLOW LIMIT.
!
!***ROUTINES CALLED  ZABS,ZEXP,ZLOG
!***END PROLOGUE  ZS1S2
!     COMPLEX CZERO,C1,S1,S1D,S2,ZR
      DOUBLE PRECISION aa, alim, aln, ascle, as1, as2, c1i, c1r, s1di, &
      s1dr, s1i, s1r, s2i, s2r, zeroi, zeror, zri, zrr
      INTEGER iuf, idum, nz
      DATA zeror,zeroi  / 0.0d0 , 0.0d0 /
      nz = 0
      as1 = zabs(s1r,s1i)
      as2 = zabs(s2r,s2i)
      IF (s1r.EQ.0.0d0 .AND. s1i.EQ.0.0d0) GO TO 10
      IF (as1.EQ.0.0d0) GO TO 10
      aln = -zrr - zrr + dlog(as1)
      s1dr = s1r
      s1di = s1i
      s1r = zeror
      s1i = zeroi
      as1 = zeror
      IF (aln.LT.(-alim)) GO TO 10
      CALL zlog(s1dr, s1di, c1r, c1i, idum)
      c1r = c1r - zrr - zrr
      c1i = c1i - zri - zri
      CALL zexp(c1r, c1i, s1r, s1i)
      as1 = zabs(s1r,s1i)
      iuf = iuf + 1
   10 CONTINUE
      aa = dmax1(as1,as2)
      IF (aa.GT.ascle) RETURN
      s1r = zeror
      s1i = zeroi
      s2r = zeror
      s2i = zeroi
      nz = 1
      iuf = 0
      RETURN
END

SUBROUTINE zrati(ZR, ZI, FNU, N, CYR, CYI, TOL)
USE complex
!***BEGIN PROLOGUE  ZRATI
!***REFER TO  ZBESI,ZBESK,ZBESH
!
!     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
!     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
!     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
!     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
!     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
!     BY D. J. SOOKNE.
!
!***ROUTINES CALLED  ZABS,ZDIV
!***END PROLOGUE  ZRATI
!     COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU
      DOUBLE PRECISION ak, amagz, ap1, ap2, arg, az, cdfnui, cdfnur, &
      conei, coner, cyi, cyr, czeroi, czeror, dfnu, fdnu, flam, fnu, &
      fnup, pti, ptr, p1i, p1r, p2i, p2r, rak, rap1, rho, rt2, rzi,  &
      rzr, test, test1, tol, tti, ttr, t1i, t1r, zi, zr
      INTEGER i, id, idnu, inu, itime, k, kk, magz, n
      dimension cyr(1), cyi(1)
      DATA czeror,czeroi,coner,conei,rt2  &
      /0.0d0, 0.0d0, 1.0d0, 0.0d0, 1.41421356237309505d0/
      az = zabs(zr,zi)
      inu = int(sngl(fnu))
      idnu = inu + n - 1
      magz = int(sngl(az))
      amagz = dble(float(magz+1))
      fdnu = dble(float(idnu))
      fnup = dmax1(amagz,fdnu)
      id = idnu - magz - 1
      itime = 1
      k = 1
      ptr = 1.0d0/az
      rzr = ptr*(zr+zr)*ptr
      rzi = -ptr*(zi+zi)*ptr
      t1r = rzr*fnup
      t1i = rzi*fnup
      p2r = -t1r
      p2i = -t1i
      p1r = coner
      p1i = conei
      t1r = t1r + rzr
      t1i = t1i + rzi
      IF (id.GT.0) id = 0
      ap2 = zabs(p2r,p2i)
      ap1 = zabs(p1r,p1i)
!-----------------------------------------------------------------------
!     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
!     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 = dsqrt(arg)
      test = test1
      rap1 = 1.0d0/ap1
      p1r = p1r*rap1
      p1i = p1i*rap1
      p2r = p2r*rap1
      p2i = p2i*rap1
      ap2 = ap2*rap1
   10 CONTINUE
      k = k + 1
      ap1 = ap2
      ptr = p2r
      pti = p2i
      p2r = p1r - (t1r*ptr-t1i*pti)
      p2i = p1i - (t1r*pti+t1i*ptr)
      p1r = ptr
      p1i = pti
      t1r = t1r + rzr
      t1i = t1i + rzi
      ap2 = zabs(p2r,p2i)
      IF (ap1.LE.test) GO TO 10
      IF (itime.EQ.2) GO TO 20
      ak = zabs(t1r,t1i)*0.5d0
      flam = ak + dsqrt(ak*ak-1.0d0)
      rho = dmin1(ap2/ap1,flam)
      test = test1*dsqrt(rho/(rho*rho-1.0d0))
      itime = 2
      GO TO 10
   20 CONTINUE
      kk = k + 1 - id
      ak = dble(float(kk))
      t1r = ak
      t1i = czeroi
      dfnu = fnu + dble(float(n-1))
      p1r = 1.0d0/ap2
      p1i = czeroi
      p2r = czeror
      p2i = czeroi
      DO 30 i=1,kk
        ptr = p1r
        pti = p1i
        rap1 = dfnu + t1r
        ttr = rzr*rap1
        tti = rzi*rap1
        p1r = (ptr*ttr-pti*tti) + p2r
        p1i = (ptr*tti+pti*ttr) + p2i
        p2r = ptr
        p2i = pti
        t1r = t1r - coner
   30 CONTINUE
      IF (p1r.NE.czeror .OR. p1i.NE.czeroi) GO TO 40
      p1r = tol
      p1i = tol
   40 CONTINUE
      CALL zdiv(p2r, p2i, p1r, p1i, cyr(n), cyi(n))
      IF (n.EQ.1) RETURN
      k = n - 1
      ak = dble(float(k))
      t1r = ak
      t1i = czeroi
      cdfnur = fnu*rzr
      cdfnui = fnu*rzi
      DO 60 i=2,n
        ptr = cdfnur + (t1r*rzr-t1i*rzi) + cyr(k+1)
        pti = cdfnui + (t1r*rzi+t1i*rzr) + cyi(k+1)
        ak = zabs(ptr,pti)
        IF (ak.NE.czeror) GO TO 50
        ptr = tol
        pti = tol
        ak = tol*rt2
   50   CONTINUE
        rak = coner/ak
        cyr(k) = rak*ptr*rak
        cyi(k) = -rak*pti*rak
        t1r = t1r - coner
        k = k - 1
   60 CONTINUE
      RETURN
END

SUBROUTINE zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZAIRY
!***DATE WRITTEN   830501   (YYMMDD)
!***REVISION DATE  830501   (YYMMDD)
!***CATEGORY NO.  B5K
!***KEYWORDS  AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
!***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
!***PURPOSE  TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z
!***DESCRIPTION
!
!                      ***A DOUBLE PRECISION ROUTINE***
!         ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR
!         ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
!         KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)*
!         DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
!         -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN
!         PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z).
!
!         WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTI! IN
!         THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED
!         FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS.
!         DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
!         MATHEMATICAL FUNCTIONS (REF. 1).
!
!         INPUT      ZR,ZI ARE DOUBLE PRECISION
!           ZR,ZI  - Z=CMPLX(ZR,ZI)
!           ID     - ORDER OF DERIVATIVE, ID=0 OR ID=1
!           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
!                    KODE= 1  RETURNS
!                             AI=AI(Z)                ON ID=0 OR
!                             AI=DAI(Z)/DZ            ON ID=1
!                        = 2  RETURNS
!                             AI=CEXP(ZTA)*AI(Z)       ON ID=0 OR
!                             AI=CEXP(ZTA)*DAI(Z)/DZ   ON ID=1 WHERE
!                             ZTA=(2/3)*Z*CSQRT(Z)
!
!         OUTPUT     AIR,AII ARE DOUBLE PRECISION
!           AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
!                    KODE
!           NZ     - UNDERFLOW INDICATOR
!                    NZ= 0   , NORMAL RETURN
!                    NZ= 1   , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN
!                              -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
!           IERR   - ERROR FLAG
!                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
!                    IERR=1, INPUT ERROR   - NO COMPUTATION
!                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(ZTA)
!                            TOO LARGE ON KODE=1
!                    IERR=3, CABS(Z) LARGE      - COMPUTATION COMPLETED
!                            LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
!                            PRODUCE LESS THAN HALF OF MACHINE ACCURACY
!                    IERR=4, CABS(Z) TOO LARGE  - NO COMPUTATION
!                            COMPLETE LOSS OF ACCURACY BY ARGUMENT
!                            REDUCTION
!                    IERR=5, ERROR              - NO COMPUTATION,
!                            ALGORITHM TERMINATION CONDITION NOT MET
!
!***LONG DESCRIPTION
!
!         AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL
!         FUNCTIONS BY
!
!            AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA)
!                           C=1.0/(PI*SQRT(3.0))
!                            ZTA=(2/3)*Z**(3/2)
!
!         WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
!
!         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
!         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
!         OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
!         THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
!         THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
!         FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
!         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
!         ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
!         ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
!         FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
!         LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
!         MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
!         AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
!         PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
!         PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
!         ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
!         NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
!         DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
!         EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
!         NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
!         PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
!         MACHINES.
!
!         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
!         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
!         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
!         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
!         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
!         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
!         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
!         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
!         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
!         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
!         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
!         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
!         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
!         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
!         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
!         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
!         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
!         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
!         OR -PI/2+P.
!
!***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
!                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
!                 COMMERCE, 1955.
!
!               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
!                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
!
!               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
!                 1018, MAY, 1985
!
!               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
!                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
!                 MATH. SOFTWARE, 1986
!
!***ROUTINES CALLED  ZACAI,ZBKNU,ZABS,ZEXP,ZSQRT,I1MACH,D1MACH
!***END PROLOGUE  ZAIRY
!     COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3
      DOUBLE PRECISION aa, ad, aii, air, ak, alim, atrm, az, az3, bk,  &
      cc, ck, coef, conei, coner, csqi, csqr, cyi, cyr, c1, c2, dig,   &
      dk, d1, d2, elim, fid, fnu, ptr, rl, r1m5, sfac, sti, str,       &
      s1i, s1r, s2i, s2r, tol, trm1i, trm1r, trm2i, trm2r, tth, zeroi, &
      zeror, zi, zr, ztai, ztar, z3i, z3r, alaz, bb
      INTEGER id, ierr, iflag, k, kode, k1, k2, mr, nn, nz
      dimension cyr(1), cyi(1)
      DATA tth, c1, c2, coef /6.66666666666666667d-01, &
      3.55028053887817240d-01,2.58819403792806799d-01, &
      1.83776298473930683d-01/
      DATA zeror, zeroi, coner, conei /0.0d0,0.0d0,1.0d0,0.0d0/
!***FIRST EXECUTABLE STATEMENT  ZAIRY
      ierr = 0
      nz=0
      IF (id.LT.0 .OR. id.GT.1) ierr=1
      IF (kode.LT.1 .OR. kode.GT.2) ierr=1
      IF (ierr.NE.0) RETURN
      az = zabs(zr,zi)
      tol = dmax1(d1mach(4),1.0d-18)
      fid = dble(float(id))
      IF (az.GT.1.0d0) GO TO 70
!-----------------------------------------------------------------------
!     POWER SERIES FOR CABS(Z).LE.1.
!-----------------------------------------------------------------------
      s1r = coner
      s1i = conei
      s2r = coner
      s2i = conei
      IF (az.LT.tol) GO TO 170
      aa = az*az
      IF (aa.LT.tol/az) GO TO 40
      trm1r = coner
      trm1i = conei
      trm2r = coner
      trm2i = conei
      atrm = 1.0d0
      str = zr*zr - zi*zi
      sti = zr*zi + zi*zr
      z3r = str*zr - sti*zi
      z3i = str*zi + sti*zr
      az3 = az*aa
      ak = 2.0d0 + fid
      bk = 3.0d0 - fid - fid
      ck = 4.0d0 - fid
      dk = 3.0d0 + fid + fid
      d1 = ak*dk
      d2 = bk*ck
      ad = dmin1(d1,d2)
      ak = 24.0d0 + 9.0d0*fid
      bk = 30.0d0 - 9.0d0*fid
      DO 30 k=1,25
        str = (trm1r*z3r-trm1i*z3i)/d1
        trm1i = (trm1r*z3i+trm1i*z3r)/d1
        trm1r = str
        s1r = s1r + trm1r
        s1i = s1i + trm1i
        str = (trm2r*z3r-trm2i*z3i)/d2
        trm2i = (trm2r*z3i+trm2i*z3r)/d2
        trm2r = str
        s2r = s2r + trm2r
        s2i = s2i + trm2i
        atrm = atrm*az3/ad
        d1 = d1 + ak
        d2 = d2 + bk
        ad = dmin1(d1,d2)
        IF (atrm.LT.tol*ad) GO TO 40
        ak = ak + 18.0d0
        bk = bk + 18.0d0
   30 CONTINUE
   40 CONTINUE
      IF (id.EQ.1) GO TO 50
      air = s1r*c1 - c2*(zr*s2r-zi*s2i)
      aii = s1i*c1 - c2*(zr*s2i+zi*s2r)
      IF (kode.EQ.1) RETURN
      CALL zsqrt(zr, zi, str, sti)
      ztar = tth*(zr*str-zi*sti)
      ztai = tth*(zr*sti+zi*str)
      CALL zexp(ztar, ztai, str, sti)
      ptr = air*str - aii*sti
      aii = air*sti + aii*str
      air = ptr
      RETURN
   50 CONTINUE
      air = -s2r*c2
      aii = -s2i*c2
      IF (az.LE.tol) GO TO 60
      str = zr*s1r - zi*s1i
      sti = zr*s1i + zi*s1r
      cc = c1/(1.0d0+fid)
      air = air + cc*(str*zr-sti*zi)
      aii = aii + cc*(str*zi+sti*zr)
   60 CONTINUE
      IF (kode.EQ.1) RETURN
      CALL zsqrt(zr, zi, str, sti)
      ztar = tth*(zr*str-zi*sti)
      ztai = tth*(zr*sti+zi*str)
      CALL zexp(ztar, ztai, str, sti)
      ptr = str*air - sti*aii
      aii = str*aii + sti*air
      air = ptr
      RETURN
!-----------------------------------------------------------------------
!     CASE FOR CABS(Z).GT.1.0
!-----------------------------------------------------------------------
   70 CONTINUE
      fnu = (1.0d0+fid)/3.0d0
!-----------------------------------------------------------------------
!     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
!     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18.
!     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
!     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
!     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
!     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
!     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
!     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
!-----------------------------------------------------------------------
      k1 = i1mach(15)
      k2 = i1mach(16)
      r1m5 = d1mach(5)
      k = min0(iabs(k1),iabs(k2))
      elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
      k1 = i1mach(14) - 1
      aa = r1m5*dble(float(k1))
      dig = dmin1(aa,18.0d0)
      aa = aa*2.303d0
      alim = elim + dmax1(-aa,-41.45d0)
      rl = 1.2d0*dig + 3.0d0
      alaz = dlog(az)
!-----------------------------------------------------------------------
!     TEST FOR PROPER RANGE
!-----------------------------------------------------------------------
      aa=0.5d0/tol
      bb=dble(float(i1mach(9)))*0.5d0
      aa=dmin1(aa,bb)
      aa=aa**tth
      IF (az.GT.aa) GO TO 260
      aa=dsqrt(aa)
      IF (az.GT.aa) ierr=3
      CALL zsqrt(zr, zi, csqr, csqi)
      ztar = tth*(zr*csqr-zi*csqi)
      ztai = tth*(zr*csqi+zi*csqr)
!-----------------------------------------------------------------------
!     RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
!-----------------------------------------------------------------------
      iflag = 0
      sfac = 1.0d0
      ak = ztai
      IF (zr.GE.0.0d0) GO TO 80
      bk = ztar
      ck = -dabs(bk)
      ztar = ck
      ztai = ak
   80 CONTINUE
      IF (zi.NE.0.0d0) GO TO 90
      IF (zr.GT.0.0d0) GO TO 90
      ztar = 0.0d0
      ztai = ak
   90 CONTINUE
      aa = ztar
      IF (aa.GE.0.0d0 .AND. zr.GT.0.0d0) GO TO 110
      IF (kode.EQ.2) GO TO 100
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
      IF (aa.GT.(-alim)) GO TO 100
      aa = -aa + 0.25d0*alaz
      iflag = 1
      sfac = tol
      IF (aa.GT.elim) GO TO 270
  100 CONTINUE
!-----------------------------------------------------------------------
!     CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
!-----------------------------------------------------------------------
      mr = 1
      IF (zi.LT.0.0d0) mr = -1
      CALL zacai(ztar, ztai, fnu, kode, mr, 1, cyr, cyi, nn, rl, tol, &
      elim, alim)
      IF (nn.LT.0) GO TO 280
      nz = nz + nn
      GO TO 130
  110 CONTINUE
      IF (kode.EQ.2) GO TO 120
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
      IF (aa.LT.alim) GO TO 120
      aa = -aa - 0.25d0*alaz
      iflag = 2
      sfac = 1.0d0/tol
      IF (aa.LT.(-elim)) GO TO 210
  120 CONTINUE
      CALL zbknu(ztar, ztai, fnu, kode, 1, cyr, cyi, nz, tol, elim, alim)
  130 CONTINUE
      s1r = cyr(1)*coef
      s1i = cyi(1)*coef
      IF (iflag.NE.0) GO TO 150
      IF (id.EQ.1) GO TO 140
      air = csqr*s1r - csqi*s1i
      aii = csqr*s1i + csqi*s1r
      RETURN
  140 CONTINUE
      air = -(zr*s1r-zi*s1i)
      aii = -(zr*s1i+zi*s1r)
      RETURN
  150 CONTINUE
      s1r = s1r*sfac
      s1i = s1i*sfac
      IF (id.EQ.1) GO TO 160
      str = s1r*csqr - s1i*csqi
      s1i = s1r*csqi + s1i*csqr
      s1r = str
      air = s1r/sfac
      aii = s1i/sfac
      RETURN
  160 CONTINUE
      str = -(s1r*zr-s1i*zi)
      s1i = -(s1r*zi+s1i*zr)
      s1r = str
      air = s1r/sfac
      aii = s1i/sfac
      RETURN
  170 CONTINUE
      aa = 1.0d+3*d1mach(1)
      s1r = zeror
      s1i = zeroi
      IF (id.EQ.1) GO TO 190
      IF (az.LE.aa) GO TO 180
      s1r = c2*zr
      s1i = c2*zi
  180 CONTINUE
      air = c1 - s1r
      aii = -s1i
      RETURN
  190 CONTINUE
      air = -c2
      aii = 0.0d0
      aa = dsqrt(aa)
      IF (az.LE.aa) GO TO 200
      s1r = 0.5d0*(zr*zr-zi*zi)
      s1i = zr*zi
  200 CONTINUE
      air = air + c1*s1r
      aii = aii + c1*s1i
      RETURN
  210 CONTINUE
      nz = 1
      air = zeror
      aii = zeroi
      RETURN
  270 CONTINUE
      nz = 0
      ierr=2
      RETURN
  280 CONTINUE
      IF(nn.EQ.(-1)) GO TO 270
      nz=0
      ierr=5
      RETURN
  260 CONTINUE
      ierr=4
      nz=0
      RETURN
END

SUBROUTINE zseri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZSERI
!***REFER TO  ZBESI,ZBESK
!
!     ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
!     MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
!     REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
!     NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
!     DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
!     CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
!     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
!
!***ROUTINES CALLED  DGAMLN,D1MACH,ZUCHK,ZABS,ZDIV,ZLOG,ZMLT
!***END PROLOGUE  ZSERI
!     COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z
      DOUBLE PRECISION aa, acz, ak, ak1i, ak1r, alim, arm, ascle, atol, &
      az, cki, ckr, coefi, coefr, conei, coner, crscr, czi, czr, dfnu, &
      elim, fnu, fnup, hzi, hzr, raz, rs, rtr1, rzi, rzr, s, ss, sti,  &
      str, s1i, s1r, s2i, s2r, tol, yi, yr, wi, wr, zeroi, zeror, zi,  &
      zr!, DGAMLN
      INTEGER i, ib, idum, iflag, il, k, kode, l, m, n, nn, nz, nw
      dimension yr(1), yi(1), wr(2), wi(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /

      nz = 0
      az = zabs(zr,zi)
      IF (az.EQ.0.0d0) GO TO 160
      arm = 1.0d+3*d1mach(1)
      rtr1 = dsqrt(arm)
      crscr = 1.0d0
      iflag = 0
      IF (az.LT.arm) GO TO 150
      hzr = 0.5d0*zr
      hzi = 0.5d0*zi
      czr = zeror
      czi = zeroi
      IF (az.LE.rtr1) GO TO 10
      CALL zmlt(hzr, hzi, hzr, hzi, czr, czi)
   10 CONTINUE
      acz = zabs(czr,czi)
      nn = n
      CALL zlog(hzr, hzi, ckr, cki, idum)
   20 CONTINUE
      dfnu = fnu + dble(float(nn-1))
      fnup = dfnu + 1.0d0
!-----------------------------------------------------------------------
!     UNDERFLOW TEST
!-----------------------------------------------------------------------
      ak1r = ckr*dfnu
      ak1i = cki*dfnu
      ak = dgamln(fnup,idum)
      ak1r = ak1r - ak
      IF (kode.EQ.2) ak1r = ak1r - zr
      IF (ak1r.GT.(-elim)) GO TO 40
   30 CONTINUE
      nz = nz + 1
      yr(nn) = zeror
      yi(nn) = zeroi
      IF (acz.GT.dfnu) GO TO 190
      nn = nn - 1
      IF (nn.EQ.0) RETURN
      GO TO 20
   40 CONTINUE
      IF (ak1r.GT.(-alim)) GO TO 50
      iflag = 1
      ss = 1.0d0/tol
      crscr = tol
      ascle = arm*ss
   50 CONTINUE
      aa = dexp(ak1r)
      IF (iflag.EQ.1) aa = aa*ss
      coefr = aa*dcos(ak1i)
      coefi = aa*dsin(ak1i)
      atol = tol*acz/fnup
      il = min0(2,nn)
      DO 90 i=1,il
        dfnu = fnu + dble(float(nn-i))
        fnup = dfnu + 1.0d0
        s1r = coner
        s1i = conei
        IF (acz.LT.tol*fnup) GO TO 70
        ak1r = coner
        ak1i = conei
        ak = fnup + 2.0d0
        s = fnup
        aa = 2.0d0
   60   CONTINUE
        rs = 1.0d0/s
        str = ak1r*czr - ak1i*czi
        sti = ak1r*czi + ak1i*czr
        ak1r = str*rs
        ak1i = sti*rs
        s1r = s1r + ak1r
        s1i = s1i + ak1i
        s = s + ak
        ak = ak + 2.0d0
        aa = aa*acz*rs
        IF (aa.GT.atol) GO TO 60
   70   CONTINUE
        s2r = s1r*coefr - s1i*coefi
        s2i = s1r*coefi + s1i*coefr
        wr(i) = s2r
        wi(i) = s2i
        IF (iflag.EQ.0) GO TO 80
        CALL zuchk(s2r, s2i, nw, ascle, tol)
        IF (nw.NE.0) GO TO 30
   80   CONTINUE
        m = nn - i + 1
        yr(m) = s2r*crscr
        yi(m) = s2i*crscr
        IF (i.EQ.il) GO TO 90
        CALL zdiv(coefr, coefi, hzr, hzi, str, sti)
        coefr = str*dfnu
        coefi = sti*dfnu
   90 CONTINUE
      IF (nn.LE.2) RETURN
      k = nn - 2
      ak = dble(float(k))
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      IF (iflag.EQ.1) GO TO 120
      ib = 3
  100 CONTINUE
      DO 110 i=ib,nn
        yr(k) = (ak+fnu)*(rzr*yr(k+1)-rzi*yi(k+1)) + yr(k+2)
        yi(k) = (ak+fnu)*(rzr*yi(k+1)+rzi*yr(k+1)) + yi(k+2)
        ak = ak - 1.0d0
        k = k - 1
  110 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     RECUR BACKWARD WITH SCALED VALUES
!-----------------------------------------------------------------------
  120 CONTINUE
!-----------------------------------------------------------------------
!     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
!     UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3
!-----------------------------------------------------------------------
      s1r = wr(1)
      s1i = wi(1)
      s2r = wr(2)
      s2i = wi(2)
      DO 130 l=3,nn
        ckr = s2r
        cki = s2i
        s2r = s1r + (ak+fnu)*(rzr*ckr-rzi*cki)
        s2i = s1i + (ak+fnu)*(rzr*cki+rzi*ckr)
        s1r = ckr
        s1i = cki
        ckr = s2r*crscr
        cki = s2i*crscr
        yr(k) = ckr
        yi(k) = cki
        ak = ak - 1.0d0
        k = k - 1
        IF (zabs(ckr,cki).GT.ascle) GO TO 140
  130 CONTINUE
      RETURN
  140 CONTINUE
      ib = l + 1
      IF (ib.GT.nn) RETURN
      GO TO 100
  150 CONTINUE
      nz = n
      IF (fnu.EQ.0.0d0) nz = nz - 1
  160 CONTINUE
      yr(1) = zeror
      yi(1) = zeroi
      IF (fnu.NE.0.0d0) GO TO 170
      yr(1) = coner
      yi(1) = conei
  170 CONTINUE
      IF (n.EQ.1) RETURN
      DO 180 i=2,n
        yr(i) = zeror
        yi(i) = zeroi
  180 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
!     THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
!-----------------------------------------------------------------------
  190 CONTINUE
      nz = -nz
      RETURN
      END

SUBROUTINE zasyi(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZASYI
!***REFER TO  ZBESI,ZBESK

!     ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
!     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
!     REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
!     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
!
!***ROUTINES CALLED  D1MACH,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT
!***END PROLOGUE  ZASYI
!     COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
      DOUBLE PRECISION aa, aez, ak, ak1i, ak1r, alim, arg, arm, atol,   &
      az, bb, bk, cki, ckr, conei, coner, cs1i, cs1r, cs2i, cs2r, czi, &
      czr, dfnu, dki, dkr, dnu2, elim, ezi, ezr, fdn, fnu, pi, p1i,    &
      p1r, raz, rl, rtpi, rtr1, rzi, rzr, s, sgn, sqk, sti, str, s2i,  &
      s2r, tol, tzi, tzr, yi, yr, zeroi, zeror, zi, zr
      INTEGER i, ib, il, inu, j, jl, k, kode, koded, m, n, nn, nz
      dimension yr(1), yi(1)
      DATA pi, rtpi  /3.14159265358979324d0 , 0.159154943091895336d0 /
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /

      nz = 0
      az = zabs(zr,zi)
      arm = 1.0d+3*d1mach(1)
      rtr1 = dsqrt(arm)
      il = min0(2,n)
      dfnu = fnu + dble(float(n-il))
!-----------------------------------------------------------------------
!     OVERFLOW TEST
!-----------------------------------------------------------------------
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      ak1r = rtpi*str*raz
      ak1i = rtpi*sti*raz
      CALL zsqrt(ak1r, ak1i, ak1r, ak1i)
      czr = zr
      czi = zi
      IF (kode.NE.2) GO TO 10
      czr = zeror
      czi = zi
   10 CONTINUE
      IF (dabs(czr).GT.elim) GO TO 100
      dnu2 = dfnu + dfnu
      koded = 1
      IF ((dabs(czr).GT.alim) .AND. (n.GT.2)) GO TO 20
      koded = 0
      CALL zexp(czr, czi, str, sti)
      CALL zmlt(ak1r, ak1i, str, sti, ak1r, ak1i)
   20 CONTINUE
      fdn = 0.0d0
      IF (dnu2.GT.rtr1) fdn = dnu2*dnu2
      ezr = zr*8.0d0
      ezi = zi*8.0d0
!-----------------------------------------------------------------------
!     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.0d0*az
      s = tol/aez
      jl = int(sngl(rl+rl)) + 2
      p1r = zeror
      p1i = zeroi
      IF (zi.EQ.0.0d0) GO TO 30
!-----------------------------------------------------------------------
!     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
!     SIGNIFICANCE WHEN FNU OR N IS LARGE
!-----------------------------------------------------------------------
      inu = int(sngl(fnu))
      arg = (fnu-dble(float(inu)))*pi
      inu = inu + n - il
      ak = -dsin(arg)
      bk = dcos(arg)
      IF (zi.LT.0.0d0) bk = -bk
      p1r = ak
      p1i = bk
      IF (mod(inu,2).EQ.0) GO TO 30
      p1r = -p1r
      p1i = -p1i
   30 CONTINUE
      DO 70 k=1,il
        sqk = fdn - 1.0d0
        atol = s*dabs(sqk)
        sgn = 1.0d0
        cs1r = coner
        cs1i = conei
        cs2r = coner
        cs2i = conei
        ckr = coner
        cki = conei
        ak = 0.0d0
        aa = 1.0d0
        bb = aez
        dkr = ezr
        dki = ezi
        DO 40 j=1,jl
          CALL zdiv(ckr, cki, dkr, dki, str, sti)
          ckr = str*sqk
          cki = sti*sqk
          cs2r = cs2r + ckr
          cs2i = cs2i + cki
          sgn = -sgn
          cs1r = cs1r + ckr*sgn
          cs1i = cs1i + cki*sgn
          dkr = dkr + ezr
          dki = dki + ezi
          aa = aa*dabs(sqk)/bb
          bb = bb + aez
          ak = ak + 8.0d0
          sqk = sqk - ak
          IF (aa.LE.atol) GO TO 50
   40   CONTINUE
        GO TO 110
   50   CONTINUE
        s2r = cs1r
        s2i = cs1i
        IF (zr+zr.GE.elim) GO TO 60
        tzr = zr + zr
        tzi = zi + zi
        CALL zexp(-tzr, -tzi, str, sti)
        CALL zmlt(str, sti, p1r, p1i, str, sti)
        CALL zmlt(str, sti, cs2r, cs2i, str, sti)
        s2r = s2r + str
        s2i = s2i + sti
   60   CONTINUE
        fdn = fdn + 8.0d0*dfnu + 4.0d0
        p1r = -p1r
        p1i = -p1i
        m = n - il + k
        yr(m) = s2r*ak1r - s2i*ak1i
        yi(m) = s2r*ak1i + s2i*ak1r
   70 CONTINUE
      IF (n.LE.2) RETURN
      nn = n
      k = nn - 2
      ak = dble(float(k))
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      ib = 3
      DO 80 i=ib,nn
        yr(k) = (ak+fnu)*(rzr*yr(k+1)-rzi*yi(k+1)) + yr(k+2)
        yi(k) = (ak+fnu)*(rzr*yi(k+1)+rzi*yr(k+1)) + yi(k+2)
        ak = ak - 1.0d0
        k = k - 1
   80 CONTINUE
      IF (koded.EQ.0) RETURN
      CALL zexp(czr, czi, ckr, cki)
      DO 90 i=1,nn
        str = yr(i)*ckr - yi(i)*cki
        yi(i) = yr(i)*cki + yi(i)*ckr
        yr(i) = str
   90 CONTINUE
      RETURN
  100 CONTINUE
      nz = -1
      RETURN
  110 CONTINUE
      nz=-2
      RETURN
      END

SUBROUTINE zbuni(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST, FNUL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZBUNI
!***REFER TO  ZBESI,ZBESK
!
!     ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
!     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
!     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
!     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
!     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
!
!***ROUTINES CALLED  ZUNI1,ZUNI2,ZABS,D1MACH
!***END PROLOGUE  ZBUNI
!     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z
      DOUBLE PRECISION alim, ax, ay, csclr, cscrr, cyi, cyr, dfnu,    &
      elim, fnu, fnui, fnul, gnu, raz, rzi, rzr, sti, str, s1i, s1r,  &
      s2i, s2r, tol, yi, yr, zi, zr, ascle, bry, c1r, c1i, c1m
      INTEGER i, iflag, iform, k, kode, n, nl, nlast, nui, nw, nz
      dimension yr(1), yi(1), cyr(2), cyi(2), bry(3)
      nz = 0
      ax = dabs(zr)*1.7321d0
      ay = dabs(zi)
      iform = 1
      IF (ay.GT.ax) iform = 2
      IF (nui.EQ.0) GO TO 60
      fnui = dble(float(nui))
      dfnu = fnu + dble(float(n-1))
      gnu = dfnu + fnui
      IF (iform.EQ.2) GO TO 10
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3.LE.ARG(Z).LE.PI/3
!-----------------------------------------------------------------------
      CALL zuni1(zr, zi, gnu, kode, 2, cyr, cyi, nw, nlast, fnul, tol, elim, alim)
      GO TO 20
   10 CONTINUE
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!-----------------------------------------------------------------------
      CALL zuni2(zr, zi, gnu, kode, 2, cyr, cyi, nw, nlast, fnul, tol, elim, alim)
   20 CONTINUE
      IF (nw.LT.0) GO TO 50
      IF (nw.NE.0) GO TO 90
      str = zabs(cyr(1),cyi(1))
!----------------------------------------------------------------------
!     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
!----------------------------------------------------------------------
      bry(1)=1.0d+3*d1mach(1)/tol
      bry(2) = 1.0d0/bry(1)
      bry(3) = bry(2)
      iflag = 2
      ascle = bry(2)
      csclr = 1.0d0
      IF (str.GT.bry(1)) GO TO 21
      iflag = 1
      ascle = bry(1)
      csclr = 1.0d0/tol
      GO TO 25
   21 CONTINUE
      IF (str.LT.bry(2)) GO TO 25
      iflag = 3
      ascle=bry(3)
      csclr = tol
   25 CONTINUE
      cscrr = 1.0d0/csclr
      s1r = cyr(2)*csclr
      s1i = cyi(2)*csclr
      s2r = cyr(1)*csclr
      s2i = cyi(1)*csclr
      raz = 1.0d0/zabs(zr,zi)
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      DO 30 i=1,nui
        str = s2r
        sti = s2i
        s2r = (dfnu+fnui)*(rzr*str-rzi*sti) + s1r
        s2i = (dfnu+fnui)*(rzr*sti+rzi*str) + s1i
        s1r = str
        s1i = sti
        fnui = fnui - 1.0d0
        IF (iflag.GE.3) GO TO 30
        str = s2r*cscrr
        sti = s2i*cscrr
        c1r = dabs(str)
        c1i = dabs(sti)
        c1m = dmax1(c1r,c1i)
        IF (c1m.LE.ascle) GO TO 30
        iflag = iflag+1
        ascle = bry(iflag)
        s1r = s1r*cscrr
        s1i = s1i*cscrr
        s2r = str
        s2i = sti
        csclr = csclr*tol
        cscrr = 1.0d0/csclr
        s1r = s1r*csclr
        s1i = s1i*csclr
        s2r = s2r*csclr
        s2i = s2i*csclr
   30 CONTINUE
      yr(n) = s2r*cscrr
      yi(n) = s2i*cscrr
      IF (n.EQ.1) RETURN
      nl = n - 1
      fnui = dble(float(nl))
      k = nl
      DO 40 i=1,nl
        str = s2r
        sti = s2i
        s2r = (fnu+fnui)*(rzr*str-rzi*sti) + s1r
        s2i = (fnu+fnui)*(rzr*sti+rzi*str) + s1i
        s1r = str
        s1i = sti
        str = s2r*cscrr
        sti = s2i*cscrr
        yr(k) = str
        yi(k) = sti
        fnui = fnui - 1.0d0
        k = k - 1
        IF (iflag.GE.3) GO TO 40
        c1r = dabs(str)
        c1i = dabs(sti)
        c1m = dmax1(c1r,c1i)
        IF (c1m.LE.ascle) GO TO 40
        iflag = iflag+1
        ascle = bry(iflag)
        s1r = s1r*cscrr
        s1i = s1i*cscrr
        s2r = str
        s2i = sti
        csclr = csclr*tol
        cscrr = 1.0d0/csclr
        s1r = s1r*csclr
        s1i = s1i*csclr
        s2r = s2r*csclr
        s2i = s2i*csclr
   40 CONTINUE
      RETURN
   50 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
   60 CONTINUE
      IF (iform.EQ.2) GO TO 70
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
!     -PI/3.LE.ARG(Z).LE.PI/3
!-----------------------------------------------------------------------
      CALL zuni1(zr, zi, fnu, kode, n, yr, yi, nw, nlast, fnul, tol, elim, alim)
      GO TO 80
   70 CONTINUE
!-----------------------------------------------------------------------
!     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
!     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
!     AND HPI=PI/2
!-----------------------------------------------------------------------
      CALL zuni2(zr, zi, fnu, kode, n, yr, yi, nw, nlast, fnul, tol, elim, alim)
   80 CONTINUE
      IF (nw.LT.0) GO TO 50
      nz = nw
      RETURN
   90 CONTINUE
      nlast = n
      RETURN
      END

SUBROUTINE zuni1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZUNI1
!***REFER TO  ZBESI,ZBESK
!
!     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
!     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
!
!     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
!     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
!     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
!     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
!     Y(I)=CZERO FOR I=NLAST+1,N
!
!***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,ZABS
!***END PROLOGUE  ZUNI1
!     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
!    *S2,Y,Z,ZETA1,ZETA2
      DOUBLE PRECISION alim, aphi, ascle, bry, conei, coner, crsc,  &
      cscl, csrr, cssr, cwrki, cwrkr, c1r, c2i, c2m, c2r, elim, fn, &
      fnu, fnul, phii, phir, rast, rs1, rzi, rzr, sti, str, sumi,   &
      sumr, s1i, s1r, s2i, s2r, tol, yi, yr, zeroi, zeror, zeta1i,  &
      zeta1r, zeta2i, zeta2r, zi, zr, cyr, cyi
      INTEGER i, iflag, init, k, kode, m, n, nd, nlast, nn, nuf, nw, nz
      dimension bry(3), yr(1), yi(1), cwrkr(16), cwrki(16), cssr(3), &
      csrr(3), cyr(2), cyi(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /

      nz = 0
      nd = n
      nlast = 0
!-----------------------------------------------------------------------
!     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
!     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
!     EXP(ALIM)=EXP(ELIM)*TOL
!-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      crsc = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = crsc
      csrr(1) = crsc
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = 1.0d+3*d1mach(1)/tol
!-----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!-----------------------------------------------------------------------
      fn = dmax1(fnu,1.0d0)
      init = 0
      CALL zunik(zr, zi, fn, 1, 1, tol, init, phir, phii, zeta1r, &
      zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
      IF (kode.EQ.1) GO TO 10
      str = zr + zeta2r
      sti = zi + zeta2i
      rast = fn/zabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = -zeta1r + str
      s1i = -zeta1i + sti
      GO TO 20
   10 CONTINUE
      s1r = -zeta1r + zeta2r
      s1i = -zeta1i + zeta2i
   20 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 130
   30 CONTINUE
      nn = min0(2,nd)
      DO 80 i=1,nn
        fn = fnu + dble(float(nd-i))
        init = 0
        CALL zunik(zr, zi, fn, 1, 0, tol, init, phir, phii, zeta1r, &
        zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
        IF (kode.EQ.1) GO TO 40
        str = zr + zeta2r
        sti = zi + zeta2i
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zeta1r + str
        s1i = -zeta1i + sti + zi
        GO TO 50
   40   CONTINUE
        s1r = -zeta1r + zeta2r
        s1i = -zeta1i + zeta2i
   50   CONTINUE
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 110
        IF (i.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 60
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phir,phii)
        rs1 = rs1 + dlog(aphi)
        IF (dabs(rs1).GT.elim) GO TO 110
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 60
        IF (i.EQ.1) iflag = 3
   60   CONTINUE
!-----------------------------------------------------------------------
!     SCALE S1 IF CABS(S1).LT.ASCLE
!-----------------------------------------------------------------------
        s2r = phir*sumr - phii*sumi
        s2i = phir*sumi + phii*sumr
        str = dexp(s1r)*cssr(iflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s2r*s1i + s2i*s1r
        s2r = str
        IF (iflag.NE.1) GO TO 70
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 110
   70   CONTINUE
        cyr(i) = s2r
        cyi(i) = s2i
        m = nd - i + 1
        yr(m) = s2r*csrr(iflag)
        yi(m) = s2i*csrr(iflag)
   80 CONTINUE
      IF (nd.LE.2) GO TO 100
      rast = 1.0d0/zabs(zr,zi)
      str = zr*rast
      sti = -zi*rast
      rzr = (str+str)*rast
      rzi = (sti+sti)*rast
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(iflag)
      ascle = bry(iflag)
      k = nd - 2
      fn = dble(float(k))
      DO 90 i=3,nd
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(k) = c2r
        yi(k) = c2i
        k = k - 1
        fn = fn - 1.0d0
        IF (iflag.GE.3) GO TO 90
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 90
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        c1r = csrr(iflag)
   90 CONTINUE
  100 CONTINUE
      RETURN
!-----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!-----------------------------------------------------------------------
  110 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 120
      yr(nd) = zeror
      yi(nd) = zeroi
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) GO TO 100
      CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
      IF (nuf.LT.0) GO TO 120
      nd = nd - nuf
      nz = nz + nuf
      IF (nd.EQ.0) GO TO 100
      fn = fnu + dble(float(nd-1))
      IF (fn.GE.fnul) GO TO 30
      nlast = nd
      RETURN
  120 CONTINUE
      nz = -1
      RETURN
  130 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 120
      nz = n
      DO 140 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  140 CONTINUE
      RETURN
      END

SUBROUTINE zuni2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZUNI2
!***REFER TO  ZBESI,ZBESK
!
!     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
!     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
!     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
!
!     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
!     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
!     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
!     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
!     Y(I)=CZERO FOR I=NLAST+1,N
!
!***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,ZABS
!***END PROLOGUE  ZUNI2
!     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
!    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
      DOUBLE PRECISION aarg, aic, aii, air, alim, ang, aphi, argi,     &
      argr, ascle, asumi, asumr, bry, bsumi, bsumr, cidi, cipi, cipr,  &
      conei, coner, crsc, cscl, csrr, cssr, c1r, c2i, c2m, c2r, daii,  &
      dair, elim, fn, fnu, fnul, hpi, phii, phir, rast, raz, rs1, rzi, &
      rzr, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, zbi, zbr, zeroi, &
      zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi, zni, znr, zr, cyr, cyi
      INTEGER i, iflag, in, inu, j, k, kode, n, nai, nd, ndai, nlast,  &
      nn, nuf, nw, nz, idum
      dimension bry(3), yr(1), yi(1), cipr(4), cipi(4), cssr(3),       &
      csrr(3), cyr(2), cyi(2)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
      DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),    &
      cipi(4)/ 1.0d0,0.0d0, 0.0d0,1.0d0, -1.0d0,0.0d0, 0.0d0,-1.0d0/
      DATA hpi, aic  /1.57079632679489662d+00,  1.265512123484645396d+00/

      nz = 0
      nd = n
      nlast = 0
!-----------------------------------------------------------------------
!     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
!     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
!     EXP(ALIM)=EXP(ELIM)*TOL
!-----------------------------------------------------------------------
      cscl = 1.0d0/tol
      crsc = tol
      cssr(1) = cscl
      cssr(2) = coner
      cssr(3) = crsc
      csrr(1) = crsc
      csrr(2) = coner
      csrr(3) = cscl
      bry(1) = 1.0d+3*d1mach(1)/tol
!-----------------------------------------------------------------------
!     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
!-----------------------------------------------------------------------
      znr = zi
      zni = -zr
      zbr = zr
      zbi = zi
      cidi = -coner
      inu = int(sngl(fnu))
      ang = hpi*(fnu-dble(float(inu)))
      c2r = dcos(ang)
      c2i = dsin(ang)
      in = inu + n - 1
      in = mod(in,4) + 1
      str = c2r*cipr(in) - c2i*cipi(in)
      c2i = c2r*cipi(in) + c2i*cipr(in)
      c2r = str
      IF (zi.GT.0.0d0) GO TO 10
      znr = -znr
      zbi = -zbi
      cidi = -cidi
      c2i = -c2i
   10 CONTINUE
!-----------------------------------------------------------------------
!     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
!-----------------------------------------------------------------------
      fn = dmax1(fnu,1.0d0)
      CALL zunhj(znr, zni, fn, 1, tol, phir, phii, argr, argi, zeta1r, &
      zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
      IF (kode.EQ.1) GO TO 20
      str = zbr + zeta2r
      sti = zbi + zeta2i
      rast = fn/zabs(str,sti)
      str = str*rast*rast
      sti = -sti*rast*rast
      s1r = -zeta1r + str
      s1i = -zeta1i + sti
      GO TO 30
   20 CONTINUE
      s1r = -zeta1r + zeta2r
      s1i = -zeta1i + zeta2i
   30 CONTINUE
      rs1 = s1r
      IF (dabs(rs1).GT.elim) GO TO 150
   40 CONTINUE
      nn = min0(2,nd)
      DO 90 i=1,nn
        fn = fnu + dble(float(nd-i))
        CALL zunhj(znr, zni, fn, 0, tol, phir, phii, argr, argi, &
        zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
        IF (kode.EQ.1) GO TO 50
        str = zbr + zeta2r
        sti = zbi + zeta2i
        rast = fn/zabs(str,sti)
        str = str*rast*rast
        sti = -sti*rast*rast
        s1r = -zeta1r + str
        s1i = -zeta1i + sti + dabs(zi)
        GO TO 60
   50   CONTINUE
        s1r = -zeta1r + zeta2r
        s1i = -zeta1i + zeta2i
   60   CONTINUE
!-----------------------------------------------------------------------
!     TEST FOR UNDERFLOW AND OVERFLOW
!-----------------------------------------------------------------------
        rs1 = s1r
        IF (dabs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 2
        IF (dabs(rs1).LT.alim) GO TO 70
!-----------------------------------------------------------------------
!     REFINE  TEST AND SCALE
!-----------------------------------------------------------------------
        aphi = zabs(phir,phii)
        aarg = zabs(argr,argi)
        rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
        IF (dabs(rs1).GT.elim) GO TO 120
        IF (i.EQ.1) iflag = 1
        IF (rs1.LT.0.0d0) GO TO 70
        IF (i.EQ.1) iflag = 3
   70   CONTINUE
!-----------------------------------------------------------------------
!     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
!     EXPONENT EXTREMES
!-----------------------------------------------------------------------
        CALL zairy(argr, argi, 0, 2, air, aii, nai, idum)
        CALL zairy(argr, argi, 1, 2, dair, daii, ndai, idum)
        str = dair*bsumr - daii*bsumi
        sti = dair*bsumi + daii*bsumr
        str = str + (air*asumr-aii*asumi)
        sti = sti + (air*asumi+aii*asumr)
        s2r = phir*str - phii*sti
        s2i = phir*sti + phii*str
        str = dexp(s1r)*cssr(iflag)
        s1r = str*dcos(s1i)
        s1i = str*dsin(s1i)
        str = s2r*s1r - s2i*s1i
        s2i = s2r*s1i + s2i*s1r
        s2r = str
        IF (iflag.NE.1) GO TO 80
        CALL zuchk(s2r, s2i, nw, bry(1), tol)
        IF (nw.NE.0) GO TO 120
   80   CONTINUE
        IF (zi.LE.0.0d0) s2i = -s2i
        str = s2r*c2r - s2i*c2i
        s2i = s2r*c2i + s2i*c2r
        s2r = str
        cyr(i) = s2r
        cyi(i) = s2i
        j = nd - i + 1
        yr(j) = s2r*csrr(iflag)
        yi(j) = s2i*csrr(iflag)
        str = -c2i*cidi
        c2i = c2r*cidi
        c2r = str
   90 CONTINUE
      IF (nd.LE.2) GO TO 110
      raz = 1.0d0/zabs(zr,zi)
      str = zr*raz
      sti = -zi*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      bry(2) = 1.0d0/bry(1)
      bry(3) = d1mach(2)
      s1r = cyr(1)
      s1i = cyi(1)
      s2r = cyr(2)
      s2i = cyi(2)
      c1r = csrr(iflag)
      ascle = bry(iflag)
      k = nd - 2
      fn = dble(float(k))
      DO 100 i=3,nd
        c2r = s2r
        c2i = s2i
        s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
        s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
        s1r = c2r
        s1i = c2i
        c2r = s2r*c1r
        c2i = s2i*c1r
        yr(k) = c2r
        yi(k) = c2i
        k = k - 1
        fn = fn - 1.0d0
        IF (iflag.GE.3) GO TO 100
        str = dabs(c2r)
        sti = dabs(c2i)
        c2m = dmax1(str,sti)
        IF (c2m.LE.ascle) GO TO 100
        iflag = iflag + 1
        ascle = bry(iflag)
        s1r = s1r*c1r
        s1i = s1i*c1r
        s2r = c2r
        s2i = c2i
        s1r = s1r*cssr(iflag)
        s1i = s1i*cssr(iflag)
        s2r = s2r*cssr(iflag)
        s2i = s2i*cssr(iflag)
        c1r = csrr(iflag)
  100 CONTINUE
  110 CONTINUE
      RETURN
  120 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 140
!-----------------------------------------------------------------------
!     SET UNDERFLOW AND UPDATE PARAMETERS
!-----------------------------------------------------------------------
      yr(nd) = zeror
      yi(nd) = zeroi
      nz = nz + 1
      nd = nd - 1
      IF (nd.EQ.0) GO TO 110
      CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
      IF (nuf.LT.0) GO TO 140
      nd = nd - nuf
      nz = nz + nuf
      IF (nd.EQ.0) GO TO 110
      fn = fnu + dble(float(nd-1))
      IF (fn.LT.fnul) GO TO 130
      fn = cidi
      j = nuf + 1
      k = mod(j,4) + 1
      s1r = cipr(k)
      s1i = cipi(k)
      IF (fn.LT.0.0d0) s1i = -s1i
      str = c2r*s1r - c2i*s1i
      c2i = c2r*s1i + c2i*s1r
      c2r = str
      GO TO 40
  130 CONTINUE
      nlast = nd
      RETURN
  140 CONTINUE
      nz = -1
      RETURN
  150 CONTINUE
      IF (rs1.GT.0.0d0) GO TO 140
      nz = n
      DO 160 i=1,n
        yr(i) = zeror
        yi(i) = zeroi
  160 CONTINUE
      RETURN
      END

SUBROUTINE zwrsk(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI, TOL, ELIM, ALIM)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZWRSK
!***REFER TO  ZBESI,ZBESK
!
!     ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
!     NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
!
!***ROUTINES CALLED  D1MACH,ZBKNU,ZRATI,ZABS
!***END PROLOGUE  ZWRSK
!     COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
      DOUBLE PRECISION act, acw, alim, ascle, cinui, cinur, csclr, cti, &
      ctr, cwi, cwr, c1i, c1r, c2i, c2r, elim, fnu, pti, ptr, ract,     &
      sti, str, tol, yi, yr, zri, zrr
      INTEGER i, kode, n, nw, nz
      dimension yr(1), yi(1), cwr(2), cwi(2)
!-----------------------------------------------------------------------
!     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 zbknu(zrr, zri, fnu, kode, 2, cwr, cwi, nw, tol, elim, alim)
      IF (nw.NE.0) GO TO 50
      CALL zrati(zrr, zri, fnu, n, yr, yi, tol)
!-----------------------------------------------------------------------
!     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
!     R(FNU+J-1,Z)=Y(J),  J=1,...,N
!-----------------------------------------------------------------------
      cinur = 1.0d0
      cinui = 0.0d0
      IF (kode.EQ.1) GO TO 10
      cinur = dcos(zri)
      cinui = dsin(zri)
   10 CONTINUE
!-----------------------------------------------------------------------
!     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 = zabs(cwr(2),cwi(2))
      ascle = 1.0d+3*d1mach(1)/tol
      csclr = 1.0d0
      IF (acw.GT.ascle) GO TO 20
      csclr = 1.0d0/tol
      GO TO 30
   20 CONTINUE
      ascle = 1.0d0/ascle
      IF (acw.LT.ascle) GO TO 30
      csclr = tol
   30 CONTINUE
      c1r = cwr(1)*csclr
      c1i = cwi(1)*csclr
      c2r = cwr(2)*csclr
      c2i = cwi(2)*csclr
      str = yr(1)
      sti = yi(1)
!-----------------------------------------------------------------------
!     CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS
!     UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
!-----------------------------------------------------------------------
      ptr = str*c1r - sti*c1i
      pti = str*c1i + sti*c1r
      ptr = ptr + c2r
      pti = pti + c2i
      ctr = zrr*ptr - zri*pti
      cti = zrr*pti + zri*ptr
      act = zabs(ctr,cti)
      ract = 1.0d0/act
      ctr = ctr*ract
      cti = -cti*ract
      ptr = cinur*ract
      pti = cinui*ract
      cinur = ptr*ctr - pti*cti
      cinui = ptr*cti + pti*ctr
      yr(1) = cinur*csclr
      yi(1) = cinui*csclr
      IF (n.EQ.1) RETURN
      DO 40 i=2,n
        ptr = str*cinur - sti*cinui
        cinui = str*cinui + sti*cinur
        cinur = ptr
        str = yr(i)
        sti = yi(i)
        yr(i) = cinur*csclr
        yi(i) = cinui*csclr
   40 CONTINUE
      RETURN
   50 CONTINUE
      nz = -1
      IF(nw.EQ.(-2)) nz=-2
      RETURN
END

SUBROUTINE zmlri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
USE utilit
USE complex
!***BEGIN PROLOGUE  ZMLRI
!***REFER TO  ZBESI,ZBESK
!
!     ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
!     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
!
!***ROUTINES CALLED  DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT
!***END PROLOGUE  ZMLRI
!     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
      DOUBLE PRECISION ack, ak, ap, at, az, bk, cki, ckr, cnormi,     &
      cnormr, conei, coner, fkap, fkk, flam, fnf, fnu, pti, ptr, p1i, &
      p1r, p2i, p2r, raz, rho, rho2, rzi, rzr, scle, sti, str, sumi,  &
      sumr, tfnf, tol, tst, yi, yr, zeroi, zeror, zi, zr!, DGAMLN
      INTEGER i, iaz, idum, ifnu, inu, itime, k, kk, km, kode, m, n, nz
      dimension yr(1), yi(1)
      DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
      scle = d1mach(1)/tol
      nz=0
      az = zabs(zr,zi)
      iaz = int(sngl(az))
      ifnu = int(sngl(fnu))
      inu = ifnu + n - 1
      at = dble(float(iaz)) + 1.0d0
      raz = 1.0d0/az
      str = zr*raz
      sti = -zi*raz
      ckr = str*at*raz
      cki = sti*at*raz
      rzr = (str+str)*raz
      rzi = (sti+sti)*raz
      p1r = zeror
      p1i = zeroi
      p2r = coner
      p2i = conei
      ack = (at+1.0d0)*raz
      rho = ack + dsqrt(ack*ack-1.0d0)
      rho2 = rho*rho
      tst = (rho2+rho2)/((rho2-1.0d0)*(rho-1.0d0))
      tst = tst/tol
!-----------------------------------------------------------------------
!     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
!-----------------------------------------------------------------------
      ak = at
      DO 10 i=1,80
        ptr = p2r
        pti = p2i
        p2r = p1r - (ckr*ptr-cki*pti)
        p2i = p1i - (cki*ptr+ckr*pti)
        p1r = ptr
        p1i = pti
        ckr = ckr + rzr
        cki = cki + rzi
        ap = zabs(p2r,p2i)
        IF (ap.GT.tst*ak*ak) GO TO 20
        ak = ak + 1.0d0
   10 CONTINUE
      GO TO 110
   20 CONTINUE
      i = i + 1
      k = 0
      IF (inu.LT.iaz) GO TO 40
!-----------------------------------------------------------------------
!     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
!-----------------------------------------------------------------------
      p1r = zeror
      p1i = zeroi
      p2r = coner
      p2i = conei
      at = dble(float(inu)) + 1.0d0
      str = zr*raz
      sti = -zi*raz
      ckr = str*at*raz
      cki = sti*at*raz
      ack = at*raz
      tst = dsqrt(ack/tol)
      itime = 1
      DO 30 k=1,80
        ptr = p2r
        pti = p2i
        p2r = p1r - (ckr*ptr-cki*pti)
        p2i = p1i - (ckr*pti+cki*ptr)
        p1r = ptr
        p1i = pti
        ckr = ckr + rzr
        cki = cki + rzi
        ap = zabs(p2r,p2i)
        IF (ap.LT.tst) GO TO 30
        IF (itime.EQ.2) GO TO 40
        ack = zabs(ckr,cki)
        flam = ack + dsqrt(ack*ack-1.0d0)
        fkap = ap/zabs(p1r,p1i)
        rho = dmin1(flam,fkap)
        tst = tst*dsqrt(rho/(rho*rho-1.0d0))
        itime = 2
   30 CONTINUE
      GO TO 110
   40 CONTINUE
!-----------------------------------------------------------------------
!     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
!-----------------------------------------------------------------------
      k = k + 1
      kk = max0(i+iaz,k+inu)
      fkk = dble(float(kk))
      p1r = zeror
      p1i = zeroi
!-----------------------------------------------------------------------
!     SCALE P2 AND SUM BY SCLE
!-----------------------------------------------------------------------
      p2r = scle
      p2i = zeroi
      fnf = fnu - dble(float(ifnu))
      tfnf = fnf + fnf
      bk = dgamln(fkk+tfnf+1.0d0,idum) - dgamln(fkk+1.0d0,idum) - &
      dgamln(tfnf+1.0d0,idum)
      bk = dexp(bk)
      sumr = zeror
      sumi = zeroi
      km = kk - inu
      DO 50 i=1,km
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
   50 CONTINUE
      yr(n) = p2r
      yi(n) = p2i
      IF (n.EQ.1) GO TO 70
      DO 60 i=2,n
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
        m = n - i + 1
        yr(m) = p2r
        yi(m) = p2i
   60 CONTINUE
   70 CONTINUE
      IF (ifnu.LE.0) GO TO 90
      DO 80 i=1,ifnu
        ptr = p2r
        pti = p2i
        p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
        p2i = p1i + (fkk+fnf)*(rzr*pti+rzi*ptr)
        p1r = ptr
        p1i = pti
        ak = 1.0d0 - tfnf/(fkk+tfnf)
        ack = bk*ak
        sumr = sumr + (ack+bk)*p1r
        sumi = sumi + (ack+bk)*p1i
        bk = ack
        fkk = fkk - 1.0d0
   80 CONTINUE
   90 CONTINUE
      ptr = zr
      pti = zi
      IF (kode.EQ.2) ptr = zeror
      CALL zlog(rzr, rzi, str, sti, idum)
      p1r = -fnf*str + ptr
      p1i = -fnf*sti + pti
      ap = dgamln(1.0d0+fnf,idum)
      ptr = p1r - ap
      pti = p1i
!-----------------------------------------------------------------------
!     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
!     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
!-----------------------------------------------------------------------
      p2r = p2r + sumr
      p2i = p2i + sumi
      ap = zabs(p2r,p2i)
      p1r = 1.0d0/ap
      CALL zexp(ptr, pti, str, sti)
      ckr = str*p1r
      cki = sti*p1r
      ptr = p2r*p1r
      pti = -p2i*p1r
      CALL zmlt(ckr, cki, ptr, pti, cnormr, cnormi)
      DO 100 i=1,n
        str = yr(i)*cnormr - yi(i)*cnormi
        yi(i) = yr(i)*cnormi + yi(i)*cnormr
        yr(i) = str
  100 CONTINUE
      RETURN
  110 CONTINUE
      nz=-2
      RETURN
END
end module bessely
!end of file tzbesy.f90