program main c*********************************************************************72 c cc MAIN is the main program for TOMS644_PRB. c c Discussion: c c TOMS644_PRB tests the TOMS644 library. c c Modified: c c 13 February 2010 c implicit none write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS644_PRB:' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test the TOMS644 library.' c c Single precision c call cqcbh ( ) call cqcbi ( ) call cqcbj ( ) call cqcbk ( ) call cqcby ( ) call cqcai ( ) c c Double precision c call zqcbh ( ) call zqcbi ( ) call zqcbj ( ) call zqcbk ( ) call zqcby ( ) call zqcai ( ) c c Terminate. c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOMS644_PRB:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine CQCBH c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) 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 KIND 2 FROM C CBESH AND CHECKS THEM AGAINST ANALYTIC CONTINUATION FORMULAS C IN THE (Z,FNU) SPACE: C C KODE = 1 TESTS (ANALYTIC CONTINUATION FORMULAE, I**2 = -1): C C H(FNU,2,Z)=-EXP(I*PI*FNU)*H(FNU,1,-Z), -PI.LT.ARG(Z).LE.0 C C = 2*COS(PI*FNU)*H(FNU,2,-Z) + EXP(I*PI*FNU)*H(FNU,1,-Z), C C 0.LT.ARG(Z).LE.PI C C KODE = 2 TESTS FOR KINDS 1 AND 2: C C EXP(-I*Z)*H(FNU,1,Z) = [EXP(-I*Z)*H(FNU,1,Z)] C C EXP( I*Z)*H(FNU,2,Z) = [EXP( I*Z)*H(FNU,2,Z)] C C WHERE THE LEFT SIDE IS COMPUTED WITH KODE = 1 AND THE RIGHT SIDE C WITH KODE = 2. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CW, CI, U, V, W, Y, Z, ZN, CSGN REAL AA, AB, AER, ALIM, ATOL, AV, CT, DIG, ERR, XX, YY, * ELIM, EPS, ER, ERTOL, FNU, FNUL, PI, R, RL, * RM, R1M4, R1M5, R2, ST, T, TOL, TS, XNU, R1MACH, SLAK, FILM INTEGER I, ICASE, IERR, IHP, IL, IR, IRB, IT, ITL, K, KODE, KK, *K1, K2, LFLG, LUN, MFLG, M, N, NU, NZ1, NZ2, NZ3, I1MACH, KEPS, *MQC, NL, NUL, KDO DIMENSION T(20), AER(20), XNU(20), U(20), V(20), W(20), Y(20), *KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = AMAX1(R1M4,1.0E-18) AA = -ALOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + AMAX1(-AB,-41.45E0) DIG = AMIN1(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0 SLAK = AMAX1(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = AMIN1(RM,200.0E0) RM = AMAX1(RM,RL+10.0E0) R2 = AMIN1(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM CBES *H'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) ATOL = 100.0E0*TOL PI = 4.0E0*ATAN(1.0E0) CI = CMPLX(0.0E0,1.0E0) write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 170 KODE=1,2 DO 160 N=1,NL DO 150 NU=1,NUL 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.GE.RM) GO TO 140 R = (R2*FLOAT(3-IR)+RM*FLOAT(IR-1))/2.0E0 80 CONTINUE DO 120 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 XX = R*CT YY = R*ST Z = CMPLX(XX,YY) IF (KODE.EQ.1) THEN M=2 CALL CBESH(Z,FNU,KODE,M,N,Y,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 IF (ST.LT.0.0E0 .OR. (ST.EQ.0.0E0.AND.CT.GT.0.0E0)) * THEN IHP = 1 ZN=-Z M=1 CALL CBESH(ZN,FNU,KODE,M,N,W,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ELSE IHP = 2 ZN=-Z M=2 CALL CBESH(ZN,FNU,KODE,M,N,W,NZ3,IERR) IF (IERR.NE.0.OR.NZ3.NE.0) GO TO 120 M=1 CALL CBESH(ZN,FNU,KODE,M,N,V,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ENDIF AB=AMOD(FNU,2.0E0)*PI CSGN = CMPLX(COS(AB),SIN(AB)) MFLG = 0 DO 100 I=1,N AB = FNU+FLOAT(I-1) AA = MAX(0.5E0,AB) IF (IHP.EQ.1) THEN V(I) = -CSGN*W(I) CW = Y(I) - V(I) ELSE V(I) = 2.0E0*REAL(CSGN)*W(I) + CSGN*V(I) CW = Y(I) - V(I) ENDIF AV = CABS(Y(I)) ER = CABS(CW) IF (YY.EQ.0.0E0) THEN IF (ABS(XX).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 CSGN = -CSGN 100 CONTINUE ELSE M=1 KK=1 CALL CBESH(Z,FNU,KK,M,N,U,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 CALL CBESH(Z,FNU,KODE,M,N,V,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 M=2 KK=1 CALL CBESH(Z,FNU,KK,M,N,W,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 CALL CBESH(Z,FNU,KODE,M,N,Y,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ZN=CI*Z ZN=EXP(ZN) MFLG = 0 DO 105 I=1,N AB = FNU+FLOAT(I-1) AA = MAX(0.5E0,AB) CW = U(I)/ZN-V(I) ER = CABS(CW) AV = CABS(V(I)) IF (YY.EQ.0.0E0) THEN IF (ABS(XX).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF ERR = ER IF (ER.GT.ERTOL) MFLG = 1 CW=ZN*W(I)-Y(I) ER = CABS(CW) AV = CABS(Y(I)) IF (YY.EQ.0.0E0) THEN IF (ABS(XX).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF IF (ER.GT.ERTOL) MFLG = 1 AER(I) = ER+ERR 105 CONTINUE ENDIF IF (MFLG.EQ.0) GO TO 120 IF (LFLG.EQ.1) GO TO 110 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,V(1),Y(1)') LFLG = 1 110 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE 99992 FORMAT (5I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) Z, FNU, V(1), Y(1) 99991 FORMAT (7E12.4) 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine CQCBI c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C CQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION C GENERATED BY SUBROUTINE CBESI. C C CQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM C CBESI AND CBESK AND CHECKS THE WRONSKIAN EVALUATION C C I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z C C IN THE RIGHT HALF PLANE AND A MODIFIED FORM C C I(FNU+1,Z)*K(FNU,ZR) - I(FNU,Z)*K(FNU+1,ZR) = C/Z C C IN THE LEFT HALF PLANE WHERE ZR=-Z AND C=EXP(I*FNU*SGN) WITH C SGN=+1 FOR IM(Z).GE.0 AND SGN=-1 FOR IM(Z).LT.0. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CONE, CSGN, CC, CV, CW, CY, W, Y, Z, ZR REAL AA, AB, AER, ALIM, ARG, ATOL, DIG, ELIM, EPS, * ER, ERTOL, FFNU, FNU, FNUL, GNU, HPI, PI, R, RL, R1M4, R1M5, * R2, RM, T, TOL, XX, YY, R1MACH, SLAK, TS, ST, CT, FILM, XNU INTEGER I, ICASE, IFNU, IL, IPRNT, IR, IT, ITL, IRB, K, KK, KODE, * K1, K2, LFLG, LUN, MFLG, N, NU, NZ, N1, NUL, IERR, * MQC, NL, KEPS, KDO, I1MACH DIMENSION T(20), AER(20), Y(20), W(20), XNU(20), KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = MAX(R1M4,1.0E-18) AA = -ALOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + MAX(-AB,-41.45E0) DIG = MIN(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0 SLAK = MAX(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = MIN(RM,200.0E0) RM = MAX(RM,RL+10.0E0) R2 = MIN(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM CBESI *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) CONE = CMPLX(1.0E0,0.0E0) ATOL = 100.0E0*TOL HPI = 2.0E0*ATAN(1.0E0) PI = HPI + HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 200 KODE=1,2 DO 190 N=1,NL N1 = N + 1 DO 180 NU=1,NUL FNU = XNU(NU) IFNU = INT(FNU) FFNU = FNU - FLOAT(IFNU) ARG = PI*FFNU CSGN = CMPLX(COS(ARG),SIN(ARG)) IF (MOD(IFNU,2).EQ.1) CSGN = -CSGN DO 170 ICASE=1,3 IRB = MIN(2,ICASE) DO 160 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (0.2E0*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.GE.RM) GO TO 170 R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 80 CONTINUE DO 150 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 XX = R*CT YY = R*ST Z = CMPLX(XX,YY) IF (CT.GE.0.0E0) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE RIGHT HALF PLANE C----------------------------------------------------------------------- CALL CBESI(Z, FNU, KODE, N1, Y, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL CBESK(Z, FNU, KODE, N1, W, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 C----------------------------------------------------------------------- C ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS C ON KODE=2 C----------------------------------------------------------------------- CV = CONE/Z IF (KODE.EQ.2) THEN CW = CMPLX(COS(YY),SIN(YY)) CV = CW*CV ENDIF CC = CONE ELSE C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE LEFT HALF PLANE C----------------------------------------------------------------------- ZR = -Z CALL CBESI(Z, FNU, KODE, N1, Y, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL CBESK(ZR, FNU, KODE, N1, W, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CV = CSGN IF (YY.LT.0.0E0) THEN CV = CONJG(CV) ENDIF CV = CV/Z IF (KODE.EQ.2) THEN C----------------------------------------------------------------------- C ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS C ON KODE=2. SCALE FACTOR = EXP(-I*YY) FOR RE(Z).LT.0 C----------------------------------------------------------------------- CW = CMPLX(COS(YY),-SIN(YY)) CV = CV*CW ENDIF CC = -CONE ENDIF MFLG = 0 KK=0 DO 130 I=1,N CW = W(I)*Y(I+1) CY = CC*W(I+1)*Y(I) CY = CY + CW - CV ER = CABS(CY)/CABS(CV) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF IF (CT.LT.0.0E0) CV = -CV 130 CONTINUE IF (MFLG.EQ.0) GO TO 150 IF (LFLG.EQ.1) GO TO 140 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK) KK=INDEX O *F FIRST NON-ZERO PAIR'/) LFLG = 1 140 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE, KK 99992 FORMAT (6I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) Z, FNU, Y(KK) 99991 FORMAT (6E12.4) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) IF (MQC.EQ.1) then return end if C----------------------------------------------------------------------- C CHECKS NEAR UNDERFLOW LIMITS ON SERIES(I=1) AND UNIFORM C ASYMPTOTIC EXPANSION(I=2) C----------------------------------------------------------------------- write ( *,99989) 99989 FORMAT (/' CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS'/) Z = CMPLX(1.4E0,1.4E0) IPRNT = 0 DO 280 I=1,2 FNU = 10.2E0 KODE = 1 N = 20 230 CONTINUE CALL CBESI(Z, FNU, KODE, N, Y, NZ, IERR) IF (NZ.NE.0) GO TO 240 FNU = FNU + 5.0E0 GO TO 230 240 CONTINUE IF (NZ.LT.10) GO TO 250 FNU = FNU - 1.0E0 GO TO 230 250 CONTINUE CALL CBESK(Z, FNU, KODE, 2, W, NZ, IERR) CV = CONE/Z CY = W(1)*Y(2) CW = W(2)*Y(1) CW = CW + CY - CV ER = CABS(CW)/CABS(CV) IF (ER.LT.ERTOL) GO TO 270 IF (IPRNT.EQ.1) GO TO 260 write ( *,99988) 99988 FORMAT (/' OUTPUT FORMAT'/' ERROR,Z,FNU,KODE,N'/) IPRNT = 1 260 CONTINUE write ( *,99987) ER, Z, FNU, KODE, N 99987 FORMAT (4E12.4, 2I5) 270 CONTINUE XX = RL + RL Z = CMPLX(XX,0.0E0) 280 CONTINUE C----------------------------------------------------------------------- C CHECK NEAR OVERFLOW LIMITS C----------------------------------------------------------------------- Z = CMPLX(ELIM,0.0E0) FNU = 0.0 290 CONTINUE CALL CBESK(Z, FNU, KODE, N, Y, NZ, IERR) IF (NZ.LT.10) GO TO 300 IF (NZ.EQ.N) FNU = FNU + 3.0E0 FNU = FNU + 2.0E0 GO TO 290 300 CONTINUE GNU = FNU + FLOAT(N-2) CALL CBESI(Z, GNU, KODE, 2, W, NZ, IERR) CV = CONE/Z CY = Y(N-1)*W(2) CW = Y(N)*W(1) CW = CW + CY - CV ER = CABS(CW)/CABS(CV) IF (ER.LT.ERTOL) GO TO 320 IF (IPRNT.EQ.1) GO TO 310 write ( *,99988) IPRNT = 1 310 CONTINUE write ( *,99987) ER, Z, FNU, KODE, N 320 CONTINUE IF (IPRNT.EQ.0) write ( *,99990) return END subroutine CQCBJ c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) 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 AND H BESSEL FUNCTIONS FROM CBESJ C AND CBESH AND CHECKS THE WRONSKIANS C C J(FNU,Z)*H(FNU+1,1,Z)-J(FNU+1,Z)*H(FNU,1,Z)=2/(PI*I*Z) Y.GE.0 C C J(FNU,Z)*H(FNU+1,2,Z)-J(FNU+1,Z)*H(FNU,2,Z)=-2/(PI*I*Z) Y.LT.0 C C IN THEIR RESPECTIVE HALF PLANES. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 Z, WR, CJ, CH, CON, T1, T2, CER REAL AA, AB, AER, ALIM, DIG, ELIM, EPS, ER, ERTOL, FNU, FNUL, * GNU, HPI, R, RL, RM, R1M4, R1M5, R2, T, TOL, XNU, XX, YY, * R1MACH, SLAK, FILM, ST, TS, CT, SGN INTEGER I, ICASE, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, K2, * LFLG, LUN, M, MFLG, N, NU, NZJ, NZH, IERRJ, IERRH, I1MACH, KEPS, * KDO, NL, NUL, MQC DIMENSION T(20), AER(20), XNU(20), CJ(20), CH(20), KEPS(20), * KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = MAX(R1M4,1.0E-18) AA = -LOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + MAX(-AB,-41.45E0) DIG = MIN(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-LOG10(TOL)-7.0E0)/11.0E0 SLAK = MAX(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = MIN(RM,200.0E0) RM = MAX(RM,RL+10.0E0) R2 = MIN(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM CBESJ *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) ATOL=100.0E0*TOL HPI = 2.0E0*ATAN(1.0E0) PI = HPI + HPI CON=CMPLX(0.0,-1.0/HPI) write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 260 KODE=1,2 DO 250 N=1,NL NP=N+1 DO 240 NU=1,NUL FNU = XNU(NU) GNU = FNU + FLOAT(N-1) + 1.0E0 GNU=SQRT(GNU) GNU=MIN(GNU,0.5*RL) DO 230 ICASE=1,3 IRB = MIN(2,ICASE) DO 220 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (GNU*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.GE.RM) GO TO 230 R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 80 CONTINUE DO 210 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 Z = CMPLX(R*CT,R*ST) XX = REAL(Z) YY = AIMAG(Z) IF(XX.EQ.0.0.AND.YY.EQ.0.0) GO TO 210 WR=CON/Z M=1 IF(YY.LT.0.0) THEN M=2 WR=-WR ENDIF CALL CBESJ(Z,FNU,KODE,NP,CJ,NZJ,IERRJ) CALL CBESH(Z,FNU,KODE,M,NP,CH,NZH,IERRH) IF(NZJ.NE.0.OR.NZH.NE.0) GO TO 210 IF(IERRJ.NE.0.OR.IERRH.NE.0) GO TO 210 IF(KODE.EQ.2) THEN SGN=3.0-2.0*FLOAT(M) WR=WR*CMPLX(COS(XX),-SGN*SIN(XX)) ENDIF KK = 0 MFLG = 0 DO 190 I=1,N T1=CJ(I)*CH(I+1) T2=CJ(I+1)*CH(I) CER=T1-T2-WR ER=CABS(CER)/CABS(WR) IF (ER.GT.ERTOL) THEN IF(MFLG.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF AER(I)=ER 190 CONTINUE IF (MFLG.EQ.0) GO TO 210 IF (LFLG.EQ.1) GO TO 200 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZJ,NZH,ICASE *') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,CJ(KK),CH(KK), KK=INDEX * OF FIRST NON-ZERO CJ,CH PAIR'/) LFLG = 1 200 CONTINUE write ( *,99992) KODE, N, IR, IT, NZJ, NZH, ICASE 99992 FORMAT (8I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) Z, FNU, CJ(KK), CH(KK) 99991 FORMAT (9E12.4) 210 CONTINUE 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine CQCBK c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C CQCBK IS A QUICK CHECK ROUTINE FOR THE COMPLEX K BESSEL FUNCTION C GENERATED BY SUBROUTINE CBESK. C C CQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM C CBESI AND CBESK AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION C C I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z C C IN THE RIGHT HALF PLANE AND THE ANALYTIC CONTINUATION FORMULA C FOR H(FNU,2,Z) IN TERMS OF THE K FUNCTION C C K(FNU,Z) = C3*H(FNU,2,ZR) + C4*H(FNU,1,ZR) IM(Z).GE.0 C C = CONJG(K(FNU,CONJG(Z))) IM(Z).LT.0 C C IN THE LEFT HALF PLANE WHERE C3=C1*CONJG(C2)*C5, C4 = C2*C5 C C1=2*COS(PI*FNU), C2=EXP(PI*FNU*I/2), C5 =-PI*I/2 AND C ZR = Z*EXP(-3*PI*I/2) = Z*I C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CONE, CSGN, CV, CW, CY, C1, C2, C3, C4, V, W, Y, Z, ZR, * ZZ, CIP, COE REAL AA, AB, AER, ALIM, ARG, ATOL, DIG, ELIM, EPS, ER, * ERTOL, FFNU, FNU, FNUL, HPI, PI, R, RL, RM, R1M4, R1M5, R2, * T, TOL, XNU, XX, R1MACH, FILM, ST, CT, TS, SLAK, YY INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, I4, K, KK, KODE, * K1, K2, LFLG, LUN, M, MFLG, N, NU, NZ, N1, IERR, I1MACH, * KEPS, KDO, MQC, NL, NUL DIMENSION T(20), AER(20), XNU(20), V(20), Y(20), W(20), KEPS(20), * KDO(20), CIP(4) DATA CIP(1),CIP(2),CIP(3),CIP(4) / (1.0E0,0.0E0),(0.0E0,1.0E0), * (-1.0E0,0.0E0),(0.0E0,-1.0E0) / PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = MAX(R1M4,1.0E-18) AA = -ALOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + MAX(-AB,-41.45E0) DIG = MIN(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0 SLAK = MAX(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = MIN(RM,200.0E0) RM = MAX(RM,RL+10.0E0) R2 = MIN(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE K BESSEL FUNCTION FROM CBESK *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) CONE = CMPLX(1.0E0,0.0E0) ATOL = 100.0E0*TOL HPI = 2.0E0*ATAN(1.0E0) PI = HPI + HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 200 KODE=1,2 DO 190 N=1,NL N1 = N + 1 DO 180 NU=1,NUL FNU = XNU(NU) IFNU = INT(FNU) FFNU = FNU - FLOAT(IFNU) ARG = HPI*FFNU CSGN = CMPLX(COS(ARG),SIN(ARG)) I4 = MOD(IFNU,4)+1 CSGN = CSGN*CIP(I4) DO 170 ICASE=1,3 IRB = MIN(2,ICASE) DO 160 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (0.2E0*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.GE.RM) GO TO 170 R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 80 CONTINUE DO 150 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 XX = R*CT YY = R*ST Z = CMPLX(XX,YY) IF (XX.GE.0.0E0) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE RIGHT HALF PLANE C----------------------------------------------------------------------- CALL CBESI(Z, FNU, KODE, N1, W, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL CBESK(Z, FNU, KODE, N1, Y, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 C----------------------------------------------------------------------- C ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS C ON KODE=2 C----------------------------------------------------------------------- CV = CONE/Z IF (KODE.EQ.2) THEN CV = CV*CMPLX(COS(YY),SIN(YY)) ENDIF MFLG = 0 KK=0 DO 130 I=1,N CW = W(I)*Y(I+1) CY = W(I+1)*Y(I) CY = CY + CW - CV ER = CABS(CY)/CABS(CV) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF 130 CONTINUE ELSE C----------------------------------------------------------------------- C ANALYTIC CONTINUATION FORMULA CHECKS FOR LEFT HALF PLANE IN TERMS C OF H(FNU,1,Z) AND H(FNU,2,Z) C----------------------------------------------------------------------- ZZ = Z IF (YY.LT.0.0E0) THEN ZZ=CONJG(Z) ENDIF ZR = ZZ*CMPLX(0.0E0,1.0E0) M=1 CALL CBESH(ZR, FNU, KODE, M, N, W, NZ, IERR) IF (IERR.NE.0) GO TO 150 M=2 CALL CBESH(ZR, FNU, KODE, M, N, V, NZ, IERR) IF (IERR.NE.0) GO TO 150 CALL CBESK(Z, FNU, KODE, N, Y, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 COE=CMPLX(0.0E0,-HPI) MFLG = 0 KK = 0 AA = 2.0E0*COS(PI*FFNU) IF(MOD(IFNU,2).NE.0) AA = -AA C1 = CMPLX(AA,0.0E0) C2 = CSGN DO 135 I=1,N C3 = C1 C4 = C2 IF (KODE.EQ.2) THEN C----------------------------------------------------------------------- C ADJUSTMENTS TO COEFICIENTS DUE TO SCALING OF H(FNU,1,Z) AND C H(FNU,2,Z) FUNCTIONS ON KODE = 2. C----------------------------------------------------------------------- AB=CABS(V(I)) AA=ALOG(AB)+XX+XX IF (AA.GT.ELIM) GO TO 150 IF (AA.LT.-ELIM) THEN C3 = CMPLX(0.0E0,0.0E0) ELSE CW = ZZ+ZZ C3 = C3*CEXP(CW) ENDIF ENDIF CY = (C3*CONJG(C2)*V(I)+C4*W(I))*COE IF (YY.LT.0.0E0) THEN CY=CONJG(CY) ENDIF ER = CABS(CY-Y(I))/CABS(Y(I)) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF C2 = C2*CMPLX(0.0E0,1.0E0) C1 = -C1 135 CONTINUE ENDIF IF (MFLG.EQ.0) GO TO 150 IF (LFLG.EQ.1) GO TO 140 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK) KK=INDEX O *F FIRST NON-ZERO PAIR'/) LFLG = 1 140 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE, KK 99992 FORMAT (6I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) Z, FNU, Y(KK) 99991 FORMAT (6E12.4) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine CQCBY c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C CQCBY IS A QUICK CHECK ROUTINE FOR THE COMPLEX Y BESSEL FUNCTION C GENERATED BY SUBROUTINE CBESY. C C CQCBY GENERATES SEQUENCES OF Y BESSEL FUNCTIONS FROM CBESY AND C CBESYH AND COMPARES THEM FOR A VARIETY OF VALUES IN THE (Z,FNU) C SPACE. CBESYH IS AN OLD VERSION OF CBESY WHICH COMPUTES THE Y C FUNCTION FROM THE H FUNCTIONS OF KINDS 1 AND 2. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CW, CWRK, V, W, Z REAL AA, AB, AER, ALIM, ATOL, AV, CABS, DIG, ELIM, EPS, ER, * ERTOL, FFNU, FNU, FNUL, PI, R, RL, RM, R1M4, R1M5, R2, * T, TOL, XNU, XX, YY, R1MACH, CT, ST, TS, SLAK, FILM INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, K, KK, * KODE, K1, K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, IERR, I1MACH, * KEPS, KDO, NL, NUL, MQC DIMENSION T(20), AER(20), XNU(20), W(20), V(20), * CWRK(20), KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = AMAX1(R1M4,1.0E-18) AA = -ALOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + AMAX1(-AB,-41.45E0) DIG = AMIN1(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0 SLAK = AMAX1(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = AMIN1(RM,200.0E0) RM = AMAX1(RM,RL+10.0E0) R2 = AMIN1(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM CBESY * '/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) ATOL = 100.0E0*TOL PI = 4.0E0*ATAN(1.0E0) write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI/2.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 190 KODE=1,2 DO 180 N=1,NL DO 170 NU=1,NUL FNU = XNU(NU) IFNU = INT(FNU) FFNU = FNU-FLOAT(IFNU) DO 160 ICASE=1,3 IRB = MIN0(2,ICASE) DO 150 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 160 R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 80 CONTINUE DO 140 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 Z = CMPLX(R*CT,R*ST) XX = REAL(Z) YY = AIMAG(Z) CALL CBESY(Z, FNU, KODE, N, V, NZ2, CWRK, IERR) IF (NZ2.NE.0.OR.IERR.NE.0) GO TO 140 CALL CBESYH(Z, FNU, KODE, N, W, NZ1, CWRK, IERR) IF (NZ1.NE.0.OR.IERR.NE.0) GO TO 140 MFLG = 0 DO 120 I=1,N AB = FNU+FLOAT(I-1) AA = MAX(0.5E0,AB) CW = W(I) - V(I) AV = CABS(V(I)) ER = CABS(CW) IF (AV.NE.0.0E0) THEN IF (YY.EQ.0.0E0) THEN IF (XX.GT.0.0E0) THEN IF (ABS(XX).LT.AA) ER = ER/AV ELSE IF (ABS(FFNU-0.5E0).LT.0.125E0) THEN IF (ABS(XX).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF ENDIF ELSE ER = ER/AV ENDIF ENDIF AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 120 CONTINUE IF (MFLG.EQ.0) GO TO 140 IF (LFLG.EQ.1) GO TO 130 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL = ', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZ1,NZ2,ICASE *') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,W(KK),V(KK), KK=INDEX O *F FIRST NON-ZERO W,V PAIR'/) LFLG = 1 130 CONTINUE KK = MAX0(NZ1,NZ2) + 1 KK = MIN0(N,KK) write ( *,99992) KODE, N, IR, IT, NZ1, NZ2, ICASE 99992 FORMAT (8I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) Z, FNU, W(KK), V(KK) 99991 FORMAT (7E12.4) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END SUBROUTINE CBESYH(Z, FNU, KODE, N, CY, NZ, CWRK, IERR) c*********************************************************************72 c C***BEGIN PROLOGUE CBESYH C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, C BESSEL FUNCTION OF SECOND KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT C***DESCRIPTION C C ON KODE=1, CBESYH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESYH RETURNS THE SCALED C FUNCTIONS C C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) C C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS C (REF. 1). C C INPUT C Z - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(I)=Y(FNU+I-1,Z), I=1,...,N C = 2 RETURNS C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N C WHERE Y=AIMAG(Z) C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C CWRK - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N C C OUTPUT C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN C VALUES FOR THE SEQUENCE C CY(I)=Y(FNU+I-1,Z) OR C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N C DEPENDING ON KODE. C NZ - NZ=0 , A NORMAL RETURN C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO C UNDERFLOW (GENERALLY ON KODE=2) C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, FNU+N-1 IS C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT BY THE FORMULA C C Y(FNU,Z) = 0.5*(H(1,FNU,Z) - H(2,FNU,Z))/I C C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) C AND H(2,FNU,Z) ARE CALCULATED IN CBESH. C C FOR NEGATIVE ORDERS,THE FORMULA C C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) C C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. C MATH. SOFTWARE, 12, NO. 3, SEPTEMBER 1986, PP 265-273. C C***ROUTINES CALLED CBESH,I1MACH,R1MACH C***END PROLOGUE CBESYH C COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL, * ATOL, AA, BB, R1M5, TOL INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH DIMENSION CY(N), CWRK(N) C***FIRST EXECUTABLE STATEMENT CBESYH XX = REAL(Z) YY = AIMAG(Z) IERR = 0 NZ=0 IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1 IF (FNU.LT.0.0E0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN HCI = CMPLX(0.0E0,0.5E0) CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 NZ = MIN0(NZ1,NZ2) IF (KODE.EQ.2) GO TO 60 DO 50 I=1,N CY(I) = HCI*(CWRK(I)-CY(I)) 50 CONTINUE RETURN 60 CONTINUE TOL = AMAX1(R1MACH(4),1.0E-18) K1 = I1MACH(12) K2 = I1MACH(13) K = MIN0(IABS(K1),IABS(K2)) R1M5 = R1MACH(5) C----------------------------------------------------------------------- C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT C----------------------------------------------------------------------- ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) R1 = COS(XX) R2 = SIN(XX) EX = CMPLX(R1,R2) EY = 0.0E0 TAY = ABS(YY+YY) IF (TAY.LT.ELIM) EY = EXP(-TAY) IF (YY.LT.0.0E0) GO TO 90 C1 = EX*CMPLX(EY,0.0E0) C2 = CONJG(EX) 70 CONTINUE NZ = 0 RTOL = 1.0E0/TOL ASCLE = R1MACH(1)*RTOL*1.0E+3 DO 80 I=1,N C CY(I) = HCI*(C2*CWRK(I)-C1*CY(I)) ZV = CWRK(I) AA=REAL(ZV) BB=AIMAG(ZV) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75 ZV = ZV*CMPLX(RTOL,0.0E0) ATOL = TOL 75 CONTINUE ZV = ZV*C2*HCI ZV = ZV*CMPLX(ATOL,0.0E0) ZU=CY(I) AA=REAL(ZU) BB=AIMAG(ZU) ATOL=1.0E0 IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85 ZU = ZU*CMPLX(RTOL,0.0E0) ATOL = TOL 85 CONTINUE ZU = ZU*C1*HCI ZU = ZU*CMPLX(ATOL,0.0E0) CY(I) = ZV - ZU IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1 80 CONTINUE RETURN 90 CONTINUE C1 = EX C2 = CONJG(EX)*CMPLX(EY,0.0E0) GO TO 70 170 CONTINUE NZ = 0 RETURN END subroutine CQCAI c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C CQCAI IS A QUICK CHECK ROUTINE FOR THE COMPLEX AIRY FUNCTIONS C GENERATED BY SUBROUTINES CAIRY AND CBIRY. C C CQCAI GENERATES AIRY FUNCTIONS AND THEIR DERIVATIVES FROM CAIRY C AND CBIRY AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION IN THE C REGION -PI/3 .LE. ARG(Z) .LE. PI/3: C C AI(Z)*BI'(Z)-AI'(Z)*BI(Z)=1/PI. C C IN THE REMAINDER OF THE CUT PLANE, THE IDENTITIES C C AI(Z) = SQRT(-Z)*( J(-1/3,ZR) + J(1/3,ZR) )/3 C C AI'(Z) = Z*( J(-2/3,ZR) - J(2/3,ZR) )/3 C C BI(Z) = I*SQRT(-Z/3)*( C1*H(1/3,1,ZR) - C2*H(1/3,2,ZR) )/2 C C BI'(Z) = I*(-Z)/SQRT(3)*( C2*H(2/3,1,ZR) - C1*H(2/3,2,ZR) )/2 C C ARE CHECKED WHERE ZR = (2/3)(-Z)**(3/2) WITH C1 = EXP(PI*I/6), C C2 = CONJG(C1) AND I**2 = -1. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CA, CAV, CHI, CI, CONA, CONB, CONC, COND, CON1, CON2, CV, * CW, CY, W, Y, YY, Z, ZR, ZW, SC REAL AA, AB, ALIM, ATOL, AV, CT, C13, C23, C43, ZX, ZY, * DIG, ELIM, EPS, ER, ERTOL, FNUL, FPI, HPI, PI, R, RL, RPI, * R1M4, R1M5, SPI, ST, SLAK, T, TOL, RM, FILM, PI3, R1MACH, TS INTEGER I, ICASE, ICL, IL, IR, IRSET, IT, ITL, J, JB, JL, K, * KEPS, KODE, K1, K2, LFLG, LUN, NZ, IERR, I1MACH, MQC, IRB, KDO DIMENSION ER(5), T(20), Y(20), YY(20), W(20), KDO(20), KEPS(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = R1MACH(4) TOL = MAX(R1M4,1.0E-18) AA = -ALOG10(R1M4) K1 = I1MACH(12) K2 = I1MACH(13) R1M5 = R1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0) AB = AA*2.303E0 ALIM = ELIM + MAX(-AB,-41.45E0) DIG = MIN(AA,18.0E0) FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0) RL = 1.2E0*DIG + 3.0E0 SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0 SLAK = MAX(SLAK,3.0E0) ERTOL = TOL*10.0E0**SLAK RM = 0.5E0*(ALIM + ELIM) RM = MIN(RM,200.0E0) RM = MAX(RM,RL+10.0E0) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM CAIRY AN *D CBIRY'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6E12.4/) ATOL = 100.0E0*TOL FPI = ATAN(1.0E0) HPI = FPI + FPI PI = HPI + HPI RPI = 1.0E0/PI SPI = PI/6.0E0 CON1 = CMPLX(COS(SPI),SIN(SPI)) CON2 = CONJG(CON1) PI3 = SPI+SPI C13 = 1.0E0/3.0E0 C23 = C13+C13 C43 = C23+C23 AV = SQRT(C13) CAV = CMPLX(AV,0.0E0) CHI = CMPLX(0.0E0,0.5E0) CI = CMPLX(0.0E0,1.0E0) write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN ICL=1 IL=5 DO 5 I=1,IL KDO(I)=0 KEPS(I)=0 5 CONTINUE ELSE ICL=2 IL=7 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KEPS(2)=1 KEPS(3)=1 KEPS(5)=1 KEPS(6)=1 ENDIF I = 2 EPS = 0.01E0 FILM=FLOAT(IL-1) T(1) = -PI + EPS DO 30 K=2,IL IF(KDO(K).EQ.0) THEN T(I) = PI*FLOAT(-IL+2*K-1)/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 180 ICASE=1,ICL DO 170 KODE=1,2 DO 160 IRSET=1,3 IRB = MIN(IRSET,2) DO 150 IR=IRB,4 GO TO (40, 50, 60), IRSET 40 CONTINUE R = (0.2E0*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0 GO TO 70 50 CONTINUE R = (2.0E0*FLOAT(4-IR)+RL*FLOAT(IR-1))/3.0E0 GO TO 70 60 CONTINUE R = (RL*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0 70 CONTINUE DO 140 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0E0 IF (ABS(ST).LT.ATOL) ST = 0.0E0 ZX = R*CT ZY = R*ST Z = CMPLX(ZX,ZY) IF(ABS(T(IT)).LE.PI3) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECK IN -PI/3.LT.ARG(Z).LT.PI/3, TEST #1 C----------------------------------------------------------------------- CALL CAIRY(Z, 0, KODE, Y(1), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL CAIRY(Z, 1, KODE, Y(2), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL CBIRY(Z, 0, KODE, W(1), IERR) CALL CBIRY(Z, 1, KODE, W(2), IERR) CW = Y(1)*W(2) CY = Y(2)*W(1) CV = CMPLX(RPI,0.0E0) IF (KODE.EQ.2) THEN ZR=Z CA=CSQRT(ZR) ZR=ZR*CA*CMPLX(C23,0.0E0) AA=REAL(ZR) AA=ABS(AA) CA = ZR - CMPLX(AA,0.0E0) CV = CEXP(CA)*CV ENDIF CY = CW - CY - CV ER(1) = CABS(CY)/CABS(CV) JB = 1 JL = 1 ELSE C----------------------------------------------------------------------- C CHECKS IN -PI.LT.ARG(Z).LT.-PI/3 AND PI/3.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C CHECK AI TEST #2 C----------------------------------------------------------------------- CALL CAIRY(Z, 0, KODE, Y(2), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 ZR=-Z CV=CSQRT(ZR) ZR=ZR*CV*CMPLX(C23,0.0E0) CALL CBESJ(ZR,C23,KODE,2,YY,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CY=CMPLX(C43,0.0E0)*YY(1)/ZR - YY(2) CA = YY(1) CALL CBESJ(ZR,C13,KODE,2,YY,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 IF (KODE.EQ.2) THEN AB = AIMAG(ZR) AB = ABS(AB) CW = CSQRT(Z) ZW = Z*CW*CMPLX(C23,0.0E0) CW = ZW+CMPLX(AB,0.0E0) CW = CEXP(CW) YY(1) = YY(1)*CW YY(2) = YY(2)*CW CY = CY*CW CA = CA*CW SC = CW ENDIF CW = CV*CMPLX(C13,0.0E0) W(2) = CW*(YY(1)+CY) ER(2) = CABS(Y(2)-W(2)) IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN ER(2)=ER(2)/CABS(Y(2)) ELSE IF (KODE.EQ.2) THEN ER(2) = ER(2)/CABS(SC) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK AI' TEST #3 C----------------------------------------------------------------------- CY=CMPLX(C23,0.0E0)*YY(1)/ZR - YY(2) W(3) = Z*CMPLX(C13,0.0E0)*(CY-CA) CALL CAIRY(Z,1,KODE,Y(3),NZ,IERR) ER(3) = CABS(Y(3)-W(3)) IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN ER(3)=ER(3)/CABS(Y(3)) ELSE IF (KODE.EQ.2) THEN ER(3) = ER(3)/CABS(SC) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK BI TEST #4 C----------------------------------------------------------------------- CALL CBESH(ZR,C13,KODE,1,1,Y,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL CBESH(ZR,C13,KODE,2,1,YY,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CONA = CON1 CONB = CON2 CONC = CON2 COND = CON1 IF (KODE.EQ.2) THEN AA = REAL(ZW) AA = ABS(AA) ZW = CI*ZR - CMPLX(AA,0.0E0) CW = CEXP(ZW) CONA = CONA*CW CONC = CONC*CW ZW = -CI*ZR - CMPLX(AA,0.0E0) CW = CEXP(ZW) CONB = CONB*CW COND = COND*CW SC = CW ENDIF CW = CONA*Y(1)-CONB*YY(1) CW = CV*CAV*CW W(4) = CW*CHI CALL CBIRY(Z,0,KODE,Y(4),IERR) ER(4) = CABS(Y(4)-W(4)) IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN ER(4)=ER(4)/CABS(Y(4)) ELSE IF (KODE.EQ.2) THEN ER(4) = ER(4)/CABS(SC) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK BI' TEST #5 C----------------------------------------------------------------------- CALL CBESH(ZR,C23,KODE,1,1,Y,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL CBESH(ZR,C23,KODE,2,1,YY,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CW = CONC*Y(1)-COND*YY(1) CW = -Z*CAV*CW W(5) = CW*CHI CALL CBIRY(Z,1,KODE,Y(5),IERR) ER(5) = CABS(Y(5)-W(5)) IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN ER(5)=ER(5)/CABS(Y(5)) ELSE IF (KODE.EQ.2) THEN ER(5) = ER(5)/CABS(SC) ENDIF ENDIF JB = 2 JL = 5 ENDIF DO 190 J=JB,JL IF (ER(J).LT.ERTOL) GO TO 190 IF (LFLG.EQ.1) GO TO 130 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST W *ITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,IR,IT,IRSET,ICASE') write ( *,99993) 99993 FORMAT (' ER'/' I, Z, Y(I), W(I), ON THE I-TH TEST, I=1, *5'/) LFLG = 1 130 CONTINUE write ( *,99992) KODE, IR, IT, IRSET, ICASE 99992 FORMAT (5I5) write ( *,99991) ER(J) write ( *,99989) J, Z, Y(J), W(J) 99991 FORMAT (E12.4) 99989 FORMAT (I5,6E12.4) 190 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine ZQCBH c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) 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 KIND 2 FROM C ZBESH AND CHECKS THEM AGAINST ANALYTIC CONTINUATION FORMULAS C IN THE (Z,FNU) SPACE: C C KODE = 1 TESTS (ANALYTIC CONTINUATION FORMULAE, I**2 = -1): C C H(FNU,2,Z)=-EXP(I*PI*FNU)*H(FNU,1,-Z), -PI.LT.ARG(Z).LE.0 C C = 2*COS(PI*FNU)*H(FNU,2,-Z) + EXP(I*PI*FNU)*H(FNU,1,-Z), C C 0.LT.ARG(Z).LE.PI C C KODE = 2 TESTS FOR KINDS 1 AND 2: C C EXP(-I*Z)*H(FNU,1,Z) = [EXP(-I*Z)*H(FNU,1,Z)] C C EXP( I*Z)*H(FNU,2,Z) = [EXP( I*Z)*H(FNU,2,Z)] C C WHERE THE LEFT SIDE IS COMPUTED WITH KODE = 1 AND THE RIGHT SIDE C WITH KODE = 2. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CW, CI, U, V, W, Y, Z, ZN, CSGN EXTERNAL ZABS DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, AV, CT, DIG, ERR, ELIM, * EPS, ER, ERTOL, FNU, FNUL, PI, R, RL, RM, D1M4, D1M5, R2, ST, * T, TOL, TS, XNU, D1MACH, SLAK, FILM, STR, STI, UR, UI, VR, VI, * WR, WI, YR, YI, CWR, CWI, CSGNR, CSGNI, ZR, ZI, ZNR, ZNI, ZABS INTEGER I, ICASE, IERR, IHP, IL, IR, IRB, IT, ITL, K, KODE, KK, *K1, K2, LFLG, LUN, MFLG, M, N, NU, NZ1, NZ2, NZ3, I1MACH, KEPS, *MQC, NL, NUL, KDO DIMENSION T(20), AER(20), XNU(20), UR(20), UI(20), VR(20), VI(20), * WR(20), WI(20), YR(20), YI(20), KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- D1M4 = D1MACH(4) TOL = MAX(D1M4,1.0D-18) AA = -DLOG10(D1M4) K1 = I1MACH(15) K2 = I1MACH(16) D1M5 = D1MACH(5) K = MIN(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*D1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) R2 = MIN(FNUL,RM) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM ZBES *H'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) ATOL = 100.0D0*TOL PI = 4.0D0*DATAN(1.0D0) write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.1D0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0D0 XNU(2) = 0.6D0 XNU(3) = 1.3D0 XNU(4) = 2.0D0 XNU(5) = 0.5D0*FNUL XNU(6) = FNUL + 1.1D0 ENDIF I = 2 EPS = 0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 170 KODE=1,2 DO 160 N=1,NL DO 150 NU=1,NUL FNU = XNU(NU) DO 140 ICASE=1,3 IRB = MIN(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.GE.RM) GO TO 140 R = (R2*DBLE(FLOAT(3-IR))+RM*DBLE(FLOAT(IR-1)))/2.0D0 80 CONTINUE DO 120 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0D0 IF (ABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF (KODE.EQ.1) THEN M=2 CALL ZBESH(ZR,ZI,FNU,KODE,M,N,YR,YI,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 IF (ST.LT.0.0D0 .OR. (ST.EQ.0.0D0.AND.CT.GT.0.0D0)) * THEN IHP = 1 ZNR = -ZR ZNI = -ZI M=1 CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,WR,WI,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ELSE IHP = 2 ZNR = -ZR ZNI = -ZI M=2 CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,WR,WI,NZ3,IERR) IF (IERR.NE.0.OR.NZ3.NE.0) GO TO 120 M=1 CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,VR,VI,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ENDIF AB=MOD(FNU,2.0D0)*PI CSGNR = COS(AB) CSGNI = SIN(AB) MFLG = 0 DO 100 I=1,N AB = FNU+DBLE(FLOAT(I-1)) AA = MAX(0.5D0,AB) IF (IHP.EQ.1) THEN VR(I) = -(CSGNR*WR(I)-CSGNI*WI(I)) VI(I) = -(CSGNR*WI(I)+CSGNI*WR(I)) CWR = YR(I) - VR(I) CWI = YI(I) - VI(I) ELSE CWR = CSGNR+CSGNR STR = CWR*WR(I) + CSGNR*VR(I)-CSGNI*VI(I) VI(I) = CWR*WI(I) + CSGNR*VI(I)+CSGNI*VR(I) VR(I) = STR CWR = YR(I) - VR(I) CWI = YI(I) - VI(I) ENDIF AV = ZABS(YR(I),YI(I)) ER = ZABS(CWR,CWI) IF(ZI.EQ.0.0D0) THEN IF(ABS(ZR).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 CSGNR = -CSGNR CSGNI = -CSGNI 100 CONTINUE ELSE M=1 KK=1 CALL ZBESH(ZR,ZI,FNU,KK,M,N,UR,UI,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 CALL ZBESH(ZR,ZI,FNU,KODE,M,N,VR,VI,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 M=2 KK=1 CALL ZBESH(ZR,ZI,FNU,KK,M,N,WR,WI,NZ1,IERR) IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120 CALL ZBESH(ZR,ZI,FNU,KODE,M,N,YR,YI,NZ2,IERR) IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120 ZNR = -ZI ZNI = ZR CALL ZEXP(ZNR,ZNI,ZNR,ZNI) MFLG = 0 DO 105 I=1,N AB = FNU+DBLE(FLOAT(I-1)) AA = MAX(0.5D0,AB) CALL ZDIV(UR(I),UI(I),ZNR,ZNI,STR,STI) CWR = STR - VR(I) CWI = STI - VI(I) AV = ZABS(VR(I),VI(I)) ER = ZABS(CWR,CWI) IF(ZI.EQ.0.0D0) THEN IF(ABS(ZR).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF ERR = ER IF (ER.GT.ERTOL) MFLG = 1 CWR = ZNR*WR(I) - ZNI*WI(I) - YR(I) CWI = ZNR*WI(I) + ZNI*WR(I) - YI(I) AV = ZABS(YR(I),YI(I)) ER = ZABS(CWR,CWI) IF(ZI.EQ.0.0D0) THEN IF(ABS(ZR).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF IF (ER.GT.ERTOL) MFLG = 1 AER(I) = ER+ERR 105 CONTINUE ENDIF IF (MFLG.EQ.0) GO TO 120 IF (LFLG.EQ.1) GO TO 110 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', D12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,V(1),Y(1)') LFLG = 1 110 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE 99992 FORMAT (5I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) ZR,ZI,FNU,VR(1),VI(1),YR(1),YI(1) 99991 FORMAT (7D12.4) 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine ZQCBI c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION C GENERATED BY SUBROUTINE ZBESI. C C ZQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM C ZBESI AND CBESK AND CHECKS THE WRONSKIAN EVALUATION C C I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z C C IN THE RIGHT HALF PLANE AND A MODIFIED FORM C C I(FNU+1,Z)*K(FNU,ZR) - I(FNU,Z)*K(FNU+1,ZR) = C/Z C C IN THE LEFT HALF PLANE WHERE ZR=-Z AND C=EXP(I*FNU*SGN) WITH C SGN=+1 FOR IM(Z).GE.0 AND SGN=-1 FOR IM(Z).LT.0. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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, CC, CV, CW, CY, W, Y, Z, ZR EXTERNAL ZABS DOUBLE PRECISION AA, AB, AER, ALIM, ARG, ATOL, CCR, CONER, CONEI, * CSGNI, CSGNR, CWI, CWR, CYI, CYR, CVR, CVI, DIG, ELIM, EPS, ER, * ERTOL, FNU, FNUL, GNU, HPI, PI, R, RL, R1M4, R1M5, R2, RM, * STI, STR, T, TOL, WI, WR, YI, YR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS, * FILM, SLAK, TS, ST, CT, FFNU, XNU INTEGER I, ICASE, IL, IFNU, IPRNT, IR, IT, ITL, K, KK, KODE, K1, * K2, LFLG, LUN, MFLG, N, NZ, N1, NL, NU, NUL, MQC, IERR, I1MACH, * KEPS, KDO, IRB DIMENSION T(20), AER(20), YR(20), YI(20), WR(20), WI(20), XNU(20), * KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = MAX(R1M4,1.0D-18) AA = -DLOG10(R1M4) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) R2 = MIN(RM,FNUL) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM ZBESI *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) CONER = 1.0D0 CONEI = 0.0D0 ATOL = 100.0D0*TOL HPI = 2.0D0*DATAN(1.0D0) PI = HPI + HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0E0 XNU(2) = 1.0E0 XNU(3) = 2.0E0 XNU(4) = 0.5E0*FNUL XNU(5) = FNUL + 1.1E0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0E0 XNU(2) = 0.6E0 XNU(3) = 1.3E0 XNU(4) = 2.0E0 XNU(5) = 0.5E0*FNUL XNU(6) = FNUL + 1.1E0 ENDIF I = 2 EPS = 0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 200 KODE=1,2 DO 190 N=1,NL N1 = N + 1 DO 180 NU=1,NUL 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 170 ICASE=1,3 IRB = MIN(2,ICASE) DO 160 IR=IRB,4 GO TO (60, 70, 80), ICASE 60 CONTINUE R = (0.2D0*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 150 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (DABS(CT).LT.ATOL) CT = 0.0D0 IF (DABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF (CT.GE.0.0E0) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE RIGHT HALF PLANE C----------------------------------------------------------------------- CALL ZBESI(ZR, ZI, FNU, KODE, N1, WR, WI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL ZBESK(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 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.2) THEN STR = COS(ZI) STI = SIN(ZI) AA = STR*CVR - STI*CVI CVI = STR*CVI + STI*CVR CVR = AA ENDIF CCR = 1.0D0 ELSE C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE LEFT HALF PLANE C----------------------------------------------------------------------- ZRR = -ZR ZRI = -ZI CALL ZBESI(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL ZBESK(ZRR,ZRI, FNU, KODE, N1, WR, WI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CVR = CSGNR CVI = CSGNI IF (ZI.LT.0.0E0) THEN CVI = -CVI ENDIF CALL ZDIV(CVR,CVI,ZR,ZI,CVR,CVI) IF (KODE.EQ.2) THEN C----------------------------------------------------------------------- C ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS C ON KODE=2. SCALE FACTOR = EXP(-I*YY) FOR RE(Z).LT.0 C----------------------------------------------------------------------- CWR = COS(ZI) CWI = -SIN(ZI) STR = CVR*CWR-CVI*CWI CVI = CVR*CWI+CVI*CWR CVR = STR ENDIF CCR = -1.0D0 ENDIF MFLG = 0 KK = 0 DO 130 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 = YR(I)*WR(I+1)-YI(I)*WI(I+1) CYI = YR(I)*WI(I+1)+YI(I)*WR(I+1) CYR = CCR*CYR CYI = CCR*CYI CYR = CYR + CWR - CVR CYI = CYI + CWI - CVI ER = ZABS(CYR,CYI)/ZABS(CVR,CVI) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF IF (CT.LT.0.0E0) THEN CVR = -CVR CVI = -CVI ENDIF 130 CONTINUE IF (MFLG.EQ.0) GO TO 150 IF (LFLG.EQ.1) GO TO 140 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK) KK=INDEX O *F FIRST NON-ZERO PAIR'/) LFLG = 1 140 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE, KK 99992 FORMAT (6I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) ZR, ZI, FNU, YR(KK), YI(KK) 99991 FORMAT (6D12.4) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) IF (MQC.EQ.1) then return end if C----------------------------------------------------------------------- C CHECKS NEAR UNDERFLOW LIMITS ON SERIES(I=1) AND UNIFORM C ASYMPTOTIC EXPANSION(I=2) C----------------------------------------------------------------------- write ( *,99989) 99989 FORMAT (/' CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS'/) ZR = 1.4D0 ZI = 1.4D0 IPRNT = 0 DO 280 I=1,2 FNU = 10.2D0 KODE = 1 N = 20 230 CONTINUE CALL ZBESI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR) IF (NZ.NE.0) GO TO 240 FNU = FNU + 5.0D0 GO TO 230 240 CONTINUE IF (NZ.LT.10) GO TO 250 FNU = FNU - 1.0D0 GO TO 230 250 CONTINUE CALL ZBESK(ZR, ZI, FNU, KODE, 2, WR, WI, NZ, IERR) CALL ZDIV(CONER, CONEI, ZR, ZI, STR, STI) CYR = WR(1)*YR(2) - WI(1)*YI(2) CYI = WR(1)*YI(2) + WI(1)*YR(2) CWR = WR(2)*YR(1) - WI(2)*YI(1) CWI = WR(2)*YI(1) + WI(2)*YR(1) CWR = CWR + CYR - STR CWI = CWI + CYI - STI ER = ZABS(CWR,CWI)/ZABS(STR,STI) IF (ER.LT.ERTOL) GO TO 270 IF (IPRNT.EQ.1) GO TO 260 write ( *,99988) 99988 FORMAT (/' OUTPUT FORMAT/19H ERROR,Z,FNU,KODE,N'/) IPRNT = 1 260 CONTINUE write ( *,99987) ER, ZR, ZI, FNU, KODE, N 99987 FORMAT (4D12.4, 2I5) 270 CONTINUE ZR = RL +RL ZI = 0.0D0 280 CONTINUE C----------------------------------------------------------------------- C CHECK NEAR OVERFLOW LIMITS C----------------------------------------------------------------------- ZR = ELIM ZI = 0.0D0 FNU = 0.0D0 290 CONTINUE CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR) IF (NZ.LT.10) GO TO 300 IF (NZ.EQ.N) FNU = FNU + 3.0D0 FNU = FNU + 2.0D0 GO TO 290 300 CONTINUE GNU = FNU + DBLE(FLOAT(N-2)) CALL ZBESI(ZR, ZI, GNU, KODE, 2, WR, WI, NZ, IERR) CALL ZDIV(CONER, CONEI, ZR, ZI, STR, STI) CYR = YR(N-1)*WR(2) - YI(N-1)*WI(2) CYI = YR(N-1)*WI(2) + YI(N-1)*WR(2) CWR = YR(N)*WR(1) - YI(N)*WI(1) CWI = YR(N)*WI(1) + YI(N)*WR(1) CWR = CWR + CYR - STR CWI = CWI + CYI - STI ER = ZABS(CWR,CWI)/ZABS(STR,STI) IF (ER.LT.ERTOL) GO TO 320 IF (IPRNT.EQ.1) GO TO 310 write ( *,99988) IPRNT = 1 310 CONTINUE write ( *,99987) ER, ZR, ZI, FNU, KODE, N 320 CONTINUE IF (IPRNT.EQ.0) write ( *,99990) return END subroutine ZQCBJ c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBJ IS A QUICK CHECK ROUTINE FOR THE COMPLEX J BESSEL FUNCTION C GENERATED BY SUBROUTINE ZBESJ. C C ZQCBJ GENERATES SEQUENCES OF J AND H BESSEL FUNCTIONS FROM ZBESJ C AND ZBESH AND CHECKS THE WRONSKIANS C C J(FNU,Z)*H(FNU+1,1,Z)-J(FNU+1,Z)*H(FNU,1,Z)=2/(PI*I*Z) Y.GE.0 C C J(FNU,Z)*H(FNU+1,2,Z)-J(FNU+1,Z)*H(FNU,2,Z)=-2/(PI*I*Z) Y.LT.0 C C IN THEIR RESPECTIVE HALF PLANES. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 Z, WR, CJ, CH, CON, T1, T2, CER EXTERNAL ZABS DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, CONR, CONI, WRR, WRI, * CERR, CERI, CT, DIG, ELIM, EPS, ER, ERTOL, FNU, FNUL, GNU, HPI, * PI, R, RL, RM, R1M4, R1M5, R2, ST, STR, STI, TOL, T1R, T1I, T2R, * T2I, T, XNU, ZI, ZR, D1MACH, ZABS, FILM, SLAK, TS, SGN, CJR, CJI, * CHR, CHI INTEGER I, ICASE, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, K2, * LFLG, LUN, M, MFLG, N, NU, NZJ, NZH, I1MACH, KEPS, KDO, NL, *NUL, MQC, IERRJ, IERRH DIMENSION T(20), AER(25), XNU(20), CJR(20), CJI(20), CHR(20), * CHI(20), KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = DMAX1(R1M4,1.0D-18) AA = -DLOG10(R1M4) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN0(ABS(K1),ABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + DMAX1(-AB,-41.45D0) DIG = DMIN1(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = DMAX1(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = DMIN1(RM,200.0D0) RM = DMAX1(RM,RL+10.0D0) R2 = DMIN1(RM,FNUL) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM ZBESJ *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) ATOL = 100.0D0*TOL HPI = 2.0D0*DATAN(1.0D0) PI = HPI + HPI CONR = 0.0D0 CONI = -1.0D0/HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.1D0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0D0 XNU(2) = 0.6D0 XNU(3) = 1.3D0 XNU(4) = 2.0D0 XNU(5) = 0.5D0*FNUL XNU(6) = FNUL + 1.1D0 ENDIF I = 2 EPS = 0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 260 KODE=1,2 DO 250 N=1,NL NP = N+1 DO 240 NU=1,NUL FNU = XNU(NU) GNU = FNU + DBLE(FLOAT(N-1)) + 1.0D0 GNU = SQRT(GNU) GNU = MIN(GNU, 0.5D0*RL) DO 230 ICASE=1,3 IRB = MIN0(2,ICASE) DO 220 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (GNU*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/ * 3.0D0 GO TO 80 60 CONTINUE R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0 GO TO 80 70 CONTINUE IF (R2.GE.RM) GO TO 230 R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0 80 CONTINUE DO 210 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (ABS(CT).LT.ATOL) CT = 0.0D0 IF (ABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF(ZR.EQ.0.0.AND.ZI.EQ.0.0) GO TO 210 CALL ZDIV(CONR,CONI,ZR,ZI,WRR,WRI) M=1 IF(ZI.LT.0.0) THEN M=2 WRR = -WRR WRI = -WRI ENDIF CALL ZBESJ(ZR,ZI,FNU,KODE,NP,CJR,CJI,NZJ,IERRJ) CALL ZBESH(ZR,ZI,FNU,KODE,M,NP,CHR,CHI,NZH,IERRH) IF(NZJ.NE.0.OR.NZH.NE.0) GO TO 210 IF(IERRJ.NE.0.OR.IERRH.NE.0) GO TO 210 IF(KODE.EQ.2) THEN SGN = 3.0D0-2.0D0*DBLE(FLOAT(M)) STR = COS(ZR) STI = -SGN*SIN(ZR) T1R = WRR*STR - WRI*STI WRI = WRR*STI + WRI*STR WRR = T1R ENDIF KK = 0 MFLG = 0 DO 190 I=1,N STR = CJR(I)*CHR(I+1) - CJI(I)*CHI(I+1) T1I = CJR(I)*CHI(I+1) + CJI(I)*CHR(I+1) T1R = STR STR = CJR(I+1)*CHR(I) - CJI(I+1)*CHI(I) T2I = CJR(I+1)*CHI(I) + CJI(I+1)*CHR(I) T2R = STR CERR = T1R - T2R -WRR CERI = T1I - T2I -WRI ER=ZABS(CERR,CERI)/ZABS(WRR,WRI) IF (ER.GT.ERTOL) THEN IF(MFLG.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF AER(I)=ER 190 CONTINUE IF (MFLG.EQ.0) GO TO 210 IF (LFLG.EQ.1) GO TO 200 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', D12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZJ,NZH,ICASE *') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,CJ(KK),CH(KK), KK=INDEX * OF FIRST NON-ZERO Y,W PAIR'/) LFLG = 1 200 CONTINUE write ( *,99992) KODE, N, IR, IT, NZJ, NZH, ICASE 99992 FORMAT (8I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) ZR, ZI, FNU, CJR(KK), CJI(KK), * CHR(KK), CHI(KK) 99991 FORMAT (9D12.4) 210 CONTINUE 220 CONTINUE 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine ZQCBK c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) 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 C I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z C C IN THE RIGHT HALF PLANE AND THE ANALYTIC CONTINUATION FORMULA C FOR H(FNU,2,Z) IN TERMS OF THE K FUNCTION C C K(FNU,Z) = C3*H(FNU,2,ZR) + C4*H(FNU,1,ZR) IM(Z).GE.0 C C = CONJG(K(FNU,CONJG(Z))) IM(Z).LT.0 C C IN THE LEFT HALF PLANE WHERE C3=C1*CONJG(C2)*C5, C4 = C2*C5 C C1=2*COS(PI*FNU), C2=EXP(PI*FNU*I/2), C5 =-PI*I/2 AND C ZR = Z*EXP(-3*PI*I/2) = Z*I C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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, C1, C2, C3, C4, V, W, Y, Z, ZR, C * CIP EXTERNAL ZABS DOUBLE PRECISION AA, AB, AER, ALIM, ARG, ATOL, CONEI, CONER, * CSGNI, CSGNR, CVI, CVR, CWI, CWR, CYI, CYR, DIG, ELIM, EPS, * ER, ERTOL, FFNU, FNU, FNUL, HPI, PI, R, RL, RM, R1M4, R1M5, R2, * STI, STR, T, TOL, WI, WR, XNU, YI, YR, ZI, C1R, C1I, C2R, C2I, * ZRI, ZRR, ZR, D1MACH, ZABS, FILM, CT, ST, TS, SLAK, C3R, C3I, * C4R, C4I, VR, VI, CIPR, CIPI, COEI, ZZR, ZZI INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, I4, K, KK, KODE, * K1, K2, LFLG, LUN, M, MFLG, N, NU, N1, IERR, I1MACH, *KEPS, KDO, NL, NUL, NZ, MQC DIMENSION T(20), AER(25), XNU(20), VR(20), VI(20), YR(20), YI(20), * WR(20), WI(20), KEPS(20), KDO(20), CIPR(4), CIPI(4) DATA CIPR(1), CIPI(1), CIPR(2), CIPI(2), CIPR(3),CIPI(3),CIPR(4), * CIPI(4) / 1.0D0,0.0D0,0.0D0,1.0D0,-1.0D0,0.0D0,0.0D0,-1.0D0/ PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = MAX(R1M4,1.0D-18) AA = -DLOG10(R1M4) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) R2 = MIN(RM,FNUL) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE K BESSEL FUNCTION FROM ZB, * 3HESK/) write ( *,99998) 99998 FORMAT (37H PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG) write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) CONER = 1.0D0 CONEI = 0.0D0 ATOL = 100.0D0*TOL HPI = 2.0D0*DATAN(1.0D0) PI = HPI + HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE NUL=5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.1D0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(12)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 KEPS(10)=1 KEPS(11)=1 NUL=6 XNU(1) = 0.0D0 XNU(2) = 0.6D0 XNU(3) = 1.3D0 XNU(4) = 2.0D0 XNU(5) = 0.5D0*FNUL XNU(6) = FNUL + 1.1D0 ENDIF I = 2 EPS = 0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 200 KODE=1,2 DO 190 N=1,NL N1 = N + 1 DO 180 NU=1,NUL FNU = XNU(NU) IFNU = INT(SNGL(FNU)) FFNU = FNU - DBLE(FLOAT(IFNU)) ARG = HPI*FFNU CSGNR = DCOS(ARG) CSGNI = DSIN(ARG) I4 = MOD(IFNU,4)+1 STR = CSGNR*CIPR(I4) - CSGNI*CIPI(I4) CSGNI = CSGNR*CIPI(I4) + CSGNI*CIPR(I4) CSGNR = STR 50 CONTINUE DO 170 ICASE=1,3 IRB = MIN(2,ICASE) DO 160 IR=IRB,4 GO TO (60, 70, 80), ICASE 60 CONTINUE R = (0.2D0*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 150 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (DABS(CT).LT.ATOL) CT = 0.0D0 IF (DABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF (ZR.GE.0.0E0) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECKS IN THE RIGHT HALF PLANE C----------------------------------------------------------------------- CALL ZBESI(ZR, ZI, FNU, KODE, N1, WR, WI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 CALL ZBESK(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 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.2) THEN STR = COS(ZI) STI = SIN(ZI) AA = STR*CVR - STI*CVI CVI = STR*CVI + STI*CVR CVR = AA ENDIF MFLG = 0 KK = 0 DO 130 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 = YR(I)*WR(I+1)-YI(I)*WI(I+1) CYI = YR(I)*WI(I+1)+YI(I)*WR(I+1) CYR = CYR + CWR - CVR CYI = CYI + CWI - CVI ER = ZABS(CYR,CYI)/ZABS(CVR,CVI) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF 130 CONTINUE ELSE C----------------------------------------------------------------------- C ANALYTIC CONTINUATION FORMULA CHECKS FOR LEFT HALF PLANE IN TERMS C OF H(FNU,1,Z) AND H(FNU,2,Z) C----------------------------------------------------------------------- ZZR = ZR ZZI = ZI IF (ZI.LT.0.0E0) THEN ZZI = -ZZI ENDIF ZRR = -ZZI ZRI = ZZR M=1 CALL ZBESH(ZRR,ZRI,FNU,KODE,M,N,WR,WI,NZ,IERR) IF (IERR.NE.0) GO TO 150 M=2 CALL ZBESH(ZRR,ZRI,FNU,KODE,M,N,VR,VI,NZ,IERR) IF (IERR.NE.0) GO TO 150 CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150 COEI = -HPI MFLG = 0 KK = 0 AA = 2.0D0*COS(PI*FFNU) IF(MOD(IFNU,2).NE.0) AA = -AA C1R = AA C1I = 0.0D0 C2R = CSGNR C2I = CSGNI DO 135 I=1,N C3R = C1R C3I = C1I C4R = C2R C4I = C2I IF (KODE.EQ.2) THEN C----------------------------------------------------------------------- C ADJUSTMENTS TO COEFICIENTS DUE TO SCALING OF H(FNU,1,Z) AND C H(FNU,2,Z) FUNCTIONS ON KODE = 2. C----------------------------------------------------------------------- AB=ZABS(VR(I),VI(I)) AA=LOG(AB)+ZR+ZR IF (AA.GT.ELIM) GO TO 150 IF (AA.LT.-ELIM) THEN C3R = 0.0D0 C3I = 0.0D0 ELSE STR = ZZR+ZZR STI = ZZI+ZZI CALL ZEXP(STR,STI,CVR,CVI) STR = C3R*CVR - C3I*CVI C3I = C3R*CVI + C3I*CVR C3R = STR ENDIF ENDIF C CY = (C3*CONJG(C2)*V(I)+C4*W(I))*COE STR = C3R*C2R + C3I*C2I STI = -C3R*C2I + C3I*C2R CYR = STR*VR(I) - STI*VI(I) CYI = STR*VI(I) + STI*VR(I) CYR = CYR + (C4R*WR(I) - C4I*WI(I)) CYI = CYI + (C4R*WI(I) + C4I*WR(I)) STR = -CYI*COEI CYI = CYR*COEI CYR = STR IF (ZI.LT.0.0D0) THEN CYI = -CYI ENDIF STR = CYR - YR(I) STI = CYI - YI(I) ER = ZABS(STR,STI)/ZABS(YR(I),YI(I)) AER(I) = ER IF (ER.GT.ERTOL) THEN IF(KK.EQ.0) THEN MFLG = 1 KK=I ENDIF ENDIF STR = -C2I C2I = C2R C2R = STR C1R = -C1R C1I = -C1I 135 CONTINUE ENDIF IF (MFLG.EQ.0) GO TO 150 IF (LFLG.EQ.1) GO TO 140 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK) KK=INDEX O *F FIRST NON-ZERO PAIR'/) LFLG = 1 140 CONTINUE write ( *,99992) KODE, N, IR, IT, ICASE, KK 99992 FORMAT (6I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) ZR, ZI, FNU, YR(KK), YI(KK) 99991 FORMAT (6D12.4) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE 200 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END subroutine ZQCBY c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCBY IS A QUICK CHECK ROUTINE FOR THE COMPLEX Y BESSEL FUNCTION C GENERATED BY SUBROUTINE ZBESY. C C ZQCBY GENERATES SEQUENCES OF Y BESSEL FUNCTIONS FROM ZBESY AND C ZBESYH AND COMPARES THEM FOR A VARIETY OF VALUES IN THE (Z,FNU) C SPACE. ZBESYH IS AN OLD VERSION OF ZBESY WHICH COMPUTES THE Y C FUNCTION FROM THE H FUNCTIONS OF KINDS 1 AND 2. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CW,CWRK,V,W,Z EXTERNAL ZABS DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, AV, CWI, CWR, CWRKI, * CWRKR, DIG, ELIM, EPS, ER, ERTOL, FFNU, FNU, FNUL, HPI, PI, R, * RL, RM, R1M4, R1M5, R2, T, TOL, VI, VR, WI, WR, XNU, ZI, ZR, * D1MACH, ZABS, ST, CT, TS, SLAK, FILM INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, K, KK, * KODE, K1, K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, IERR, I1MACH, * MQC, KDO, KEPS, NL, NUL DIMENSION T(20), AER(20), XNU(20), WR(20), WI(20), VR(20), VI(20), * CWRKR(20), CWRKI(20), KDO(20), KEPS(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = MAX(R1M4,1.0D-18) AA = -DLOG10(R1M4) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM = MIN(RM,200.0D0) RM = MAX(RM,RL+10.0D0) R2 = MIN(RM,FNUL) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ZBESY *'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) ATOL = 100.0D0*TOL HPI = 2.0D0*DATAN(1.0D0) PI = HPI + HPI write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI/2.LT.ARG(Z).LE.PI NEAR BOUNDARIES C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN NL=2 IL=5 DO 5 I=1,IL KEPS(I)=0 KDO(I)=0 5 CONTINUE KDO(5)=1 NUL=5 XNU(1) = 0.0D0 XNU(2) = 1.0D0 XNU(3) = 2.0D0 XNU(4) = 0.5D0*FNUL XNU(5) = FNUL + 1.2D0 ELSE NL=4 IL=13 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KDO(2)=1 KDO(6)=1 KDO(8)=1 KDO(11)=1 KDO(12)=1 KDO(13)=1 KEPS(3)=1 KEPS(4)=1 KEPS(5)=1 KEPS(9)=1 NUL=6 XNU(1) = 0.0D0 XNU(2) = 0.6D0 XNU(3) = 1.3D0 XNU(4) = 2.0D0 XNU(5) = 0.5D0*FNUL XNU(6) = FNUL + 1.2D0 ENDIF I = 2 EPS = 0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF (KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 190 KODE=1,2 DO 180 N=1,NL DO 170 NU=1,NUL FNU = XNU(NU) IFNU = INT(SNGL(FNU)) FFNU = FNU - DBLE(FLOAT(IFNU)) DO 160 ICASE=1,3 IRB = MIN(2,ICASE) DO 150 IR=IRB,4 GO TO (50, 60, 70), ICASE 50 CONTINUE R = (EPS*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/ * 3.0D0 GO TO 80 60 CONTINUE R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0 GO TO 80 70 CONTINUE IF (RM.EQ.R2) GO TO 160 R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0 80 CONTINUE DO 140 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (DABS(CT).LT.ATOL) CT = 0.0D0 IF (DABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST CALL ZBESY(ZR, ZI, FNU, KODE, N, VR, VI, NZ2, CWRKR, & CWRKI, IERR) IF (NZ2.NE.0.OR.IERR.NE.0) GO TO 140 CALL ZBESYH(ZR, ZI, FNU, KODE, N, WR, WI, NZ1, CWRKR, & CWRKI, IERR) IF (NZ1.NE.0.OR.IERR.NE.0) GO TO 140 MFLG = 0 DO 120 I=1,N AB = FNU+DBLE(FLOAT(I-1)) AA = MAX(0.5D0,AB) CWR = WR(I) - VR(I) CWI = WI(I) - VI(I) AV = ZABS(VR(I),VI(I)) ER = ZABS(CWR,CWI) IF (AV.NE.0.0D0) THEN IF (ZI.EQ.0.0D0) THEN IF (ZR.GT.0.0D0) THEN IF (DABS(ZR).LT.AA) ER = ER/AV ELSE IF (DABS(FFNU-0.5D0).LT.0.125D0) THEN IF (DABS(ZR).LT.AA) ER = ER/AV ELSE ER = ER/AV ENDIF ENDIF ELSE ER = ER/AV ENDIF ENDIF AER(I) = ER IF (ER.GT.ERTOL) MFLG = 1 120 CONTINUE IF (MFLG.EQ.0) GO TO 140 IF (LFLG.EQ.1) GO TO 130 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST * WITH ERTOL = ', D12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZ1,NZ2,ICASE *') write ( *,99993) 99993 FORMAT (' ER(K),K=1,N'/' Z,FNU,W(KK),V(KK), KK=INDEX O *F FIRST NON-ZERO W,V PAIR'/) LFLG = 1 130 CONTINUE KK = MAX(NZ1,NZ2) + 1 KK = MIN(N,KK) write ( *,99992) KODE, N, IR, IT, NZ1, NZ2, ICASE 99992 FORMAT (8I5) write ( *,99991) (AER(K),K=1,N) write ( *,99991) ZR, ZI, FNU, WR(KK), WI(KK), * VR(KK), VI(KK) 99991 FORMAT (7D12.4) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END SUBROUTINE ZBESYH(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR, * CWRKI, IERR) c*********************************************************************72 c C***BEGIN PROLOGUE ZBESYH C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C***CATEGORY NO. B5K C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, C BESSEL FUNCTION OF SECOND KIND C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT C***DESCRIPTION C C ***A DOUBLE PRECISION ROUTINE*** C C ON KODE=1, ZBESYH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE C -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESYH RETURNS THE SCALED C FUNCTIONS C C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) C C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS C (REF. 1). C C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), C -PI.LT.ARG(Z).LE.PI C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION C KODE= 1 RETURNS C CY(I)=Y(FNU+I-1,Z), I=1,...,N C = 2 RETURNS C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N C WHERE Y=AIMAG(Z) C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 C CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT C CWRKI AT LEAST N C C OUTPUT CYR,CYI ARE DOUBLE PRECISION C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE C CY(I)=Y(FNU+I-1,Z) OR C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N C DEPENDING ON KODE. C NZ - NZ=0 , A NORMAL RETURN C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO C UNDERFLOW (GENERALLY ON KODE=2) C IERR - ERROR FLAG C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED C IERR=1, INPUT ERROR - NO COMPUTATION C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT C REDUCTION PRODUCE LESS THAN HALF OF MACHINE C ACCURACY C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- C CANCE BY ARGUMENT REDUCTION C IERR=5, ERROR - NO COMPUTATION, C ALGORITHM TERMINATION CONDITION NOT MET C C***LONG DESCRIPTION C C THE COMPUTATION IS CARRIED OUT BY THE FORMULA C C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I C C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) C AND H(2,FNU,Z) ARE CALCULATED IN ZBESH. C C FOR NEGATIVE ORDERS,THE FORMULA C C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU) C C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z). C C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. C C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, C OR -PI/2+P. C C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF C COMMERCE, 1955. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 C C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- C 1018, MAY, 1985 C C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. C MATH. SOFTWARE, 12, NO. 3, SEPTEMBER 1986, PP 265-273. C C***ROUTINES CALLED ZBESH,I1MACH,D1MACH C***END PROLOGUE ZBESYH C C COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R, * ELIM, EXI, EXR, EY, FNU, HCII, STI, STR, TAY, ZI, ZR, DEXP, * D1MACH, ASCLE, RTOL, ATOL, AA, BB, TOL, R1M5 INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH DIMENSION CYR(N), CYI(N), CWRKR(N), CWRKI(N) C***FIRST EXECUTABLE STATEMENT ZBESYH IERR = 0 NZ=0 IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1 IF (FNU.LT.0.0D0) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (N.LT.1) IERR=1 IF (IERR.NE.0) RETURN HCII = 0.5D0 CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, NZ1, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, NZ2, IERR) IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170 NZ = MIN0(NZ1,NZ2) IF (KODE.EQ.2) GO TO 60 DO 50 I=1,N STR = CWRKR(I) - CYR(I) STI = CWRKI(I) - CYI(I) CYR(I) = -STI*HCII CYI(I) = STR*HCII 50 CONTINUE RETURN 60 CONTINUE TOL = DMAX1(D1MACH(4),1.0D-18) K1 = I1MACH(15) K2 = I1MACH(16) K = MIN0(IABS(K1),IABS(K2)) R1M5 = D1MACH(5) C----------------------------------------------------------------------- C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT C----------------------------------------------------------------------- ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) EXR = DCOS(ZR) EXI = DSIN(ZR) EY = 0.0D0 TAY = DABS(ZI+ZI) IF (TAY.LT.ELIM) EY = DEXP(-TAY) IF (ZI.LT.0.0D0) GO TO 90 C1R = EXR*EY C1I = EXI*EY C2R = EXR C2I = -EXI 70 CONTINUE NZ = 0 RTOL = 1.0D0/TOL ASCLE = D1MACH(1)*RTOL*1.0D+3 DO 80 I=1,N C CY(I) = HCI*(C2*CWRK(I)-C1*CY(I)) C STR = C1R*CYR(I) - C1I*CYI(I) C STI = C1R*CYI(I) + C1I*CYR(I) C STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I) C STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I) C CYR(I) = -STI*HCII C CYI(I) = STR*HCII AA = CWRKR(I) BB = CWRKI(I) ATOL = 1.0D0 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 75 AA = AA*RTOL BB = BB*RTOL ATOL = TOL 75 CONTINUE STR = (AA*C2R - BB*C2I)*ATOL STI = (AA*C2I + BB*C2R)*ATOL AA = CYR(I) BB = CYI(I) ATOL = 1.0D0 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 85 AA = AA*RTOL BB = BB*RTOL ATOL = TOL 85 CONTINUE STR = STR - (AA*C1R - BB*C1I)*ATOL STI = STI - (AA*C1I + BB*C1R)*ATOL CYR(I) = -STI*HCII CYI(I) = STR*HCII IF (STR.EQ.0.0D0 .AND. STI.EQ.0.0D0 .AND. EY.EQ.0.0D0) NZ = NZ * + 1 80 CONTINUE RETURN 90 CONTINUE C1R = EXR C1I = EXI C2R = EXR*EY C2I = -EXI*EY GO TO 70 170 CONTINUE NZ = 0 RETURN END subroutine ZQCAI c*********************************************************************72 c C***DATE WRITTEN 830501 (YYMMDD) C***REVISION DATE 890801, 930101 (YYMMDD) C C *** A DOUBLE PRECISION ROUTINE *** C C ZQCAI IS A QUICK CHECK ROUTINE FOR THE COMPLEX AIRY FUNCTIONS C GENERATED BY SUBROUTINES ZAIRY AND ZBIRY. C C ZQCAI GENERATES AIRY FUNCTIONS AND THEIR DERIVATIVES FROM ZAIRY C AND ZBIRY AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION IN THE C REGION -PI/3 .LE. ARG(Z) .LE. PI/3: C C AI(Z)*BI'(Z)-AI'(Z)*BI(Z)=1/PI. C C IN THE REMAINDER OF THE CUT PLANE, THE IDENTITIES C C AI(Z) = SQRT(-Z)*( J(-1/3,ZR) + J(1/3,ZR) )/3 C C AI'(Z) = Z*( J(-2/3,ZR) - J(2/3,ZR) )/3 C C BI(Z) = I*SQRT(-Z/3)*( C1*H(1/3,1,ZR) - C2*H(1/3,2,ZR) )/2 C C BI'(Z) = I*(-Z)/SQRT(3)*( C2*H(2/3,1,ZR) - C1*H(2/3,2,ZR) )/2 C C ARE CHECKED WHERE ZR = (2/3)(-Z)**(3/2) WITH C1 = EXP(PI*I/6), C C2 = CONJG(C1) AND I**2 = -1. C C THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER, C LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST. 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 CA, CAV, CHI, CI, CONA, CONB, CONC, COND, CON1, CON2, CV, C * CW, CY, W, Y, YY, Z, ZR, ZW, SC EXTERNAL ZABS DOUBLE PRECISION AA, ALIM, ATOL, CAR, CAI, CVR, CVI, * CON1I, CON1R, CON2I, CON2R, CONAI, CONAR, CONBI, CONBR, CONCI, * CONCR, CONDI, CONDR, C13, C23, C43, CAVR, CAVI, CIR, CII, CHIR, * CHII, CWI, CWR, CYI, CYR, CT, DIG, ELIM, EPS, ER, ERTOL, FNUL, * FPI, HPI, PI, PTR, R, RL, RPI, R1M4, R1M5, SPI, STI, STR, ST, * T, TOL, RM, WI, WR, YI, YR, YYR, YYI, ZI, ZR, ZRI, ZRR, ZWR, ZWI, * D1MACH, ZABS, PI3, SLAK, FILM, AB, TS, SCR, SCI INTEGER I, ICASE, IL, IR, IRB, IRSET, IT, ITL, K, KEPS, KODE, J, * JB, JL, K1, K2, LFLG, LUN, NZ, IERR, I1MACH, MQC, ICL, KDO DIMENSION ER(5), T(20), YR(20), YI(20), YYR(20), YYI(20), WR(20), * WI(20), KEPS(20), KDO(20) PARAMETER (MQC=1) 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----------------------------------------------------------------------- R1M4 = D1MACH(4) TOL = MAX(R1M4,1.0D-18) AA = -DLOG10(R1M4) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) AB = AA*2.303D0 ALIM = ELIM + MAX(-AB,-41.45D0) DIG = MIN(AA,18.0D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0 SLAK = MAX(SLAK,3.0D0) ERTOL = TOL*10.0D0**SLAK RM = 0.5D0*(ALIM + ELIM) RM=MIN(RM,200.0D0) RM=MAX(RM,RL+10.0D0) C----------------------------------------------------------------------- write ( *,99999) 99999 FORMAT (' QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM ZAIRY AN *D ZBIRY'/) write ( *,99998) 99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG') write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) ATOL = 100.0D0*TOL FPI = DATAN(1.0D0) HPI = FPI + FPI PI = HPI + HPI RPI = 1.0D0/PI SPI = PI/6.0D0 CON1R = COS(SPI) CON1I = SIN(SPI) CON2R = CON1R CON2I = -CON1I PI3 = SPI+SPI C13 = 1.0D0/3.0D0 C23 = C13+C13 C43 = C23+C23 CAVR = SQRT(C13) CAVI = 0.0D0 CHIR = 0.0D0 CHII = 0.5D0 CIR = 0.0D0 CII = 1.0D0 write ( *,99996) MQC 99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/) C----------------------------------------------------------------------- C TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C KDO(K), K=1,IL DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI C ARE USE TO COMPUTE VALUES OF Z C KDO(K) = 0 MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO C VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K) C = 1 MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE C WILL BE SKIPPED C KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED C UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT C FORMULAE ARE USED. C KEPS(K) =0 MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE C =1 MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND C DOWN BY EPS C THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL C----------------------------------------------------------------------- IF (MQC.NE.2) THEN ICL=1 IL=5 DO 5 I=1,IL KDO(I)=0 KEPS(I)=0 5 CONTINUE ELSE ICL=2 IL=7 DO 6 I=1,IL KDO(I)=0 KEPS(I)=0 6 CONTINUE KEPS(2)=1 KEPS(3)=1 KEPS(5)=1 KEPS(6)=1 ENDIF I = 2 EPS=0.01D0 FILM=DBLE(FLOAT(IL-1)) T(1) = -PI + EPS DO 30 K=2,IL IF(KDO(K).EQ.0) THEN T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM IF (KEPS(K).EQ.0) GO TO 20 TS=T(I) T(I) = TS - EPS I = I + 1 T(I) = TS + EPS 20 CONTINUE I = I + 1 ENDIF 30 CONTINUE ITL = I - 1 LFLG = 0 DO 180 ICASE=1,ICL DO 170 KODE=1,2 DO 160 IRSET=1,3 IRB = MIN(IRSET,2) DO 150 IR=IRB,4 GO TO (40, 50, 60), IRSET 40 CONTINUE R =(0.2D0*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/3.0D0 GO TO 70 50 CONTINUE R = (2.0D0*DBLE(FLOAT(4-IR))+RL*DBLE(FLOAT(IR-1)))/3.0D0 GO TO 70 60 CONTINUE R = (RL*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0 70 CONTINUE DO 140 IT=1,ITL CT = COS(T(IT)) ST = SIN(T(IT)) IF (DABS(CT).LT.ATOL) CT = 0.0D0 IF (DABS(ST).LT.ATOL) ST = 0.0D0 ZR = R*CT ZI = R*ST IF(ABS(T(IT)).LE.PI3) THEN C----------------------------------------------------------------------- C WRONSKIAN CHECK IN -PI/3.LT.ARG(Z).LT.PI/3, TEST #1 C----------------------------------------------------------------------- CALL ZAIRY(ZR, ZI, 0, KODE, YR(1), YI(1), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL ZAIRY(ZR, ZI, 1, KODE, YR(2), YI(2), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL ZBIRY(ZR, ZI, 0, KODE, WR(1), WI(1), IERR) CALL ZBIRY(ZR, ZI, 1, KODE, WR(2), WI(2), IERR) CWR = YR(1)*WR(2)-YI(1)*WI(2) CWI = YR(1)*WI(2)+YI(1)*WR(2) CYR = YR(2)*WR(1)-YI(2)*WI(1) CYI = YR(2)*WI(1)+YI(2)*WR(1) CVR = RPI CVI = 0.0D0 IF (KODE.EQ.2) THEN CALL ZSQRT(ZR,ZI,CAR,CAI) ZRR = (ZR*CAR-ZI*CAI)*C23 ZRI = (ZR*CAI+ZI*CAR)*C23 AA=ABS(ZRR) CAR = ZRR - AA CAI = ZRI CALL ZEXP(CAR,CAI,STR,STI) PTR = STR*CVR-STR*CVI CVI = STR*CVI+STI*CVR CVR = PTR ENDIF CYR = CWR - CYR - CVR CYI = CWI - CYI - CVI ER(1) = ZABS(CYR,CYI)/ZABS(CVR,CVI) JB = 1 JL = 1 ELSE C----------------------------------------------------------------------- C CHECKS IN -PI.LT.ARG(Z).LT.-PI/3 AND PI/3.LT.ARG(Z).LE.PI C----------------------------------------------------------------------- C----------------------------------------------------------------------- C CHECK AI TEST #2 C----------------------------------------------------------------------- CALL ZAIRY(ZR, ZI, 0, KODE, YR(2), YI(2), NZ, IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 ZRR = -ZR ZRI = -ZI CALL ZSQRT(ZRR,ZRI,CVR,CVI) PTR = (ZRR*CVR-ZRI*CVI)*C23 ZRI = (ZRR*CVI+ZRI*CVR)*C23 ZRR = PTR CALL ZBESJ(ZRR,ZRI,C23,KODE,2,YYR,YYI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 STR = YYR(1)*C43 STI = YYI(1)*C43 CALL ZDIV(STR,STI,ZRR,ZRI,STR,STI) CYR = STR - YYR(2) CYI = STI - YYI(2) CAR = YYR(1) CAI = YYI(1) CALL ZBESJ(ZRR,ZRI,C13,KODE,2,YYR,YYI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 IF (KODE.EQ.2) THEN AB = ABS(ZRI) CALL ZSQRT(ZR,ZI,CWR,CWI) ZWR = (ZR*CWR-ZI*CWI)*C23 ZWI = (ZR*CWI+ZI*CWR)*C23 CWR = ZWR+AB CWI = ZWI CALL ZEXP(CWR,CWI,STR,STI) CWR = STR CWI = STI STR = YYR(1)*CWR-YYI(1)*CWI YYI(1) = YYR(1)*CWI+YYI(1)*CWR YYR(1) = STR STR = YYR(2)*CWR-YYI(2)*CWI YYI(2) = YYR(2)*CWI+YYI(2)*CWR YYR(2) = STR STR = CYR*CWR-CYI*CWI CYI = CYR*CWI+CYI*CWR CYR = STR STR = CAR*CWR-CAI*CWI CAI = CAR*CWI+CAI*CWR CAR = STR SCR = CWR SCI = CWI ENDIF CWR = CVR*C13 CWI = CVI*C13 WR(2) = CWR*(YYR(1)+CYR)-CWI*(YYI(1)+CYI) WI(2) = CWR*(YYI(1)+CYI)+CWI*(YYR(1)+CYR) STR = YR(2)-WR(2) STI = YI(2)-WI(2) ER(2) = ZABS(STR,STI) IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN ER(2) = ER(2)/ZABS(YR(2),YI(2)) ELSE IF (KODE.EQ.2) THEN ER(2) = ER(2)/ZABS(SCR,SCI) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK AI' TEST #3 C----------------------------------------------------------------------- STR = YYR(1)*C23 STI = YYI(1)*C23 CALL ZDIV(STR,STI,ZRR,ZRI,CYR,CYI) CYR = CYR-YYR(2)-CAR CYI = CYI-YYI(2)-CAI WR(3) = (ZR*CYR-ZI*CYI)*C13 WI(3) = (ZR*CYI+ZI*CYR)*C13 CALL ZAIRY(ZR, ZI, 1, KODE, YR(3), YI(3), NZ, IERR) STR = YR(3)-WR(3) STI = YI(3)-WI(3) ER(3) = ZABS(STR,STI) IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN ER(3) = ER(3)/ZABS(YR(3),YI(3)) ELSE IF (KODE.EQ.2) THEN ER(3) = ER(3)/ZABS(SCR,SCI) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK BI TEST #4 C----------------------------------------------------------------------- CALL ZBESH(ZRR,ZRI,C13,KODE,1,1,YR,YI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL ZBESH(ZRR,ZRI,C13,KODE,2,1,YYR,YYI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CONAR = CON1R CONAI = CON1I CONBR = CON2R CONBI = CON2I CONCR = CON2R CONCI = CON2I CONDR = CON1R CONDI = CON1I IF (KODE.EQ.2) THEN AA = ABS(ZWR) ZWR = CIR*ZRR-CII*ZRI-AA ZWI = CIR*ZRI+CII*ZRR CALL ZEXP(ZWR,ZWI,CWR,CWI) STR = CONAR*CWR-CONAI*CWI CONAI =CONAR*CWI+CONAI*CWR CONAR = STR STR = CONCR*CWR-CONCI*CWI CONCI =CONCR*CWI+CONCI*CWR CONCR = STR ZWR = -(CIR*ZRR-CII*ZRI)-AA ZWI = -(CIR*ZRI+CII*ZRR) CALL ZEXP(ZWR,ZWI,CWR,CWI) STR = CONBR*CWR-CONBI*CWI CONBI =CONBR*CWI+CONBI*CWR CONBR = STR STR = CONDR*CWR-CONDI*CWI CONDI =CONDR*CWI+CONDI*CWR CONDR = STR SCR = CWR SCI = CWI ENDIF CWR = CONAR*YR(1)-CONAI*YI(1) CWI = CONAR*YI(1)+CONAI*YR(1) CWR = CWR - (CONBR*YYR(1)-CONBI*YYI(1)) CWI = CWI - (CONBR*YYI(1)+CONBI*YYR(1)) STR = CVR*CAVR-CVI*CAVI STI = CVR*CAVI+CVI*CAVR PTR = STR*CWR-STI*CWI CWI = STR*CWI+STI*CWR CWR = PTR WR(4) = CWR*CHIR-CWI*CHII WI(4) = CWR*CHII+CWI*CHIR CALL ZBIRY(ZR,ZI,0,KODE,YR(4),YI(4),IERR) STR = YR(4)-WR(4) STI = YI(4)-WI(4) ER(4) = ZABS(STR,STI) IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN ER(4) = ER(4)/ZABS(YR(4),YI(4)) ELSE IF (KODE.EQ.2) THEN ER(4) = ER(4)/ZABS(SCR,SCI) ENDIF ENDIF C----------------------------------------------------------------------- C CHECK BI' TEST #5 C----------------------------------------------------------------------- CALL ZBESH(ZRR,ZRI,C23,KODE,1,1,YR,YI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CALL ZBESH(ZRR,ZRI,C23,KODE,2,1,YYR,YYI,NZ,IERR) IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140 CWR = CONCR*YR(1)-CONCI*YI(1) CWI = CONCR*YI(1)+CONCI*YR(1) CWR = CWR - (CONDR*YYR(1)-CONDI*YYI(1)) CWI = CWI - (CONDR*YYI(1)+CONDI*YYR(1)) STR = -(ZR*CAVR-ZI*CAVI) STI = -(ZR*CAVI+ZI*CAVR) PTR = STR*CWR-STI*CWI CWI = STR*CWI+STI*CWR CWR = PTR WR(5) = CWR*CHIR-CWI*CHII WI(5) = CWR*CHII+CWI*CHIR CALL ZBIRY(ZR,ZI,1,KODE,YR(5),YI(5),IERR) STR = YR(5)-WR(5) STI = YI(5)-WI(5) ER(5) = ZABS(STR,STI) IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN ER(5) = ER(5)/ZABS(YR(5),YI(5)) ELSE IF (KODE.EQ.2) THEN ER(5) = ER(5)/ZABS(SCR,SCI) ENDIF ENDIF JB = 2 JL = 5 ENDIF DO 190 J=JB,JL IF (ER(J).LT.ERTOL) GO TO 190 IF (LFLG.EQ.1) GO TO 130 write ( *,99995) ERTOL 99995 FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST W *ITH ERTOL =', E12.4/) write ( *,99994) 99994 FORMAT (/' OUTPUT FORMAT'/' KODE,IR,IT,IRSET,ICASE') write ( *,99993) 99993 FORMAT (' ER'/' I, Z, Y(I), W(I), ON THE I-TH TEST, I=1, *5'/) LFLG = 1 130 CONTINUE write ( *,99992) KODE, IR, IT, IRSET, ICASE 99992 FORMAT (5I5) write ( *,99991) ER(J) write ( *,99989) J, ZR, ZI, YR(J), YI(J), WR(J), WI(J) 99991 FORMAT (D12.4) 99989 FORMAT (I5,6D12.4) 190 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE IF (LFLG.EQ.0) write ( *,99990) 99990 FORMAT (/' QUICK CHECKS OK'/) return END