besselj.f90

module besselj
    implicit none
    private
    public zbesj
contains
!*****************************************************************
!* EVALUATE A J-BESSEL FUNCTION OF COMPLEX ARGUMENT (FIRST KIND) *
!* ------------------------------------------------------------- *
!* SAMPLE RUN:                                                   *
!* (Evaluate J0 to J4 for argument Z=(1.0,2.0) ).                *
!*                                                               *
!* zr(0) =   1.586259                                            *
!* zi(0) =  -1.391602                                            *
!* zr(1) =   1.291848                                            *
!* zi(1) =   1.010488                                            *
!* zr(2) =  -0.261130                                            *
!* zi(2) =   0.762320                                            *
!* zr(3) =  -0.281040                                            *
!* zi(3) =   0.017175                                            *
!* zr(4) =  -0.034898                                            *
!* zi(4) =  -0.067215                                            *
!* 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)             *
!*****************************************************************
!Note: to link with Complex.f90 and Utilit.f90.
!---------------------------------------------
!PROGRAM TEST_ZBESJ
!real*8 zr,zi, cyr(10), cyi(10)

!  n=5
!  zr=1.d0; zi=2.d0
!  call ZBESJ(zr,zi,0,1,n,cyr,cyi,nz,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 *,' '

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

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

!***FIRST EXECUTABLE STATEMENT  ZBESJ
      ierr = 0
      nz=0
      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
!-----------------------------------------------------------------------
!     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)
      rl = 1.2d0*dig + 3.0d0
      fnul = 10.0d0 + 6.0d0*(dig-3.0d0)
!-----------------------------------------------------------------------
!     TEST FOR PROPER RANGE
!-----------------------------------------------------------------------
      az = zabs(zr,zi)
      fn = fnu+dble(float(n-1))
      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
!-----------------------------------------------------------------------
!     CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
!     WHEN FNU IS LARGE
!-----------------------------------------------------------------------
      cii = 1.0d0
      inu = int(sngl(fnu))
      inuh = inu/2
      ir = inu - 2*inuh
      arg = (fnu-dble(float(inu-ir)))*hpi
      csgnr = dcos(arg)
      csgni = dsin(arg)
      IF (mod(inuh,2).EQ.0) GO TO 40
      csgnr = -csgnr
      csgni = -csgni
   40 CONTINUE
!-----------------------------------------------------------------------
!     ZN IS IN THE RIGHT HALF PLANE
!-----------------------------------------------------------------------
      znr = zi
      zni = -zr
      IF (zi.GE.0.0d0) GO TO 50
      znr = -znr
      zni = -zni
      csgni = -csgni
      cii = -cii
   50 CONTINUE
      CALL zbinu(znr, zni, fnu, kode, n, cyr, cyi, nz, rl, fnul, tol, &
      elim, alim)
      IF (nz.LT.0) GO TO 130
      nl = n - nz
      IF (nl.EQ.0) RETURN
      DO 60 i=1,nl
        str = cyr(i)*csgnr - cyi(i)*csgni
        cyi(i) = cyr(i)*csgni + cyi(i)*csgnr
        cyr(i) = str
        str = -csgni*cii
        csgni = csgnr*cii
        csgnr = str
   60 CONTINUE
      RETURN
  130 CONTINUE
      IF(nz.EQ.(-2)) GO TO 140
      nz = 0
      ierr = 2
      RETURN
  140 CONTINUE
      nz=0
      ierr=5
      RETURN
  260 CONTINUE
      nz=0
      ierr=4
      RETURN
      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 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 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 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 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 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 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

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

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 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 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 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,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 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
end module besselj
!end of file tzbesj.f90