PROGRAM CQCBI(INPUT,OUTPUT,TAPE7=OUTPUT) C C CQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION C GENERATED BY SUBROUTINE CBESI. C C CQCBI 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 COMPLEX CK, CONE, CSGN, CW, CY, W, Y, Z, ZN, ZSC, ZT, ZZ REAL AA, AER, ALIM, ARG, ATOL, AW, C, CARG, C1, DIG, ELIM, EPS, * ER, ERTOL, FAC, FNU, FNUL, GNU, HPI, PI, R, RL, RLT, R1, R1M5, * R2, RM, S, SARG, SEPS, S1, T, TOL, TPI, XX, YY, R1MACH 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), Y(20), W(20), CK(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 = 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.2*DIG + 3.0 RM = 0.5E0*(ALIM+ELIM) RM=AMIN1(RM,650.0E0) R2 = AMIN1(RM,FNUL) R1 = 2.0*SQRT(FNUL+1.0) C----------------------------------------------------------------------- WRITE (LUN,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM CB, * 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 (6E12.4/) ERTOL = TOL*1.0E+5 ZZ = CMPLX(1.0E0/TOL,0.0E0) CONE = CMPLX(1.0,0.0) FAC = 100. ATOL = FAC*TOL HPI = 2.0*ATAN(1.0E0) PI = HPI + HPI TPI = PI + PI I = 1 IA = 4 EPS = 0.01 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*FLOAT(K-1)/FLOAT(IL-1) C1 = COS(T) S1 = SIN(T) IF (ABS(C1).LT.ATOL) C1 = 0.0 IF (ABS(S1).LT.ATOL) S1 = 0.0 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.0*FLOAT(3-IR)+RL*FLOAT(IR-1))/2.0 GNU = R*R/4.0 - 0.2 - FLOAT(N-1) FNU = AMAX1(0.0,GNU) GO TO 110 60 CONTINUE R = (RL*FLOAT(3-IR)+R2*FLOAT(IR-1))/2.0 GNU = SQRT(R+R) - 0.2 - FLOAT(N-1) FNU = AMAX1(0.0,GNU) GO TO 110 70 CONTINUE IF (R2.EQ.RM) GO TO 220 R = (R2*FLOAT(3-IR)+RM*FLOAT(IR-1))/2.0 GNU = SQRT(R+R) - 0.2 - FLOAT(N-1) FNU = AMAX1(0.0,GNU) GO TO 110 80 CONTINUE IF (R1.GE.RL) GO TO 220 R = (R1*FLOAT(3-IR)+RL*FLOAT(IR-1))/2.0 FNU = FNUL - 0.2 - FLOAT(N-1) GO TO 110 90 CONTINUE R = (RL*FLOAT(3-IR)+R2*FLOAT(IR-1))/2.0 FNU = FNUL - 0.2 - FLOAT(N-1) GO TO 110 100 CONTINUE IF (R2.EQ.RM) GO TO 220 R = (R2*FLOAT(3-IR)+RM*FLOAT(IR-1))/2.0 FNU = FNUL - 0.2 - FLOAT(N-1) 110 CONTINUE DO 180 IT=1,ITL Z = CMPLX(R*C(IT),R*S(IT)) CALL CBESI(Z, FNU, KODE, N1, Y, 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----------------------------------------------------------------------- ZN = Z XX = REAL(Z) YY = AIMAG(Z) IF (XX.GE.0.0) GO TO 120 ZN = -Z INU = INT(FNU) ARG = (FNU-FLOAT(INU))*PI IF (YY.LT.0.0) ARG = -ARG CARG = COS(ARG) SARG = SIN(ARG) CSGN = CMPLX(CARG,SARG) IF (MOD(INU,2).EQ.1) CSGN = -CSGN 120 CONTINUE CALL CWRSK(ZN, FNU, KODE, N, W, NZ2, CK, TOL, ELIM, * ALIM) IF (XX.GE.0.0) GO TO 150 DO 130 I=1,N W(I) = W(I)*CSGN CSGN = -CSGN 130 CONTINUE GO TO 150 140 CONTINUE CALL CBESI(Z, FNU, KODE, N, W, NZ2, IERR) 150 CONTINUE MFLG = 0 DO 160 I=1,N ZT = W(I) ZSC = ZZ IF (CABS(ZT).GT.1.0E0) ZSC = CMPLX(TOL,0.0E0) CW = W(I)*ZSC CY = Y(I)*ZSC ER = CABS(CY-CW) AW = CABS(CW) IF (AW.NE.0.0) ER = ER/AW IF (AW.EQ.0.0) ER = ER/CABS(ZSC) 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 =, E12.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) Z, FNU, Y(KK), W(KK) 99991 FORMAT (7E12.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/) Z = CMPLX(1.4,1.4) IPRNT = 0 DO 280 I=1,2 FNU = 10.2 KODE = 1 N = 20 230 CONTINUE CALL CBESI(Z, FNU, KODE, N, Y, NZI, IERR) IF (NZI.NE.0) GO TO 240 FNU = FNU + 5.0 GO TO 230 240 CONTINUE IF (NZI.LT.10) GO TO 250 FNU = FNU - 1.0 GO TO 230 250 CONTINUE CALL CBESK(Z, FNU, KODE, 2, W, NZK, IERR) ZT = CONE/Z CY = W(1)*Y(2) CW = W(2)*Y(1) CW = CW + CY - ZT ER = CABS(CW)/CABS(ZT) 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, Z, FNU, KODE, N 99987 FORMAT (4E12.4, 2I5) 270 CONTINUE RLT = RL + RL Z = CMPLX(RLT,0.0) 280 CONTINUE C----------------------------------------------------------------------- C CHECK NEAR OVERFLOW LIMITS C----------------------------------------------------------------------- Z = CMPLX(ELIM,0.0) FNU = 0.0 290 CONTINUE CALL CBESK(Z, FNU, KODE, N, Y, NZK, IERR) IF (NZK.LT.10) GO TO 300 IF (NZK.EQ.N) FNU = FNU + 3.0E0 FNU = FNU + 2.0D0 GO TO 290 300 CONTINUE GNU = FNU + FLOAT(N-2) CALL CBESI(Z, GNU, KODE, 2, W, NZI, IERR) ZT = CONE/Z CY = Y(N-1)*W(2) CW = Y(N)*W(1) CW = CW + CY - ZT ER = CABS(CW)/CABS(ZT) 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, Z, FNU, KODE, N 320 CONTINUE IF (IPRNT.EQ.0) WRITE (LUN,99990) STOP END .