C Subroutine FNHANK: This computes Hankel functions of the first kind
C  and of order zero and one.

C Input parameter
C complex Z    : The argument of the Hankel function (Im(Z)>=0).

C Output parameter
C H(0:1)       : The Hankel functions of order zero and one respectively.
C
      SUBROUTINE FNHANK(Z,H)
      COMPLEX*16 Z,H(0:1)
      COMPLEX  ZR, CY(2)
      INTEGER NZ(2), IERR
      ZR = Z
      CALL HANKEL01(ZR, CY, NZ, IERR)
      IF(IERR.EQ.0) THEN
        H(0)=CY(1)
        H(1)=CY(2)
      ELSE
        WRITE(*,*) ' Error in FNHANK for argument ', Z
        STOP
      END IF
      RETURN
      END
C******       FNU=0.     KODE=1       M=1       N=2      ******
C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
C                    NZ= 0   , NORMAL RETURN
C                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO
C                              DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0)
C                              J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
C                              Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
C                              HALF PLANES, NZ STATES ONLY THE NUMBER
C                              OF UNDERFLOWS.
C           IERR    -ERROR FLAG
C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
C                    IERR=1, INPUT ERROR   - NO COMPUTATION
C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU+N-1 TOO
C                            LARGE OR CABS(Z) TOO SMALL OR BOTH
C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
C                            ACCURACY
C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
C                            CANCE BY ARGUMENT REDUCTION
C                    IERR=5, ERROR              - NO COMPUTATION,
C                            ALGORITHM TERMINATION CONDITION NOT MET
C                    IERR=5, ERROR              - NO COMPUTATION,
C                            aimag(Z) < 0
      SUBROUTINE HANKEL01(Z, CY, NZ, IERR)
      COMPLEX CY, Z, ZN, ZT, CSGN
      REAL AA, ALIM, ALN, ARG, AZ, CPN, DIG, ELIM, FMM, FN, FNUL,
     * HPI, RHPI, RL, R1M5, SGN, SPN, TOL, UFL, XN, XX, YN, YY, R1MACH,
     * BB, ASCLE, RTOL, ATOL
      INTEGER I, IERR, INU, INUH, IR, K, K1, K2, M,
     * MM, MR, NN, NUF, NW, NZ, I1MACH
      DIMENSION CY(2)
      DATA HPI /1.57079632679489662E0/
      NZ=0
      XX = REAL(Z)
      if ( aimag(Z) .LE. -1.0E-7) then
          IERR=7
      else
          YY = AIMAG(Z)
          YY = AMAX1(0.0, YY)
      end if
      IF ( Z .EQ.(0.0,0.0) ) THEN
          IERR=1
          RETURN
      ELSE
          IERR = 0
      END IF
      NN = 2
      TOL = AMAX1(R1MACH(4),1.0E-18)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN0(IABS(K1),IABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      K1 = I1MACH(11) - 1
      AA = R1M5*FLOAT(K1)
      DIG = AMIN1(AA,18.0E0)
      AA = AA*2.303E0
      ALIM = ELIM + AMAX1(-AA,-41.45E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      FN = 1.0
      MM = 1
      FMM = FLOAT(MM)
      ZN = Z*CMPLX(0.0E0,-1.0)
      XN = REAL(ZN)
      YN = AIMAG(ZN)
      AZ = CABS(Z)
      AA = 0.5E0/TOL
      BB=FLOAT(I1MACH(9))*0.5E0
      AA=AMIN1(AA,BB)
      IF(AZ.GT.AA) THEN
         NZ=0
         IERR=4
         RETURN
      END IF
      AA=SQRT(AA)
      IF(AZ.GT.AA ) IERR=3
      UFL = R1MACH(1)*1.0E+3
      IF (AZ.LT.UFL) THEN
         IERR=2
         NZ=0
         RETURN
      END IF
      CALL CBKNU(ZN, NN, CY, NZ, TOL, ELIM, ALIM)
      SGN = SIGN(HPI,-FMM)
      INU = 0
      INUH = 0
      IR = 0
      ARG = 0.0
      RHPI = 1.0E0/SGN
      CPN = RHPI
      SPN = 0.0
      CSGN = CMPLX(-SPN,CPN)
      ZT = CMPLX(0.0E0,-FMM)
      RTOL = 1.0E0/TOL
      ASCLE = UFL*RTOL
      DO 120 I=1,NN
        ZN=CY(I)
        AA=REAL(ZN)
        BB=AIMAG(ZN)
        ATOL=1.0E0
        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 125
          ZN = ZN*CMPLX(RTOL,0.0E0)
          ATOL = TOL
  125   CONTINUE
        ZN = ZN*CSGN
        CY(I) = ZN*CMPLX(ATOL,0.0E0)
        CSGN = CSGN*ZT
  120 CONTINUE
      RETURN
      END

      SUBROUTINE CMLRI(Z, N, Y, NZ, TOL)
      COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z
      REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, RHO,
     * RHO2, SCLE, TFNF, TOL, TST, X, R1MACH
      INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, M, N
      DIMENSION Y(2)
      DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/
      SCLE = 1.0E+3*R1MACH(1)/TOL
      NZ=0
      AZ = CABS(Z)
      X = REAL(Z)
      IAZ = INT(AZ)
      INU = 1
      AT = FLOAT(IAZ) + 1.0E0
      CK = CMPLX(AT,0.0E0)/Z
      RZ = CTWO/Z
      P1 = CZERO
      P2 = CONE
      ACK = (AT+1.0E0)/AZ
      RHO = ACK + SQRT(ACK*ACK-1.0E0)
      RHO2 = RHO*RHO
      TST = (RHO2+RHO2)/((RHO2-1.0E0)*(RHO-1.0E0))
      TST = TST/TOL
      AK = AT
      DO 10 I=1,80
        PT = P2
        P2 = P1 - CK*P2
        P1 = PT
        CK = CK + RZ
        AP = CABS(P2)
        IF (AP.GT.TST*AK*AK) GO TO 20
        AK = AK + 1.0E0
   10 CONTINUE
      GO TO 110
   20 CONTINUE
      I = I + 1
      K = 0
      IF (INU.LT.IAZ) GO TO 40
      P1 = CZERO
      P2 = CONE
      AT = 2.0
      CK = CMPLX(AT,0.0E0)/Z
      ACK = AT/AZ
      TST = SQRT(ACK/TOL)
      ITIME = 1
      DO 30 K=1,80
        PT = P2
        P2 = P1 - CK*P2
        P1 = PT
        CK = CK + RZ
        AP = CABS(P2)
        IF (AP.LT.TST) GO TO 30
        IF (ITIME.EQ.2) GO TO 40
        ACK = CABS(CK)
        FLAM = ACK + SQRT(ACK*ACK-1.0E0)
        FKAP = AP/CABS(P1)
        RHO = AMIN1(FLAM,FKAP)
        TST = TST*SQRT(RHO/(RHO*RHO-1.0E0))
        ITIME = 2
   30 CONTINUE
      GO TO 110
   40 CONTINUE
      K = K + 1
      KK = MAX0(I+IAZ,K+INU)
      FKK = FLOAT(KK)
      P1 = CZERO
      P2 = CMPLX(SCLE,0.0E0)
      FNF = 0.0
      TFNF =0.0! FNF + FNF
      BK = 1.0
      SUM = CZERO
      KM = KK - INU
      DO 50 I=1,KM
        PT = P2
        P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
        P1 = PT
        AK = 1.0E0 - TFNF/(FKK+TFNF)
        ACK = BK*AK
        SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
        BK = ACK
        FKK = FKK - 1.0E0
   50 CONTINUE
      Y(2) = P2
        PT = P2
        P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2
        P1 = PT
        AK = 1.0E0 - TFNF/(FKK+TFNF)
        ACK = BK*AK
        SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1
        BK = ACK
        FKK = FKK - 1.0E0
        M = 1 
        Y(M) = P2
      PT = Z
      P1 = -CMPLX(FNF,0.0E0)*CLOG(RZ) + PT    ! FNF=0.0
      AP = 0.0 
      PT = P1
      P2 = P2 + SUM
      AP = CABS(P2)
      P1 = CMPLX(1.0E0/AP,0.0E0)
      CK = CEXP(PT)*P1
      PT = CONJG(P2)*P1
      CNORM = CK*PT
      DO 100 I=1,2
        Y(I) = Y(I)*CNORM
  100 CONTINUE
      RETURN
  110 CONTINUE
      NZ=-2
      RETURN
      END

      SUBROUTINE CBKNU(Z, N, Y, NZ, TOL, ELIM, ALIM)
      COMPLEX CCH, CK, COEF, CONE, CRSC, CS, CSCL, CSH, CSR, CSS, CTWO,
     * CZ, CZERO, F, P, PT, P1, P2, Q, RZ, SMU, ST, S1, S2, Y, Z,
     * ZD, CELM, CY
      REAL AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, CC, DNU,
     * DNU2, ELIM, ETEST, FC, FHS, FK, FKS, FPI, G1, G2, HPI, PI,
     * P2I, P2M, P2R, RK, RTHPI, R1, S, SPI, TM, TOL, TTH, T1, T2, XX,
     * YY, R1MACH, HELIM, ELM, XD, YD, ALAS, AS
      INTEGER I, IDUM, IFLAG, INU, K, KFLAG, KK, KMAX, KODED, N,
     * NZ, I1MACH, NW, J, IC, INUB
      DIMENSION BRY(3), CC(8), CSS(3), CSR(3), Y(2), CY(2)
      DATA KMAX / 30 /
      DATA R1 / 2.0E0 /
      DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/
      DATA PI, RTHPI, SPI ,HPI, FPI, TTH /
     1     3.14159265358979324E0,       1.25331413731550025E0,
     2     1.90985931710274403E0,       1.57079632679489662E0,
     3     1.89769999331517738E0,       6.66666666666666666E-01/
      DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
     1     5.77215664901532861E-01,    -4.20026350340952355E-02,
     2    -4.21977345555443367E-02,     7.21894324666309954E-03,
     3    -2.15241674114950973E-04,    -2.01348547807882387E-05,
     4     1.13302723198169588E-06,     6.11609510448141582E-09/
      XX = REAL(Z)
      YY = AIMAG(Z)
      CAZ = CABS(Z)
      CSCL = CMPLX(1.0E0/TOL,0.0E0)
      CRSC = CMPLX(TOL,0.0E0)
      CSS(1) = CSCL
      CSS(2) = CONE
      CSS(3) = CRSC
      CSR(1) = CRSC
      CSR(2) = CONE
      CSR(3) = CSCL
      BRY(1) = 1.0E+3*R1MACH(1)/TOL
      BRY(2) = 1.0E0/BRY(1)
      BRY(3) = R1MACH(2)
      NZ = 0
      IFLAG = 0
      KODED = 1
      RZ = CTWO/Z
      INU = 0
      DNU = 0.0 
      DNU2 = 0.0E0
      IF (CAZ.GT.R1) GO TO 110
      FC = 1.0E0
      SMU = CLOG(RZ)
      CSH=(0.0,0.0)
      CCH=(1.0,0.0)
      A2 = 1.0E0 
      T2 = 1.0
      T1 = 1.0E0/FC
      AK = 1.0E0
      S = CC(1)
      G1 = -S   
      G2 = 0.5E0*(T1+T2)*FC
      G1 = G1*FC
      F = CMPLX(G1,0.0E0)*CCH + SMU*CMPLX(G2,0.0E0)
      PT = (1.0,0.0)
      P = CMPLX(0.5E0,0.0E0)
      Q = CMPLX(0.5E0/T1,0.0E0)
      S1 = F
      S2 = P
      AK = 1.0E0
      A1 = 1.0E0
      CK = CONE
      BK = 1.0E0
      IF (CAZ.LT.TOL) GO TO 100
      CZ = Z*Z*CMPLX(0.25E0,0.0E0)
      T1 = 0.25E0*CAZ*CAZ
   90 CONTINUE
      F = (F*CMPLX(AK,0.0E0)+P+Q)*CMPLX(1.0E0/BK,0.0E0)
      P = P*CMPLX(1.0E0/(AK-DNU),0.0E0)
      Q = Q*CMPLX(1.0E0/(AK+DNU),0.0E0)
      RK = 1.0E0/AK
      CK = CK*CZ*CMPLX(RK,0.0E0)
      S1 = S1 + CK*F
      S2 = S2 + CK*(P-F*CMPLX(AK,0.0E0))
      A1 = A1*T1*RK
      BK = BK + AK + AK + 1.0E0
      AK = AK + 1.0E0
      IF (A1.GT.TOL) GO TO 90
  100 CONTINUE
      KFLAG = 2
      BK = REAL(SMU)
      A1 = 1.0E0
      AK = A1*ABS(BK)
      IF (AK.GT.ALIM) KFLAG = 3
      P2 = S2*CSS(KFLAG)
      S2 = P2*RZ
      S1 = S1*CSS(KFLAG)
      IF (KODED.EQ.1) GO TO 210
      F = CEXP(Z)
      S1 = S1*F
      S2 = S2*F
      GO TO 210
  110 CONTINUE
      COEF = CMPLX(RTHPI,0.0E0)/CSQRT(Z)
      KFLAG = 2
      IF (KODED.EQ.2) GO TO 120
      IF (XX.GT.ALIM) GO TO 290
      A1 = EXP(-XX)*REAL(CSS(KFLAG))
      PT = CMPLX(A1,0.0E0)*CMPLX(COS(YY),-SIN(YY))
      COEF = COEF*PT
  120 CONTINUE
      AK = 1.0  
      FHS = ABS(0.25E0-DNU2)
      T1 = FLOAT(I1MACH(11)-1)*R1MACH(5)*3.321928094E0
      T1 = AMAX1(T1,12.0E0)
      T1 = AMIN1(T1,60.0E0)
      T2 = TTH*T1 - 6.0E0
      IF (XX.NE.0.0E0) GO TO 130
      T1 = HPI
      GO TO 140
  130 CONTINUE
      T1 = ABS(ATAN(YY/XX))
  140 CONTINUE
      IF (T2.GT.CAZ) GO TO 170
      ETEST = AK/(PI*CAZ*TOL)
      FK = 1.0E0
      IF (ETEST.LT.1.0E0) GO TO 180
      FKS = 2.0E0
      RK = CAZ + CAZ + 2.0E0
      A1 = 0.0E0
      A2 = 1.0E0
      DO 150 I=1,KMAX
        AK = FHS/FKS
        BK = RK/(FK+1.0E0)
        TM = A2
        A2 = BK*A2 - AK*A1
        A1 = TM
        RK = RK + 2.0E0
        FKS = FKS + FK + FK + 2.0E0
        FHS = FHS + FK + FK
        FK = FK + 1.0E0
        TM = ABS(A2)*FK
        IF (ETEST.LT.TM) GO TO 160
  150 CONTINUE
      GO TO 310
  160 CONTINUE
      FK = FK + SPI*T1*SQRT(T2/CAZ)
      FHS = ABS(0.25E0-DNU2)
      GO TO 180
  170 CONTINUE
      A2 = SQRT(CAZ)
      AK = FPI*AK/(TOL*SQRT(A2))
      AA = 3.0E0*T1/(1.0E0+CAZ)
      BB = 14.7E0*T1/(28.0E0+CAZ)
      AK = (ALOG(AK)+CAZ*COS(AA)/(1.0E0+0.008E0*CAZ))/COS(BB)
      FK = 0.12125E0*AK*AK/CAZ + 1.5E0
  180 CONTINUE
      K = INT(FK)
      FK = FLOAT(K)
      FKS = FK*FK
      P1 = CZERO
      P2 = CMPLX(TOL,0.0E0)
      CS = P2
      DO 190 I=1,K
        A1 = FKS - FK
        A2 = (FKS+FK)/(A1+FHS)
        RK = 2.0E0/(FK+1.0E0)
        T1 = (FK+XX)*RK
        T2 = YY*RK
        PT = P2
        P2 = (P2*CMPLX(T1,T2)-P1)*CMPLX(A2,0.0E0)
        P1 = PT
        CS = CS + P2
        FKS = A1 - FK + 1.0E0
        FK = FK - 1.0E0
  190 CONTINUE
      TM = CABS(CS)
      PT = CMPLX(1.0E0/TM,0.0E0)
      S1 = PT*P2
      CS = CONJG(CS)*PT
      S1 = COEF*S1*CS
      TM = CABS(P2)
      PT = CMPLX(1.0E0/TM,0.0E0)
      P1 = PT*P1
      P2 = CONJG(P2)*PT
      PT = P1*P2
      S2 = S1*(CONE+(CMPLX(DNU+0.5E0,0.0E0)-PT)/Z)
  210 CONTINUE
      CK =  RZ
      ZD = Z
      IF(IFLAG.EQ.1) GO TO 270
      Y(1) = S1*CSR(KFLAG)
      Y(2) = S2*CSR(KFLAG)
      RETURN
  270 CONTINUE
      Y(1) = S1
      Y(2) = S2
      ASCLE = BRY(1)
      CALL CKSCL(ZD, 2, Y, NZ, RZ, ASCLE, TOL, ELIM)
      INU = 2 - NZ
      IF (INU.LE.0) RETURN
      KK = NZ + 1
      S1 = Y(KK)
      Y(KK) = S1*CSR(1)
      IF (INU.EQ.1) RETURN
      KK = NZ + 2
      S2 = Y(KK)
      Y(KK) = S2*CSR(1)
      RETURN
  290 CONTINUE
      KODED = 2
      IFLAG = 1
      KFLAG = 2
      GO TO 120
  310 CONTINUE
      NZ = -2
      RETURN
      END

      SUBROUTINE CKSCL(ZR, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
      COMPLEX CK, CS, CY, CZERO, RZ, S1, S2, Y, ZR, ZD, CELM
      REAL AA, ASCLE, ACS, AS, CSI, CSR, ELIM, FN, TOL, XX, ZRI,
     * ELM, ALAS, HELIM
      INTEGER I, IC, K, KK, N, NW, NZ
      DIMENSION Y(2), CY(2)
      DATA CZERO / (0.0E0,0.0E0) /
      NZ = 0
      IC = 0
      XX = REAL(ZR)
      DO 10 I=1,2
        S1 = Y(I)
        CY(I) = S1
        AS = CABS(S1)
        ACS = -XX + ALOG(AS)
        NZ = NZ + 1
        Y(I) = CZERO
        IF (ACS.LT.(-ELIM)) GO TO 10
        CS = -ZR + CLOG(S1)
        CSR = REAL(CS)
        CSI = AIMAG(CS)
        AA = EXP(CSR)/TOL
        CS = CMPLX(AA,0.0E0)*CMPLX(COS(CSI),SIN(CSI))
        CALL CUCHK(CS, NW, ASCLE, TOL)
        IF (NW.NE.0) GO TO 10
        Y(I) = CS
        NZ = NZ - 1
        IC = I
   10 CONTINUE
      IF (IC.GT.1) GO TO 20
      Y(1) = CZERO
      NZ = 2
   20 CONTINUE
      RETURN
      END

      SUBROUTINE CUCHK(Y, NZ, ASCLE, TOL)
      COMPLEX Y
      REAL ASCLE, SS, ST, TOL, YR, YI
      INTEGER NZ
      NZ = 0
      YR = REAL(Y)
      YI = AIMAG(Y)
      YR = ABS(YR)
      YI = ABS(YI)
      ST = AMIN1(YR,YI)
      IF (ST.GT.ASCLE) RETURN
      SS = AMAX1(YR,YI)
      ST=ST/TOL
      IF (SS.LT.ST) NZ = 1
      RETURN
      END

      INTEGER FUNCTION I1MACH(I)
      INTEGER IMACH(16),OUTPUT
      SAVE IMACH
      EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR THE IBM PC
C
      DATA IMACH( 1) /     5 /
      DATA IMACH( 2) /     6 /
      DATA IMACH( 3) /     0 /
      DATA IMACH( 4) /     0 /
      DATA IMACH( 5) /    32 /
      DATA IMACH( 6) /     4 /
      DATA IMACH( 7) /     2 /
      DATA IMACH( 8) /    31 /
      DATA IMACH( 9) / 2147483647 /
      DATA IMACH(10) /     2 /
      DATA IMACH(11) /    24 /
      DATA IMACH(12) /  -125 /
      DATA IMACH(13) /   127 /
      DATA IMACH(14) /    53 /
      DATA IMACH(15) / -1021 /
      DATA IMACH(16) /  1023 /

      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
C
      I1MACH = IMACH(I)
      RETURN
C
   10 CONTINUE
      WRITE (UNIT = OUTPUT, FMT = 9000)
 9000 FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
      STOP
      END

      REAL FUNCTION R1MACH(I)
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
      REAL RMACH(5)
      SAVE RMACH
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE IBM PC
C
      DATA SMALL(1) /     8420761 /
      DATA LARGE(1) /  2139081118 /
      DATA RIGHT(1) /   863997169 /
      DATA DIVER(1) /   872385777 /
      DATA LOG10(1) /  1050288283 /
C
      IF (I .LT. 1  .OR.  I .GT. 5)
     1   WRITE(*,*) 'R1MACH -- I OUT OF BOUNDS [1,5]', i
      R1MACH = RMACH(I)
      RETURN
      END