PROGRAM ZQCAI(INPUT,OUTPUT,TAPE7=OUTPUT) 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 Z PLANE. 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 CON1,CON2,CON3,CV,CW,CY,W,Y,Z,ZR DOUBLE PRECISION AA, ACW, ACY, ALIM, ARG, ATOL, AV, AZRR, A1, A2, * C, CON1I, CON1R, CON2I, CON2R, CON3I, CON3R, CVI, CVR, CWI, CWR, * CYI, CYR, C1, C23, DIG, ELIM, EPS, ER, ERTOL, FAC, FNUL, FPI, * HPI, PI, PTR, R, RL, RPI, RTPI, R1M5, S, SEPS, SPI, STI, STR, * S1, T, TOL, TPI, TPI3, RM, WI, WR, YI, YR, ZI, ZR, ZRI, ZRR, * D1MACH, ZABS INTEGER I, IA, ICASE, IL, IR, IRB, IRSET, IT, ITL, K, KODE, * K1, K2, LFLG, LUN, NZ1, NZ2, NZ3, NZ4, IERR, I1MACH DIMENSION C(20), S(20), YR(20), YI(20), WR(20), WI(20) DATA LUN /7/ C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU C----------------------------------------------------------------------- TOL = DMAX1(D1MACH(4),1.0D-18) K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN0(IABS(K1),IABS(K2)) ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) K1 = I1MACH(14) - 1 AA = R1M5*DBLE(FLOAT(K1)) DIG = DMIN1(AA,18.0D0) AA = AA*2.303D0 ALIM = ELIM + DMAX1(-AA,-41.45D0) FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) RL = 1.2D0*DIG + 3.0D0 C----------------------------------------------------------------------- WRITE (LUN,99999) 99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM ZAIRY, * 10H AND ZBIRY/) WRITE (LUN,99998) 99998 FORMAT (37H PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG) WRITE (LUN,99997) TOL, ELIM, ALIM, RL, FNUL, DIG 99997 FORMAT (6D12.4/) ERTOL = TOL*1.0D+5 FAC = 100.0D0 RM = 0.5D0*(ALIM+ELIM) RM=DMIN1(RM,650.0D0) ATOL = FAC*TOL FPI = DATAN(1.0D0) HPI = FPI + FPI PI = HPI + HPI TPI = PI + PI RPI = 1.0D0/PI TPI3 = TPI/3.0D0 CON1R = DCOS(TPI3) CON1I = DSIN(TPI3) SPI = PI/6.0D0 RTPI = 1.0D0/TPI A1 = RTPI*DCOS(SPI) A2 = RTPI*DSIN(SPI) CON2R = A1 CON2I = -A2 CON3R = RPI CON3I = 0.0D0 C23 = 2.0D0/3.0D0 I = 1 IA = 3 EPS = 0.01D0 SEPS = EPS IL = 7 C----------------------------------------------------------------------- C TEST 9 VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES C----------------------------------------------------------------------- WRITE (LUN,99996) 99996 FORMAT (/22H CHECKS IN THE Z PLANE/) DO 30 K=1,IL IF (K.NE.4) GO TO 10 IA = 5 SEPS = -SEPS 10 CONTINUE T = -PI + TPI*DBLE(FLOAT(K-1))/DBLE(FLOAT(IL-1)) C1 = DCOS(T) S1 = DSIN(T) IF (DABS(C1).LT.ATOL) C1 = 0.0D0 IF (DABS(S1).LT.ATOL) S1 = 0.0D0 C(I) = C1 S(I) = S1 IF (IABS(K-IA).GT.0) GO TO 20 C(I) = C1 - SEPS S(I) = S1 C(I+1) = C1 + SEPS S(I+1) = S1 I = I + 2 GO TO 30 20 CONTINUE I = I + 1 30 CONTINUE S(1) = -EPS ITL = I - 1 LFLG = 0 DO 180 ICASE=1,2 C----------------------------------------------------------------------- C ICASE=1 COMPUTES WRON(AI(Z),BI(Z)) =CON3 C ICASE=2 COMPUTES WRON(AI(Z),AI(Z*CON1))=CON2 C----------------------------------------------------------------------- DO 170 KODE=1,2 DO 160 IRSET=1,3 IRB = MIN0(IRSET,2) DO 150 IR=IRB,4 GO TO (40, 50, 60), IRSET 40 CONTINUE R = 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 ZR = R*C(IT) ZI = R*S(IT) CALL ZSQRT(ZR, ZI, STR, STI) PTR = (ZR*STR-ZI*STI)*C23 ZRI = (ZR*STI+ZI*STR)*C23 ZRR = PTR AZRR = DABS(ZRR) IF (AZRR.EQ.0.0D0) GO TO 80 ARG = -AZRR - 0.5D0*DLOG(AZRR) + 0.226D0 ARG = ARG + ARG IF (ARG.LT.(-ELIM)) GO TO 140 80 CONTINUE CALL ZAIRY(ZR, ZI, 0, KODE, YR(1), YI(1), NZ1, IERR) CALL ZAIRY(ZR, ZI, 1, KODE, YR(2), YI(2), NZ2, IERR) IF (ICASE.EQ.2) GO TO 100 CALL ZBIRY(ZR, ZI, 0, KODE, WR(1), WI(1), IERR) CALL ZBIRY(ZR, ZI, 1, KODE, WR(2), WI(2), IERR) IF (KODE.EQ.1) GO TO 90 CVR = AZRR - ZRR CVI = -ZRI CALL ZEXP(CVR, CVI, CVR, CVI) STR = WR(1)*CVR - WI(1)*CVI WI(1) = WR(1)*CVI + WI(1)*CVR WR(1) = STR STR = WR(2)*CVR - WI(2)*CVI WI(2) = WR(2)*CVI + WI(2)*CVR WR(2) = STR 90 CONTINUE CVR = CON3R CVI = CON3I NZ3 = 0 NZ4 = 0 GO TO 120 100 CONTINUE CVR = ZR*CON1R - ZI*CON1I CVI = ZR*CON1I + ZI*CON1R CALL ZAIRY(CVR, CVI, 0, KODE, WR(1), WI(1), NZ3, IERR) CALL ZAIRY(CVR, CVI, 1, KODE, WR(2), WI(2), NZ4, IERR) IF (KODE.EQ.1) GO TO 110 C---------------------------------------------------------------------- C WHEN KODE=2,THE SCALING FACTOR EXP(-ZETA1-ZETA2) IS 1.0 FOR C -PI.LT.ARG(Z).LE.PI/3 AND EXP(-2.0*ZETA1) FOR -PI/3.LT.ARG(Z) C .LE.PI WHERE ZETA1=ZETA2 IN THIS RANGE. THIS IS DUE TO THE FACT C THAT ARG(Z*CON1) IS TAKEN TO BE IN (-PI,PI) BY THE PRINCIPAL C SQUARE ROOT. C---------------------------------------------------------------------- IF (IT.LT.7) GO TO 110 CVR = ZRR + ZRR CVI = ZRI + ZRI CALL ZEXP(-CVR, -CVI, CVR, CVI) STR = WR(1)*CVR - WI(1)*CVI WI(1) = WR(1)*CVI + WI(1)*CVR WR(1) = STR STR = WR(2)*CVR - WI(2)*CVI WI(2) = WR(2)*CVI + WI(2)*CVR WR(2) = STR 110 CONTINUE STR = WR(2)*CON1R - WI(2)*CON1I WI(2) = WR(2)*CON1I + WI(2)*CON1R WR(2) = STR CVR = CON2R CVI = CON2I 120 CONTINUE AV = ZABS(CVR,CVI) C ERROR RELATIVE TO MAXIMUM TERM 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) CYR = CWR - CYR - CVR CYI = CWI - CYI - CVI ACY = ZABS(YR(1),YI(1))*ZABS(WR(2),WI(2)) ACW = ZABS(WR(1),WI(1))*ZABS(YR(2),YI(2)) AV = DMAX1(ACW,ACY,AV) ER = ZABS(CYR,CYI)/AV IF (ER.LT.ERTOL) GO TO 140 IF (LFLG.EQ.1) GO TO 130 WRITE (LUN,99995) ERTOL 99995 FORMAT (/43H CASES WHICH VIOLATE THE RELATIVE ERROR TES, * 14HT WITH ERTOL =, D12.4/) WRITE (LUN,99994) 99994 FORMAT (/14H OUTPUT FORMAT/25H KODE,IR,IT,NZ1,NZ2,NZ3,N, * 14HZ4,IRSET,ICASE) WRITE (LUN,99993) 99993 FORMAT (3H ER/22H Z,Y(1),Y(2),W(1),W(2)) LFLG = 1 130 CONTINUE WRITE (LUN,99992) KODE, IR, IT, NZ1, NZ2, NZ3, NZ4, * IRSET, ICASE 99992 FORMAT (9I5) WRITE (LUN,99991) ER WRITE (LUN,99991) ZR, ZI, YR(1), YI(1), YR(2), YI(2), * WR(1), WI(1), WR(2), WI(2) 99991 FORMAT (10D12.4) 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE IF (LFLG.EQ.0) WRITE (LUN,99990) 99990 FORMAT (/16H QUICK CHECKS OK/) STOP END .