PROGRAM CQCBH(INPUT,OUTPUT,TAPE7=OUTPUT) C C CQCBH IS A QUICK CHECK ROUTINE FOR THE COMPLEX H BESSEL FUNCTIONS C GENERATED BY SUBROUTINE CBESH. C C CQCBH GENERATES SEQUENCES OF H BESSEL FUNCTIONS FOR KINDS 1 AND 2 C FROM CBESH AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION C IN THE (Z,FNU) SPACE. 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 CV, CW, CY, W, Y, Z, ZN REAL AA, ACW, ACY, AER, ALIM, ATOL, AV, AW, AY, AZ, C, C1, DIG, * ELIM, EPS, ER, ERTOL, FAC, FNU, FNUL, FPI, HPI, PI, R, RFPI, RL, * RM, R1M5, R2, S, SEPS, S1, T, TOL, TPI, XNU, R1MACH INTEGER I, IA, ICASE, IERR, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, * K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, N1, I1MACH DIMENSION C(20), S(20), AER(20), XNU(7), Y(20), W(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(FNUL,RM) C----------------------------------------------------------------------- WRITE (LUN,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM C, * 4HBESH/) 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 FAC = 100.0E0 ATOL = FAC*TOL FPI = ATAN(1.0E0) HPI = FPI + FPI PI = HPI + HPI TPI = PI + PI RFPI = 1.0E0/FPI ZN = CMPLX(0.0E0,-RFPI) 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 NEAR FORMULA BOUNDARIES 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 170 KODE=1,2 DO 160 N=1,4 N1 = N + 1 DO 150 NU=1,7 FNU = XNU(NU) DO 140 ICASE=1,3 IRB = MIN0(ICASE,2) DO 130 IR=IRB,3 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (EPS*FLOAT(3-IR)+2.0E0*FLOAT(IR-1))/2.0E0 GO TO 80 60 CONTINUE R = (2.0E0*FLOAT(3-IR)+R2*FLOAT(IR-1))/2.0E0 GO TO 80 70 CONTINUE IF (R2.EQ.RM) GO TO 140 R = (R2*FLOAT(3-IR)+RM*FLOAT(IR-1))/2.0E0 80 CONTINUE DO 120 IT=1,ITL Z = CMPLX(R*C(IT),R*S(IT)) IF (FNU.LT.2.0E0) GO TO 90 CV = Z*CMPLX(0.0E0,1.0E0) CALL CUOIK(CV, FNU, KODE, 2, N1, W, NZ2, TOL, ELIM, * ALIM) IF (NZ2.EQ.(-1)) GO TO 120 CV = -CV CALL CUOIK(CV, FNU, KODE, 2, N1, W, NZ2, TOL, ELIM, * ALIM) IF (NZ2.EQ.(-1)) GO TO 120 90 CONTINUE CALL CBESH(Z, FNU, KODE, 1, N1, Y, NZ1, IERR) IF (NZ1.NE.0) GO TO 120 CALL CBESH(Z, FNU, KODE, 2, N1, W, NZ2, IERR) IF (NZ2.NE.0) GO TO 120 CV = ZN/Z MFLG = 0 DO 100 I=1,N AW = CABS(W(I+1)) AY = CABS(Y(I)) AZ = ALOG(AW) + ALOG(AY) AZ = ABS(AZ) IF (AZ.GT.ALIM) GO TO 100 AV = CABS(CV) C ERROR RELATIVE TO MAXIMUM TERM CW = W(I)*Y(I+1) CY = W(I+1)*Y(I) CY = CW - CY - CV ACY = AW*AY ACW = CABS(W(I))*CABS(Y(I+1)) AV = AMAX1(ACW,ACY,AV) ER = CABS(CY)/AV AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 100 CONTINUE IF (MFLG.EQ.0) GO TO 120 IF (LFLG.EQ.1) GO TO 110 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), KK=IND, * 25HEX OF FIRST NON-ZERO PAIR/) LFLG = 1 110 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) 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE IF (LFLG.EQ.0) WRITE (LUN,99990) 99990 FORMAT (/16H QUICK CHECKS OK/) STOP END .