PROGRAM ZQCBI(INPUT,OUTPUT,TAPE7=OUTPUT)
C
C                 *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION
C     GENERATED BY SUBROUTINE CBESI.
C
C     ZQCBI GENERATES SEQUENCES CROSSING FORMULA BOUNDARIES WHICH
C     ARE STARTED BY ONE FORMULA AND TERMINATED IN A REGION WHERE
C     ANOTHER FORMULA APPLIES. THE TERMINATED VALUE IS CHECKED BY
C     THE FORMULA APPROPRIATE TO THAT REGION.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX CK,CONE,CSGN,CW,CY,W,Y,Z,ZN,ZSC,ZT,ZZ
      DOUBLE PRECISION AA, AER, ALIM, ARG, ATOL, AW, C, CARG, CKI, CKR,
     * CONEI, CONER, CSGNI, CSGNR, CWI, CWR, CYI, CYR, C1, DIG, ELIM,
     * EPS, ER, ERTOL, FAC, FNU, FNUL, FNULT, GNU, HPI, PI, R, RL, RLT,
     * R1, R1M5, R2, RM, S, SARG, SEPS, STI, STR, S1, T, TOL, TPI, WI,
     * WR, YI, YR, ZI, ZNI, ZNR, ZR, ZSCR, ZTI, ZTR, ZZR, D1MACH, ZABS
      INTEGER I, IA, ICASE, IL, INU, IPRNT, IR, IT, ITL, K, KK, KODE,
     * K1, K2, LFLG, LUN, MFLG, N, NZI, NZK, NZ1, NZ2, N1, IERR, I1MACH
      DIMENSION C(25), S(25), AER(25), YR(20), YI(20), WR(20), WI(20),
     * CKR(2), CKI(2)
      DATA LUN /7/
C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      TOL = DMAX1(D1MACH(4),1.0D-18)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN0(IABS(K1),IABS(K2))
      ELIM = 2.303E0*(DBLE(FLOAT(K))*R1M5-3.0E0)
      K1 = I1MACH(14) - 1
      AA = R1M5*DBLE(FLOAT(K1))
      DIG = DMIN1(AA,18.0D0)
      AA = AA*2.303D0
      ALIM = ELIM + DMAX1(-AA,-41.45D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      RM = 0.5D0*(ALIM+ELIM)
      RM=DMIN1(RM,650.0D0)
      R2 = DMIN1(RM,FNUL)
      R1 = 2.0D0*DSQRT(FNUL+1.0D0)
C-----------------------------------------------------------------------
      WRITE (LUN,99999)
99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM ZB,
     * 3HESI/)
      WRITE (LUN,99998)
99998 FORMAT (37H PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG)
      WRITE (LUN,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      FNULT = FNUL + 20.0D0
      ERTOL = TOL*1.0D+5
      ZZR = 1.0D0/TOL
      CONER = 1.0D0
      CONEI = 0.0D0
      FAC = 100.0D0
      ATOL = FAC*TOL
      HPI = 2.0D0*DATAN(1.0D0)
      PI = HPI + HPI
      TPI = PI + PI
      I = 1
      IA = 4
      EPS = 0.01D0
      SEPS = EPS
      IL = 13
C-----------------------------------------------------------------------
C     TEST 17 VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
      WRITE (LUN,99996)
99996 FORMAT (/33H CHECKS ACROSS FORMULA BOUNDARIES/)
      DO 40 K=1,IL
        IF (K.NE.7) GO TO 10
        IA = 10
        SEPS = -SEPS
   10   CONTINUE
        IF (IABS(K-IA).EQ.2) GO TO 40
        T = -PI + TPI*DBLE(FLOAT(K-1))/DBLE(FLOAT(IL-1))
        C1 = DCOS(T)
        S1 = DSIN(T)
        IF (DABS(C1).LT.ATOL) C1 = 0.0D0
        IF (DABS(S1).LT.ATOL) S1 = 0.0D0
        C(I) = C1
        S(I) = S1
        IF (IABS(K-IA).GT.1) GO TO 30
        IF (IABS(K-IA).EQ.1) GO TO 20
        C(I) = C1 - SEPS
        S(I) = S1
        C(I+1) = C1
        S(I+1) = S1
        C(I+2) = C1 + SEPS
        S(I+2) = S1
        I = I + 3
        GO TO 40
   20   CONTINUE
        C(I) = C1 - SEPS
        S(I) = S1
        C(I+1) = C1 + SEPS
        S(I+1) = S1
        I = I + 2
        GO TO 40
   30   CONTINUE
        I = I + 1
   40 CONTINUE
      S(1) = -EPS
      ITL = I - 1
      LFLG = 0
      DO 220 ICASE=1,6
        DO 210 KODE=1,2
          DO 200 N=1,4
            N1 = N + 2
            DO 190 IR=1,3
              GO TO (50, 60, 70, 80, 90, 100), ICASE
   50         CONTINUE
              R = (2.0D0*DBLE(FLOAT(3-IR))+RL*DBLE(FLOAT(IR-1)))/2.0D0
              GNU = R*R/4.0D0 - 1.2D0 - DBLE(FLOAT(N-1))
              FNU = DMAX1(0.0D0,GNU)
              GO TO 110
   60         CONTINUE
              R = (RL*DBLE(FLOAT(3-IR))+R2*DBLE(FLOAT(IR-1)))/2.0D0
              GNU = DSQRT(R+R) - 0.2D0 - DBLE(FLOAT(N-1))
              FNU = DMAX1(0.0D0,GNU)
              GO TO 110
   70         CONTINUE
              IF (R2.EQ.RM) GO TO 220
              R = (R2*DBLE(FLOAT(3-IR))+RM*DBLE(FLOAT(IR-1)))/2.0D0
              GNU = DSQRT(R+R) - 0.2D0 - DBLE(FLOAT(N-1))
              FNU = DMAX1(0.0D0,GNU)
              GO TO 110
   80         CONTINUE
              IF (R1.GE.RL) GO TO 220
              R = (R1*DBLE(FLOAT(3-IR))+RL*DBLE(FLOAT(IR-1)))/2.0D0
              FNU = FNUL - 0.2D0 - DBLE(FLOAT(N-1))
              GO TO 110
   90         CONTINUE
              R = (RL*DBLE(FLOAT(3-IR))+R2*DBLE(FLOAT(IR-1)))/2.0D0
              FNU = FNUL - 0.2D0 - DBLE(FLOAT(N-1))
              GO TO 110
  100         CONTINUE
              IF (R2.EQ.RM) GO TO 220
              R = (R2*DBLE(FLOAT(3-IR))+RM*DBLE(FLOAT(IR-1)))/2.0D0
              FNU = FNUL - 0.2D0 - DBLE(FLOAT(N-1))
  110         CONTINUE
              DO 180 IT=1,ITL
                ZR = R*C(IT)
                ZI = R*S(IT)
                CALL ZBESI(ZR, ZI, FNU, KODE, N1, YR, YI, NZ1, IERR)
                IF (NZ1.NE.0) GO TO 180
                IF (ICASE.NE.6) GO TO 140
C-----------------------------------------------------------------------
C     COMPARE VALUES FOR R.GT.FNUL AND FNU.LT.FNUL WITH WRONSKIAN
C     EVALUATION. TO COMPARE ALL CASES WITH WRONSKIAN EVALUATION,
C     CHANGE ICASE.NE.6 TO ICASE.LT.0 ABOVE.
C-----------------------------------------------------------------------
                ZNR = ZR
                ZNI = ZI
                IF (ZR.GE.0.0D0) GO TO 120
                ZNR = -ZR
                ZNI = -ZI
                INU = INT(SNGL(FNU))
                ARG = (FNU-DBLE(FLOAT(INU)))*PI
                IF (ZI.LT.0.0D0) ARG = -ARG
                CARG = DCOS(ARG)
                SARG = DSIN(ARG)
                CSGNR = CARG
                CSGNI = SARG
                IF (MOD(INU,2).EQ.0) GO TO 120
                CSGNR = -CSGNR
                CSGNI = -CSGNI
  120           CONTINUE
                CALL ZWRSK(ZNR, ZNI, FNU, KODE, N, WR, WI, NZ2, CKR,
     *           CKI, TOL, ELIM, ALIM)
                IF (ZR.GE.0.0D0) GO TO 150
                DO 130 I=1,N
                  STR = WR(I)*CSGNR - WI(I)*CSGNI
                  WI(I) = WR(I)*CSGNI + WI(I)*CSGNR
                  WR(I) = STR
                  CSGNR = -CSGNR
                  CSGNI = -CSGNI
  130           CONTINUE
                GO TO 150
  140           CONTINUE
                CALL ZBESI(ZR, ZI, FNU, KODE, N, WR, WI, NZ2, IERR)
  150           CONTINUE
                MFLG = 0
                DO 160 I=1,N
                  ZTR = WR(I)
                  ZTI = WI(I)
                  ZSCR = ZZR
                  IF (DABS(ZTR).GT.1.0D0 .OR. DABS(ZTI).GT.1.0D0) ZSCR
     *             = TOL
                  CWR = WR(I)*ZSCR
                  CWI = WI(I)*ZSCR
                  CYR = YR(I)*ZSCR
                  CYI = YI(I)*ZSCR
                  STR = CYR - CWR
                  STI = CYI - CWI
                  ER = ZABS(STR,STI)
                  AW = ZABS(CWR,CWI)
                  IF (AW.NE.0.0D0) ER = ER/AW
                  IF (AW.EQ.0.0D0) ER = ER/DABS(ZSCR)
                  AER(I) = ER
                  IF (ER.GT.ERTOL .OR. NZ1.NE.0 .OR. NZ2.NE.0) MFLG = 1
  160           CONTINUE
                IF (MFLG.EQ.0) GO TO 180
                IF (LFLG.EQ.1) GO TO 170
                WRITE (LUN,99995) ERTOL
99995           FORMAT (/43H CASES WHICH UNDERFLOW OR VIOLATE THE RELAT,
     *           19HIVE ERROR TEST WITH/8H ERTOL =, D12.4/)
                WRITE (LUN,99994)
99994           FORMAT (/14H OUTPUT FORMAT/25H KODE,N,IR,IT,NZ1,NZ2,ICA,
     *           2HSE)
                WRITE (LUN,99993)
99993           FORMAT (12H ER(K),K=1,N/18H Z,FNU,Y(KK),W(KK)/)
                LFLG = 1
  170           CONTINUE
                KK = MAX0(NZ1,NZ2) + 1
                KK = MIN0(N,KK)
                WRITE (LUN,99992) KODE, N, IR, IT, NZ1, NZ2, ICASE
99992           FORMAT (8I5)
                WRITE (LUN,99991) (AER(K),K=1,N)
                WRITE (LUN,99991) ZR, ZI, FNU, YR(KK), YI(KK), WR(KK),
     *           WI(KK)
99991           FORMAT (7D12.4)
  180         CONTINUE
  190       CONTINUE
  200     CONTINUE
  210   CONTINUE
  220 CONTINUE
      IF (LFLG.EQ.0) WRITE (LUN,99990)
99990 FORMAT (/16H QUICK CHECKS OK/)
C-----------------------------------------------------------------------
C     CHECKS NEAR UNDERFLOW LIMITS ON SERIES(I=1) AND UNIFORM
C     ASYMPTOTIC EXPANSION(I=2)
C-----------------------------------------------------------------------
      WRITE (LUN,99989)
99989 FORMAT (/42H CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS/)
      ZR = 1.4D0
      ZI = 1.4D0
      IPRNT = 0
      DO 280 I=1,2
        FNU = 10.2D0
        KODE = 1
        N = 20
  230   CONTINUE
        CALL ZBESI(ZR, ZI, FNU, KODE, N, YR, YI, NZI, IERR)
        IF (NZI.NE.0) GO TO 240
        FNU = FNU + 5.0D0
        GO TO 230
  240   CONTINUE
        IF (NZI.LT.10) GO TO 250
        FNU = FNU - 1.0D0
        GO TO 230
  250   CONTINUE
        CALL ZBESK(ZR, ZI, FNU, KODE, 2, WR, WI, NZK, IERR)
        CALL ZDIV(CONER, CONEI, ZR, ZI, ZTR, ZTI)
        CYR = WR(1)*YR(2) - WI(1)*YI(2)
        CYI = WR(1)*YI(2) + WI(1)*YR(2)
        CWR = WR(2)*YR(1) - WI(2)*YI(1)
        CWI = WR(2)*YI(1) + WI(2)*YR(1)
        CWR = CWR + CYR - ZTR
        CWI = CWI + CYI - ZTI
        ER = ZABS(CWR,CWI)/ZABS(ZTR,ZTI)
        IF (ER.LT.ERTOL) GO TO 270
        IF (IPRNT.EQ.1) GO TO 260
        WRITE (LUN,99988)
99988   FORMAT (/14H OUTPUT FORMAT/19H ERROR,Z,FNU,KODE,N/)
        IPRNT = 1
  260   CONTINUE
        WRITE (LUN,99987) ER, ZR, ZI, FNU, KODE, N
99987   FORMAT (4D12.4, 2I5)
  270   CONTINUE
        RLT = RL + RL
        ZR = RLT
        ZI = 0.0D0
  280 CONTINUE
C-----------------------------------------------------------------------
C     CHECK NEAR OVERFLOW LIMITS
C-----------------------------------------------------------------------
      ZR = ELIM
      ZI = 0.0D0
      FNU = 0.0D0
  290 CONTINUE
      CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZK, IERR)
      IF (NZK.LT.10) GO TO 300
      IF (NZK.EQ.N) FNU = FNU + 3.0D0
      FNU = FNU + 2.0D0
      GO TO 290
  300 CONTINUE
      GNU = FNU + DBLE(FLOAT(N-2))
      CALL ZBESI(ZR, ZI, GNU, KODE, 2, WR, WI, NZI, IERR)
      CALL ZDIV(CONER, CONEI, ZR, ZI, ZTR, ZTI)
      CYR = YR(N-1)*WR(2) - YI(N-1)*WI(2)
      CYI = YR(N-1)*WI(2) + YI(N-1)*WR(2)
      CWR = YR(N)*WR(1) - YI(N)*WI(1)
      CWI = YR(N)*WI(1) + YI(N)*WR(1)
      CWR = CWR + CYR - ZTR
      CWI = CWI + CYI - ZTI
      ER = ZABS(CWR,CWI)/ZABS(ZTR,ZTI)
      IF (ER.LT.ERTOL) GO TO 320
      IF (IPRNT.EQ.1) GO TO 310
      WRITE (LUN,99988)
      IPRNT = 1
  310 CONTINUE
      WRITE (LUN,99987) ER, ZR, ZI, FNU, KODE, N
  320 CONTINUE
      IF (IPRNT.EQ.0) WRITE (LUN,99990)
      STOP
      END

.