'******************************************************************* '* Purpose: This program computes the first NT zeros of Airy * '* functions Ai(x) and Ai'(x), and the associated * '* values of Ai(a') and Ai'(a), and the first NT * '* zeros of Airy functions Bi(x) and Bi'(x), and * '* the associated values of Bi(b') and Bi'(b) using * '* subroutine AIRYZO * '* Input : NT --- Total number of zeros * '* KF --- Function code * '* KF=1 for Ai(x) and Ai'(x) * '* KF=2 for Bi(x) and Bi'(x) * '* Output: XA(m) --- a, the m-th zero of Ai(x) or * '* b, the m-th zero of Bi(x) * '* XB(m) --- a', the m-th zero of Ai'(x) or * '* b', the m-th zero of Bi'(x) * '* XC(m) --- Ai(a') or Bi(b') * '* XD(m) --- Ai'(a) or Bi'(b) * '* ( m --- Serial number of zeros ) * '* Example: NT=5 * '* * '* m a Ai'(a) a' Ai(a') * '* ----------------------------------------------------------- * '* 1 -2.33810741 .70121082 -1.01879297 .53565666 * '* 2 -4.08794944 -.80311137 -3.24819758 -.41901548 * '* 3 -5.52055983 .86520403 -4.82009921 .38040647 * '* 4 -6.78670809 -.91085074 -6.16330736 -.35790794 * '* 5 -7.94413359 .94733571 -7.37217726 .34230124 * '* * '* m b Bi'(b) b' Bi(b') * '* ----------------------------------------------------------- * '* 1 -1.17371322 .60195789 -2.29443968 -.45494438 * '* 2 -3.27109330 -.76031014 -4.07315509 .39652284 * '* 3 -4.83073784 .83699101 -5.51239573 -.36796916 * '* 4 -6.16985213 -.88947990 -6.78129445 .34949912 * '* 5 -7.37676208 .92998364 -7.94017869 -.33602624 * '* * '* --------------------------------------------------------------- * '* REFERENCE: "Fortran Routines for Computation of Special Func- * '* tions jin.ece.uiuc.edu/routines/routines.html". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************* DEFDBL A-H, O-Z DEFINT I-N ' PROGRAM MAIRYZO DIM XA(50), XB(50), XC(50), XD(50) PRINT PRINT " KF=1 for Ai(x) and Ai'(x); KF=2 for Bi(x) and Bi'(x)." PRINT " NT is the number of the zeros." INPUT " Please enter KF,NT: ", KF, NT PRINT IF KF = 1 THEN PRINT " m a Ai'(a) a'" PRINT " Ai(a') " ELSEIF KF = 2 THEN PRINT " m b Bi'(b) b'" PRINT " Bi(b') " END IF PRINT "------------------------------------------------------------" CALL AIRYZO(NT, KF, XA(), XB(), XC(), XD()) FOR K = 1 TO NT PRINT K; " "; XA(K); " "; XD(K); " "; XB(K) PRINT " "; XC(K) NEXT K END SUB AIRYZO (NT, KF, XA(), XB(), XC(), XD()) ' ======================================================== ' Purpose: Compute the first NT zeros of Airy functions ' Ai(x) and Ai'(x), a and a', and the associated ' values of Ai(a') and Ai'(a); and the first NT ' zeros of Airy functions Bi(x) and Bi'(x), b and ' b', and the associated values of Bi(b') and ' Bi'(b) ' Input : NT --- Total number of zeros ' KF --- Function code ' KF=1 for Ai(x) and Ai'(x) ' KF=2 for Bi(x) and Bi'(x) ' Output: XA(m) --- a, the m-th zero of Ai(x) or ' b, the m-th zero of Bi(x) ' XB(m) --- a', the m-th zero of Ai'(x) or ' b', the m-th zero of Bi'(x) ' XC(m) --- Ai(a') or Bi(b') ' XD(m) --- Ai'(a) or Bi'(b) ' ( m --- Serial number of zeros ) ' Routine called: AIRYB for computing Airy functions and ' their derivatives ' ======================================================= PI = 3.141592653589793 FOR I = 1 TO NT IF KF = 1 THEN U = 3.0 * PI * (4.0 * I - 1) / 8.0 U1 = 1 / (U * U) RT0 = -(U * U) ^ (1.0 / 3.0) * ((((-15.5902 * U1 + .929844) * U1 - .138889) * U1 + .10416667D0) * U1 + 1.0) ELSEIF KF = 2 THEN IF I = 1 THEN RT0 = -1.17371 ELSE U = 3.0 * PI * (4.0 * I - 3.0) / 8.0 U1 = 1.0 / (U * U) RT0 = -(U * U) ^ (1.0 / 3.0) * ((((-15.5902 * U1 + .929844) * U1 - .138889) * U1 + .10416667) * U1 + 1.0) END IF END IF 10 X = RT0 CALL AIRYB(X, AI, BI, AD, BD) IF (KF = 1) THEN RT = RT0 - AI / AD IF (KF = 2) THEN RT = RT0 - BI / BD IF (ABS((RT - RT0) / RT) > 1.D-9) THEN RT0 = RT GOTO 10 ELSE XA(I) = RT IF (KF = 1) THEN XD(I) = AD IF (KF = 2) THEN XD(I) = BD END IF NEXT I FOR I = 1 TO NT IF (KF = 1) THEN IF (I = 1) THEN RT0 = -1.01879 ELSE U = 3.0 * PI * (4.0 * I - 3.0) / 8.0 U1 = 1 / (U * U) RT0 = -(U * U) ^ (1.0 / 3.0) * ((((15.0168 * U1 - .873954) * U1 + .121528) * U1 - .145833) * U1 + 1.0) END IF ELSEIF (KF = 2) THEN IF (I = 1) THEN RT0 = -2.29444 ELSE U = 3.0 * PI * (4.0 * I - 1.0) / 8.0 U1 = 1.0 / (U * U) RT0 = -(U * U) ^ (1.0 / 3.0) * ((((15.0168 * U1 - .873954) * U1 + .121528) * U1 - .145833) * U1 + 1.0) END IF END IF 20 X = RT0 CALL AIRYB(X, AI, BI, AD, BD) IF (KF = 1) THEN RT = RT0 - AD / (AI * X) IF (KF = 2) THEN RT = RT0 - BD / (BI * X) IF (ABS((RT - RT0) / RT) > 1D-9) THEN RT0 = RT GOTO 20 ELSE XB(I) = RT IF (KF = 1) THEN XC(I) = AI IF (KF = 2) THEN XC(I) = BI END IF NEXT I END SUB SUB AIRYB (X, AI, BI, AD, BD) ' ======================================================= ' Purpose: Compute Airy functions and their derivatives ' Input: x --- Argument of Airy function ' Output: AI --- Ai(x) ' BI --- Bi(x) ' AD --- Ai'(x) ' BD --- Bi'(x) ' ======================================================= DIM CK(40), DK(40) EPS = 1.0D-15 PI = 3.141592653589793 C1 = 0.355028053887817 C2 = 0.258819403792807 SR3 = 1.732050807568877 XA = ABS(X) XQ = SQR(XA) IF (X > 0.0) THEN XM = 5.0 IF (X <= 0.0) THEN XM = 8.0 IF (X = 0.0) THEN AI = C1 BI = SR3 * C1 AD = -C2 BD = SR3 * C2 GOTO 50 'return END IF IF (XA <= XM) THEN FX = 1.0 R = 1.0 FOR K = 1 TO 40 R = R * X / (3.0 * K) * X / (3.0 * K - 1.0) * X FX = FX + R IF (ABS(R / FX) < EPS) THEN GOTO 15 NEXT K 15 GX = X R = X FOR K = 1 TO 40 R = R * X / (3.0 * K) * X / (3.0 * K + 1.0) * X GX = GX + R IF (ABS(R / GX) < EPS) THEN GOTO 25 NEXT K 25 AI = C1 * FX - C2 * GX BI = SR3 * (C1 * FX + C2 * GX) DF = 0.5 * X * X R = DF FOR K = 1 TO 40 R = R * X / (3.0 * K) * X / (3.0 * K + 2.0) * X DF = DF + R IF (ABS(R / DF) < EPS) THEN GOTO 35 NEXT K 35 DG = 1.0 R = 1.0 FOR K = 1 TO 40 R = R * X / (3.0 * K) * X / (3.0 * K - 2.0) * X DG = DG + R IF (ABS(R / DG) < EPS) THEN GOTO 45 NEXT K 45 AD = C1 * DF - C2 * DG BD = SR3 * (C1 * DF + C2 * DG) ELSE XE = XA * XQ / 1.5 XR1 = 1.0 / XE XAR = 1.0 / XQ XF = DSQRT(XAR) RP = 0.5641895835477563 R = 1.0 FOR K = 1 TO 40 R = R * (6.0 * K - 1.0) / 216.0 * (6.0 * K - 3.0) / K * (6.0 * K - 5.0) / (2.0 * K - 1.0) CK(K) = R DK(K) = -(6.0 * K + 1.0) / (6.0 * K - 1.0) * CK(K) NEXT K KM = INT(24.5 - XA) IF (XA < 6.0) THEN KM = 14 IF (XA > 15.0) THEN KM = 10 IF (X > 0.0) THEN SAI = 1.0 SAD = 1.0 R = 1.0 FOR K = 1 TO KM R = -R * XR1 SAI = SAI + CK(K) * R SAD = SAD + DK(K) * R NEXT K SBI = 1.0 SBD = 1.0 R = 1.0 FOR K = 1 TO KM R = R * XR1 SBI = SBI + CK(K) * R SBD = SBD + DK(K) * R NEXT K XP1 = DEXP(-XE) AI = 0.5 * RP * XF * XP1 * SAI BI = RP * XF / XP1 * SBI AD = -0.5 * RP / XF * XP1 * SAD BD = RP / XF / XP1 * SBD ELSE XCS = COS(XE + PI / 4.0) XSS = SIN(XE + PI / 4.0) SSA = 1.0 SDA = 1.0 R = 1.0 XR2 = 1.0 / (XE * XE) FOR K = 1 TO KM R = -R * XR2 SSA = SSA + CK(2 * K) * R SDA = SDA + DK(2 * K) * R NEXT K SSB = CK(1) * XR1 SDB = DK(1) * XR1 R = XR1 FOR K = 1 TO KM R = -R * XR2 SSB = SSB + CK(2 * K + 1) * R SDB = SDB + DK(2 * K + 1) * R NEXT K AI = RP * XF * (XSS * SSA - XCS * SSB) BI = RP * XF * (XCS * SSA + XSS * SSB) AD = -RP / XF * (XCS * SDA + XSS * SDB) BD = RP / XF * (XSS * SDA - XCS * SDB) END IF END IF 50 END SUB 'end of file mairyzo.bas