PROGRAM ZQCBK(INPUT,OUTPUT,TAPE7=OUTPUT) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBK IS A QUICK CHECK ROUTINE FOR THE COMPLEX K BESSEL FUNCTION C GENERATED BY SUBROUTINE ZBESK. C C ZQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM C ZBESI AND ZBESK 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 CONE,CSGN,CV,CW,CY,W,Y,Z,ZN DOUBLE PRECISION AA, AER, ALIM, ARG, ATOL, AXX, C, CONEI, CONER, * CSGNI, CSGNR, CVI, CVR, CWI, CWR, CYI, CYR, C1, DIG, ELIM, EPS, * ER, ERTOL, FAC, FFNU, FNU, FNUL, HPI, PI, R, RL, RM, R1M5, R2, * S, SEPS, STI, STR, S1, T, TOL, TPI, WI, WR, XNU, YI, YR, ZI, * ZNI, ZNR, ZR, D1MACH, ZABS INTEGER I, IA, ICASE, IFNU, IL, IR, IRB, IT, ITL, K, KK, KODE, * K1, K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, N1, IERR, I1MACH DIMENSION C(25), S(25), AER(25), 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 K BESSEL FUNCTION FROM ZB, * 3HESK/) 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 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 (/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 210 KODE=1,2 DO 200 N=1,4 N1 = N + 1 DO 190 NU=1,7 FNU = XNU(NU) IFNU = INT(SNGL(FNU)) FFNU = FNU - DBLE(FLOAT(IFNU)) ARG = PI*FFNU CSGNR = DCOS(ARG) CSGNI = DSIN(ARG) IF (MOD(IFNU,2).EQ.0) GO TO 50 CSGNR = -CSGNR CSGNI = -CSGNI 50 CONTINUE DO 180 ICASE=1,3 IRB = MIN0(2,ICASE) DO 170 IR=IRB,4 GO TO (60, 70, 80), ICASE 60 CONTINUE R = (EPS*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/ * 3.0D0 GO TO 90 70 CONTINUE R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0 GO TO 90 80 CONTINUE IF (R2.EQ.RM) GO TO 180 R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0 90 CONTINUE DO 160 IT=1,ITL ZR = R*C(IT) ZI = R*S(IT) CALL ZBESI(ZR, ZI, FNU, KODE, N1, WR, WI, NZ2, IERR) IF (NZ2.NE.0) GO TO 160 CALL ZBESK(ZR, ZI, FNU, KODE, N1, YR, YI, NZ1, IERR) IF (ICASE.EQ.1) GO TO 110 IF (IABS(IT-9).LE.4) GO TO 110 C----------------------------------------------------------------------- C IN THE LEFT HALF PLANE, THE ANALYTIC CONTINUATION FORMULA FOR K C INTRODUCES AN I FUNCTION. THE DOMINANT TERMS IN THE WRONSKIAN C I(FNU,Z)*I(FNU+1,Z) CANCEL OUT GIVING LOSSES OF SIGNIFICANCE. C THIS CANCELATION CAN BE DONE ANALYTICALLY TO GIVE A WRONSKIAN C IN TERMS OF I IN THE LEFT HALF PLANE AND K IN THE RIGHT HALF C PLANE. C----------------------------------------------------------------------- ZNR = -ZR ZNI = -ZI CALL ZBESK(ZNR, ZNI, FNU, KODE, N1, YR, YI, NZ1, IERR) ZNR = CSGNR ZNI = CSGNI IF (IT.GT.9) ZNI = -ZNI DO 100 KK=1,N1 STR = YR(KK)*ZNR - YI(KK)*ZNI YI(KK) = YR(KK)*ZNI + YI(KK)*ZNR YR(KK) = STR ZNR = -ZNR ZNI = -ZNI 100 CONTINUE 110 CONTINUE C----------------------------------------------------------------------- C ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS C ON KODE=2 C----------------------------------------------------------------------- CALL ZDIV(CONER, CONEI, ZR, ZI, CVR, CVI) IF (KODE.EQ.1) GO TO 130 AXX = DABS(ZR) ZNR = -AXX ZNI = 0.0D0 CVR = ZNR + ZR CVI = ZNI + ZI IF (ICASE.EQ.1) GO TO 120 IF (IABS(IT-9).LE.4) GO TO 120 CVR = ZNR - ZR CVI = ZNI - ZI 120 CONTINUE CALL ZEXP(CVR, CVI, STR, STI) CALL ZDIV(STR, STI, ZR, ZI, CVR, CVI) 130 CONTINUE MFLG = 0 DO 140 I=1,N 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 = CYR + CWR - CVR CYI = CYI + CWI - CVI ER = ZABS(CYR,CYI)/ZABS(CVR,CVI) AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 140 CONTINUE IF (MFLG.EQ.0) GO TO 160 IF (LFLG.EQ.1) GO TO 150 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 PUTPUT 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 150 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) 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE 210 CONTINUE IF (LFLG.EQ.0) WRITE (LUN,99990) 99990 FORMAT (/16H QUICK CHECKS OK/) STOP END .