'********************************************************************** '* Purpose: This program computes Airy functions and their * '* derivatives using subroutine AIRYA * '* Input: xstart, xend --- Arguments of Airy function * '* xstep --- x increment * '* Output: AI --- Ai(x) * '* BI --- Bi(x) * '* AD --- Ai'(x) * '* BD --- Bi'(x) * '* Example: * '* xstart =0 xend = 30 xstep = 10 * '* * '* x Ai(x) Bi(x) Ai'(x) Bi'(x) * '* ---------------------------------------------------------------- * '* 0 .35502805D+00 .61492663D+00 -.25881940D+00 .44828836D+00 * '* 10 .11047533D-09 .45564115D+09 -.35206337D-09 .14292361D+10 * '* 20 .16916729D-26 .21037650D+26 -.75863916D-26 .93818393D+26 * '* 30 .32082176D-48 .90572885D+47 -.17598766D-47 .49533045D+48 * '* * '* x Ai(-x) Bi(-x) Ai'(-x) Bi'(-x) * '* ---------------------------------------------------------------- * '* 0 .35502805 .61492663 -.25881940 .44828836 * '* 10 .04024124 -.31467983 .99626504 .11941411 * '* 20 -.17640613 -.20013931 .89286286 -.79142903 * '* 30 -.08796819 -.22444694 1.22862060 -.48369473 * '* ------------------------------------------------------------------ * '* REFERENCE: "Fortran Routines for Computation of Special Functions, * '* jin.ece.uiuc.edu/routines/routines.html". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************************** ' PROGRAM MAIRYA DEFINT I-N DEFDBL A-H, O-Z INPUT " Please enter x start, x end, x step: ", XSTART, XEND, XSTEP X = XSTART 5 CALL AIRYA(X, AI, BI, AD, BD) IF (X = XSTART) THEN PRINT " x Ai(x) Bi(x) Ai'(x) Bi'(x) " PRINT " -------------------------------------------------------------------------" END IF PRINT X; TAB(5); AI; TAB(30); BI; TAB(55); AD PRINT TAB(5); BD X = X + XSTEP IF (X <= XEND) THEN GOTO 5 X = XSTART PRINT " " 7 CALL AIRYA(-X, AI, BI, AD, BD) IF (X = XSTART) THEN PRINT " x Ai(-x) Bi(-x) Ai'(-x) Bi'(-x) " PRINT " -------------------------------------------------------------------------" END IF PRINT X; TAB(5); AI; TAB(30); BI; TAB(55); AD PRINT TAB(5); BD X = X + XSTEP IF (X <= XEND) THEN GOTO 7 PRINT END SUB AIRYA (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) ' Routine called: ' AJYIK for computing Jv(x), Yv(x), Iv(x) and ' Kv(x) with v=1/3 and 2/3 ' ====================================================== XA = ABS(X) PIR = 0.318309886183891 C1 = 0.355028053887817 C2 = 0.258819403792807 SR3 = 1.732050807568877 Z = XA ^ 1.5 / 1.5 XQ = SQR(XA) CALL AJYIK(Z, VJ1, VJ2, VY1, VY2, VI1, VI2, VK1, VK2) IF X = 0.0 THEN AI = C1 BI = SR3 * C1 AD = -C2 BD = SR3 * C2 ELSEIF X > 0.0 THEN AI = PIR * XQ / SR3 * VK1 BI = XQ * (PIR * VK1 + 2.0 / SR3 * VI1) AD = -XA / SR3 * PIR * VK2 BD = XA * (PIR * VK2 + 2.0 / SR3 * VI2) ELSEIF X < 0.0 THEN AI = 0.5 * XQ * (VJ1 - VY1 / SR3) BI = -0.5 * XQ * (VJ1 / SR3 + VY1) AD = 0.5 * XA * (VJ2 + VY2 / SR3) BD = 0.5 * XA * (VJ2 / SR3 - VY2) END IF END SUB SUB AJYIK (X, VJ1, VJ2, VY1, VY2, VI1, VI2, VK1, VK2) ' ======================================================= ' Purpose: Compute Bessel functions Jv(x) and Yv(x), ' and modified Bessel functions Iv(x) and ' Kv(x), and their derivatives with v=1/3,2/3 ' Input : x --- Argument of Jv(x),Yv(x),Iv(x) and ' Kv(x) ( x ò 0 ) ' Output: VJ1 --- J1/3(x) ' VJ2 --- J2/3(x) ' VY1 --- Y1/3(x) ' VY2 --- Y2/3(x) ' VI1 --- I1/3(x) ' VI2 --- I2/3(x) ' VK1 --- K1/3(x) ' VK2 --- K2/3(x) ' ======================================================= IF (X = 0.0) THEN VJ1 = 0.0 VJ2 = 0.0 VY1 = -1.0E+100 VY2 = 1.0E+100 VI1 = 0.0 VI2 = 0.0 VK1 = -1.0E+100 VK2 = -1.0E+100 GOTO 100 END IF PI = 3.141592653589793 RP2 = .63661977236758 GP1 = .892979511569249 GP2 = .902745292950934 GN1 = 1.3541179394264 GN2 = 2.678938534707747 VV0 = 0.444444444444444 UU0 = 1.1547005383793 X2 = X * X K0 = 12 IF (X >= 35.0) THEN K0 = 10 IF (X >= 50.0) THEN K0 = 8 IF (X <= 12.0) THEN FOR L = 1 TO 2 VL = L / 3.0 VJL = 1.0 R = 1.0 FOR K = 1 TO 40 R = -0.25 * R * X2 / (K * (K + VL)) VJL = VJL + R IF (ABS(R) < 1.0E-15) THEN GOTO 20 NEXT K 20 A0 = (0.5 * X) ^ VL IF (L = 1) THEN VJ1 = A0 / GP1 * VJL IF (L = 2) THEN VJ2 = A0 / GP2 * VJL NEXT L ELSE FOR L = 1 TO 2 VV = VV0 * L * L PX = 1.0 RP = 1.0 FOR K = 1 TO K0 RP = -0.78125D-2 * RP * (VV - (4.0 * K - 3.0) ^ 2.0) * (VV - (4.0 * K - 1.0) ^ 2.0) / (K * (2.0 * K - 1.0) * X2) PX = PX + RP NEXT K QX = 1.0 RQ = 1.0 FOR K = 1 TO K0 RQ = -0.78125E-2 * RQ * (VV - (4.0 * K - 1.0) ^ 2.0) * (VV - (4.0 * K + 1.0) ^ 2.0) / (K * (2.0 * K + 1.0) * X2) QX = QX + RQ NEXT K QX = 0.125 * (VV - 1.0) * QX / X XK = X - (0.5 * L / 3.0 + 0.25) * PI A0 = SQR(RP2 / X) CK = COS(XK) SK = SIN(XK) IF (L = 1) THEN VJ1 = A0 * (PX * CK - QX * SK) VY1 = A0 * (PX * SK + QX * CK) ELSEIF (L = 2) THEN VJ2 = A0 * (PX * CK - QX * SK) VY2 = A0 * (PX * SK + QX * CK) END IF NEXT L END IF IF (X <= 12.0) THEN FOR L = 1 TO 2 VL = L / 3.0 VJL = 1.0 R = 1.0 FOR K = 1 TO 40 R = -0.25 * R * X2 / (K * (K - VL)) VJL = VJL + R IF (ABS(R) < 1.0E-15) THEN GOTO 50 NEXT K 50 B0 = (2.0 / X) ^ VL IF (L = 1) THEN UJ1 = B0 * VJL / GN1 IF (L = 2) THEN UJ2 = B0 * VJL / GN2 NEXT L PV1 = PI / 3.0 PV2 = PI / 1.5 VY1 = UU0 * (VJ1 * COS(PV1) - UJ1) VY2 = UU0 * (VJ2 * COS(PV2) - UJ2) END IF IF (X <= 18.0) THEN FOR L = 1 TO 2 VL = L / 3.0 VIL = 1.0 R = 1.0 FOR K = 1 TO 40 R = 0.25 * R * X2 / (K * (K + VL)) VIL = VIL + R IF (ABS(R) < 1.0E-15) THEN GOTO 65 NEXT K 65 A0 = (0.5 * X) ^ VL IF (L = 1) THEN VI1 = A0 / GP1 * VIL IF (L = 2) THEN VI2 = A0 / GP2 * VIL NEXT L ELSE C0 = EXP(X) / SQR(2.0 * PI * X) FOR L = 1 TO 2 VV = VV0 * L * L VSL = 1.0 R = 1.0 FOR K = 1 TO K0 R = -0.125 * R * (VV - (2.0 * K - 1.0) ^ 2.0) / (K * X) VSL = VSL + R NEXT K IF (L = 1) THEN VI1 = C0 * VSL IF (L = 2) THEN VI2 = C0 * VSL NEXT L END IF IF (X <= 9.0) THEN FOR L = 1 TO 2 VL = L / 3.0 IF (L = 1) THEN GN = GN1 IF (L = 2) THEN GN = GN2 A0 = (2.0 / X) ^ VL / GN SUM = 1.0 R = 1.0 FOR K = 1 TO 60 R = 0.25 * R * X2 / (K * (K - VL)) SUM = SUM + R IF (ABS(R) < 1.0E-15) THEN GOTO 90 NEXT K 90 IF (L = 1) THEN VK1 = 0.5 * UU0 * PI * (SUM * A0 - VI1) IF (L = 2) THEN VK2 = 0.5 * UU0 * PI * (SUM * A0 - VI2) NEXT L ELSE C0 = EXP(-X) * SQR(0.5 * PI / X) FOR L = 1 TO 2 VV = VV0 * L * L SUM = 1.0 R = 1.0 FOR K = 1 TO K0 R = 0.125 * R * (VV - (2.0 * K - 1.0) ^ 2.0) / (K * X) SUM = SUM + R NEXT K IF (L = 1) THEN VK1 = C0 * SUM IF (L = 2) THEN VK2 = C0 * SUM NEXT L END IF 100 END SUB