PROGRAM CQCBJ(INPUT,OUTPUT,TAPE7=OUTPUT) C C CQCBJ IS A QUICK CHECK ROUTINE FOR THE COMPLEX J BESSEL FUNCTION C GENERATED BY SUBROUTINE CBESJ. C C CQCBJ GENERATES SEQUENCES OF J BESSEL FUNCTIONS FROM C CBESJ AND CHECKS THEM AGAINST THE EVALUATION FROM THE FORMULA C C J(FNU,Z)=0.5*( H(1,FNU,Z) + H(2,FNU,Z) ) C -PI.LT.ARG(Z).LE.PI C C FOR CABS(Z).GE.FNU. FOR CABS(Z).LT.FNU, THE FIRST N MEMBERS OF C A SEQUENCE OF LENGTH N+16 ARE CHECKED AGAINST A CORRESPONDING N C MEMBER SEQUENCE WHERE BOTH SEQUENCES ARE GENERATED BY CBESJ C BEGINNING AT ORDER FNU. 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 CHALF, COE1, COE2, CW, V, W, Y, Z REAL AA, AER, ALIM, ATOL, AV, C, CC, C1, DD, DIG, ELIM, EPS, ER, * ERTOL, FAC, FNU, FNUL, GNU, HPI, PI, R, RL, RM, R1M5, R2, S, * SEPS, S1, T, TOL, TPI, XNU, XX, YY, R1MACH INTEGER I, IA, ICASE, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, K2, * LFLG, LUN, M, MFLG, N, NU, NZ, NZ1, NZ2, IERR, I1MACH DIMENSION C(25), S(25), AER(25), XNU(7), Y(20), W(20), V(20) 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.2E0*DIG + 3.0E0 RM = 0.5E0*(ALIM+ELIM) RM=AMIN1(RM,650.0E0) R2 = AMIN1(RM,FNUL) C----------------------------------------------------------------------- WRITE (LUN,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM CB, * 3HESJ/) 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 CHALF = CMPLX(0.5E0,0.0E0) FAC = 100.0E0 ATOL = FAC*TOL HPI = 2.0E0*ATAN(1.0E0) PI = HPI + HPI TPI = PI + PI I = 1 IA = 4 EPS = 0.01E0 SEPS = EPS IL = 13 C----------------------------------------------------------------------- C TEST 17 VALUES OF Z IN -PI.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- WRITE (LUN,99996) 99996 FORMAT (/28H CHECKS IN THE (Z,FNU) SPACE/) 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.0E0 IF (ABS(S1).LT.ATOL) S1 = 0.0E0 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 XNU(1) = 0.0E0 XNU(2) = 0.3E0 XNU(3) = 0.7E0 XNU(4) = 1.0E0 XNU(5) = 1.3E0 XNU(6) = 1.7E0 XNU(7) = FNUL + 1.1E0 LFLG = 0 DO 260 KODE=1,2 DO 250 N=1,4 DO 240 NU=1,7 FNU = XNU(NU) DO 230 ICASE=1,3 IRB = MIN0(2,ICASE) DO 220 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (EPS*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0 GO TO 80 60 CONTINUE R = (2.0E0*FLOAT(4-IR)+R2*FLOAT(IR-1))/3.0E0 GO TO 80 70 CONTINUE IF (R2.EQ.RM) GO TO 230 R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 80 CONTINUE GNU = FNU + FLOAT(N-1) C----------------------------------------------------------------------- C 10.*AA IS APPROX. LOCAL MAX FOR REAL OR NEAR REAL ARUMENTS WHEN C GNU.LT.R C----------------------------------------------------------------------- AA = 0.079788E0/SQRT(R) DO 210 IT=1,ITL Z = CMPLX(R*C(IT),R*S(IT)) IF (R.LT.GNU) GO TO 140 C----------------------------------------------------------------------- C CASES FOR CABS(Z).GE.FNU+N-1 C----------------------------------------------------------------------- CALL CBESJ(Z, FNU, KODE, N, V, NZ, IERR) IF (NZ.NE.0) GO TO 210 CALL CBESH(Z, FNU, KODE, 1, N, W, NZ1, IERR) CALL CBESH(Z, FNU, KODE, 2, N, Y, NZ2, IERR) IF (KODE.EQ.1) GO TO 160 C----------------------------------------------------------------------- C SCALING ADJUSTMENTS IN H FUNCTIONS FOR KODE=2 C----------------------------------------------------------------------- XX = REAL(Z) YY = AIMAG(Z) CC = -YY - ABS(YY) IF (CC.GT.(-ALIM)) GO TO 90 COE1 = CMPLX(0.0E0,0.0E0) GO TO 100 90 CONTINUE CW = CMPLX(CC,XX) COE1 = CEXP(CW) 100 CONTINUE DD = YY - ABS(YY) IF (DD.GT.(-ALIM)) GO TO 110 COE2 = CMPLX(0.0E0,0.0E0) GO TO 120 110 CONTINUE CW = CMPLX(DD,-XX) COE2 = CEXP(CW) 120 CONTINUE DO 130 KK=1,N Y(KK) = Y(KK)*COE2 W(KK) = W(KK)*COE1 130 CONTINUE GO TO 160 C----------------------------------------------------------------------- C CASES FOR CABS(Z).LT.FNU+N-1 C----------------------------------------------------------------------- 140 CONTINUE M = N + 16 CALL CBESJ(Z, FNU, KODE, M, V, NZ, IERR) IF (NZ.GT.10) GO TO 210 CALL CBESJ(Z, FNU, KODE, N, W, NZ, IERR) DO 150 KK=1,N Y(KK) = W(KK) 150 CONTINUE 160 CONTINUE MFLG = 0 DO 190 I=1,N CW = (W(I)+Y(I))*CHALF AV = CABS(V(I)) ER = CABS(CW-V(I)) IF ((IT.EQ.1) .OR. (IT.EQ.9) .OR. (IT.EQ.17)) GO TO * 170 ER = ER/AV GO TO 180 170 CONTINUE IF ((GNU.LT.R) .AND. (AV.LT.AA)) GO TO 180 ER = ER/AV 180 CONTINUE AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 190 CONTINUE IF (MFLG.EQ.0) GO TO 210 IF (LFLG.EQ.1) GO TO 200 WRITE (LUN,99995) ERTOL 99995 FORMAT (/41H CASES WHICH VIOLATE THE RELATIVE ERROR T, * 16HEST WITH ERTOL =, E12.4/) WRITE (LUN,99994) 99994 FORMAT (/14H OUTPUT FORMAT/23H KODE,N,IR,IT,NZ1,NZ2,I, * 4HCASE) WRITE (LUN,99993) 99993 FORMAT (12H ER(K),K=1,N/26H Z,FNU,Y(KK),W(KK),V(KK), , * 35HKK=INDEX OF FIRST NON-ZERO Y,W PAIR/) LFLG = 1 200 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), V(KK) 99991 FORMAT (9E12.4) 210 CONTINUE 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE IF (LFLG.EQ.0) WRITE (LUN,99990) 99990 FORMAT (/16H QUICK CHECKS OK/) STOP END .