code 35140; procedure AIRY(Z,AI,AID,BI,BID,EXPON,FIRST); value Z,FIRST; boolean FIRST; real Z,AI,AID,BI,BID,EXPON; begin real S,T,U,V,SC,TC,UC,VC,X,K1,K2,K3,K4, C,ZT,SI,CO,EXPZT,SQRTZ,WWL,PL,PL1,PL2,PL3; own real C1,C2,SQRT3,SQRT1OPI,PIO4; own real array XX,WW[1:10]; integer N,L; if FIRST then begin SQRT3:= 1.73205080756887729; SQRT1OPI:= 0.56418958354775629; PIO4:= 0.78539816339744831; C1:= 0.35502 80538 87817; C2:= 0.25881 94037 92807; XX[ 1]:= 1.40830 81072 180964 "+1; XX[ 2]:= 1.02148 85479 197331 "+1; XX[ 3]:= 7.44160 18450 450930 ; XX[ 4]:= 5.30709 43061 781927 ; XX[ 5]:= 3.63401 35029 132462 ; XX[ 6]:= 2.33106 52303 052450 ; XX[ 7]:= 1.34479 70824 609268 ; XX[ 8]:= 6.41888 58369 567296 "-1; XX[ 9]:= 2.01003 45998 121046 "-1; XX[10]:= 8.05943 59172 052833 "-3; WW[ 1]:= 3.15425 15762 964787"-14; WW[ 2]:= 6.63942 10819 584921"-11; WW[ 3]:= 1.75838 89061 345669"- 8; WW[ 4]:= 1.37123 92370 435815"- 6; WW[ 5]:= 4.43509 66639 284350"- 5; WW[ 6]:= 7.15550 10917 718255"- 4; WW[ 7]:= 6.48895 66103 335381"- 3; WW[ 8]:= 3.64404 15875 773282"- 2; WW[ 9]:= 1.43997 92418 590999"- 1; WW[10]:= 8.12311 41336 261486"- 1; end; EXPON:= 0; if Z >= -5.0 and Z <= 8 then begin U:= V:= T:= UC:= VC:= TC:= 1; S:= SC:= 0.5; N:= 0; X:= Z*Z*Z; for N:= N+3 while ABS(U)+ABS(V)+ABS(S)+ABS(T) > "-18 do begin U:=U*X/(N*(N-1)); V:= V*X/(N*(N+1)); S:=S*X/(N*(N+2)); T:= T*X/(N*(N-2)); UC:= UC+U; VC:= VC+V; SC:= SC+S; TC:= TC+T end; BI:= SQRT3 * (C1*UC + C2*Z*VC); BID:=SQRT3 * (C1*Z*Z*SC +C2*TC); if Z<2.5 then begin AI:= C1*UC - C2*Z*VC; AID:= C1*SC*Z*Z - C2*TC; goto END end end; K1:= K2:= K3:= K4:= 0; SQRTZ:= SQRT(ABS(Z)); ZT:= 0.66666 66666 66667 * ABS(Z)*SQRTZ; C:= SQRT1OPI/SQRT(SQRTZ); if Z<0 then begin Z:= -Z; CO:= COS(ZT-PIO4); SI:= SIN(ZT-PIO4); for L:= 1 step 1 until 10 do begin WWL:= WW[L]; PL:= XX[L]/ZT; PL2:=PL*PL; PL1:= 1+PL2; PL3:= PL1*PL1; K1:= K1 + WWL/PL1; K2:= K2 + WWL*PL/PL1; K3:= K3 + WWL*PL*(1+PL*(2/ZT+PL))/PL3; K4:= K4 + WWL*(-1-PL*(1+PL*(ZT-PL))/ZT)/PL3; end; AI:= C*(CO*K1+SI*K2); AID:= 0.25*AI/Z - C*SQRTZ*(CO*K3+SI*K4); BI:= C*(CO*K2-SI*K1); BID:= 0.25*BI/Z - C*SQRTZ*(CO*K4-SI*K3); end else begin if Z < 9 then EXPZT:= EXP(ZT) else begin EXPZT:= 1; EXPON:= ZT end; for L:= 1 step 1 until 10 do begin WWL:= WW[L]; PL:= XX[L]/ZT; PL1:= 1+PL; PL2:= 1-PL; K1:= K1 + WWL/PL1; K2:= K2 + WWL*PL/(ZT*PL1*PL1); K3:= K3 + WWL/PL2; K4:= K4 + WWL*PL/(ZT*PL2*PL2); end; AI:= 0.5*C*K1/EXPZT; AID:= AI*(-.25/Z-SQRTZ) + 0.5*C*SQRTZ*K2/EXPZT; if Z >= 8 then begin BI:= C*K3*EXPZT; BID:= BI*(SQRTZ-0.25/Z) - C*K4*SQRTZ*EXPZT; end; end; END: end AIRY eop