PROGRAM ZQCBH(INPUT,OUTPUT,TAPE7=OUTPUT) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBH IS A QUICK CHECK ROUTINE FOR THE COMPLEX H BESSEL FUNCTIONS C GENERATED BY SUBROUTINE ZBESH. C C ZQCBH GENERATES SEQUENCES OF H BESSEL FUNCTIONS FOR KINDS 1 AND 2 C FROM ZBESH 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 C COMPLEX CV,CW,CY,W,Y,Z,ZN DOUBLE PRECISION AA, ACW, ACY, AER, ALIM, ATOL, AV, AW, AY, AZ, * C, CVI, CVR, CWI, CWR, CYI, CYR, C1, DIG, ELIM, EPS, ER, ERTOL, * FAC, FNU, FNUL, FPI, HPI, PI, R, RFPI, RL, RM, R1M5, R2, S, * SEPS, S1, T, TOL, TPI, WI, WR, XNU, YI, YR, ZI, ZNI, ZNR, ZR, * D1MACH, ZABS 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), YR(20), YI(20), WR(20), * WI(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 = DMAX1(D1MACH(4),1.0D-18) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 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) C----------------------------------------------------------------------- WRITE (LUN,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM Z, * 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 (6D12.4/) ERTOL = TOL*1.0D+5 FAC = 100.0D0 ATOL = FAC*TOL FPI = DATAN(1.0D0) HPI = FPI + FPI PI = HPI + HPI TPI = PI + PI RFPI = 1.0D0/FPI ZNR = 0.0D0 ZNI = -RFPI 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 (/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*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 XNU(1) = 0.0D0 XNU(2) = 0.3D0 XNU(3) = 0.7D0 XNU(4) = 1.0D0 XNU(5) = 1.3D0 XNU(6) = 1.7D0 XNU(7) = FNUL + 1.1D0 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*DBLE(FLOAT(3-IR))+2.0D0*DBLE(FLOAT(IR-1)))/ * 2.0D0 GO TO 80 60 CONTINUE R = (2.0D0*DBLE(FLOAT(3-IR))+R2*DBLE(FLOAT(IR-1)))/2.0D0 GO TO 80 70 CONTINUE IF (R2.EQ.RM) GO TO 140 R = (R2*DBLE(FLOAT(3-IR))+RM*DBLE(FLOAT(IR-1)))/2.0D0 80 CONTINUE DO 120 IT=1,ITL ZR = R*C(IT) ZI = R*S(IT) IF (FNU.LT.2.0D0) GO TO 90 CVR = -ZI CVI = ZR CALL ZUOIK(CVR, CVI, FNU, KODE, 2, N1, WR, WI, NZ2, * TOL, ELIM, ALIM) IF (NZ2.EQ.(-1)) GO TO 120 CVR = -CVR CVI = -CVI CALL ZUOIK(CVR, CVI, FNU, KODE, 2, N1, WR, WI, NZ2, * TOL, ELIM, ALIM) IF (NZ2.EQ.(-1)) GO TO 120 90 CONTINUE CALL ZBESH(ZR, ZI, FNU, KODE, 1, N1, YR, YI, NZ1,IERR) IF (NZ1.NE.0) GO TO 120 CALL ZBESH(ZR, ZI, FNU, KODE, 2, N1, WR, WI, NZ2,IERR) IF (NZ2.NE.0) GO TO 120 CALL ZDIV(ZNR, ZNI, ZR, ZI, CVR, CVI) MFLG = 0 DO 100 I=1,N AW = ZABS(WR(I+1),WI(I+1)) AY = ZABS(YR(I),YI(I)) AZ = DLOG(AW) + DLOG(AY) AZ = DABS(AZ) IF (AZ.GT.ALIM) GO TO 100 AV = ZABS(CVR,CVI) C ERROR RELATIVE TO MAXIMUM TERM CWR = WR(I)*YR(I+1) - WI(I)*YI(I+1) CWI = WR(I)*YI(I+1) + WI(I)*YR(I+1) CYR = WR(I+1)*YR(I) - WI(I+1)*YI(I) CYI = WR(I+1)*YI(I) + WI(I+1)*YR(I) CYR = CWR - CYR - CVR CYI = CWI - CYI - CVI ACY = AW*AY ACW = ZABS(WR(I),WI(I))*ZABS(YR(I+1),YI(I+1)) AV = DMAX1(ACW,ACY,AV) ER = ZABS(CYR,CYI)/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 =, D12.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) ZR, ZI, FNU, YR(KK), YI(KK), * WR(KK), WI(KK) 99991 FORMAT (7D12.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 .