airy ai and bi functions

[roadrunner]

Functions for calculation of airy ai and bi follow.

Code:


local c(k)
BEGIN
 LOCAL j, n;
 n:=1/((216^k)*k!);
 FOR j FROM 0 TO 2*k-1 DO
  n:=n*(2*k+2*j+1);
 END;
 RETURN n;
END;

local f(z)
begin
 local t,t2,k,p,p2;
 t:=1;
 t2:=0;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k;
  p2:=p2*(3*k-2);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

local g(z)
begin
 local t,t2,k,p,p2;
 t:=z;
 t2:=0;
 p2:=1;
 for k from 1 to 80 do
  p:=3*k+1;
  p2:=p2*(3*k-1);
  t:=t+p2*z^(p)/(p)!;
  if t == t2 then return t; end;
  t2:=t;
 end;
 return t;
end;

local aismallz(z)
begin
 return
 0.355028053887817*f(z)-
 0.258819403792807*g(z);
end;

local aibigz(z)
begin
 local s, k, old1, old2, new1, new2;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  new1:=1; 
  for k from 1 to 48 do
   old1:=new1;
   new1:=new1+(−1)^k*c(k)*s^(−k);
   if new1 == old1 then break; end;
  end; 
  return 0.282094791774/z^(1/4)/exp(s)*new1;
 end; 
 z:=−z;
 s:=2/3*z^(3/2);
 new1:=1;
 new2:=c(1)/s;
 for k from 1 to 24 do
  old1:=new1;
  old2:=new2;
  new1:=new1+(−1)^k*c(2*k)*s^(−2*k);
  new2:=new2+(−1)^k*c(2*k+1)*s^(−2*k-1);
  if (new1 == old1) and (new2 == old2) then
   break;
  end;
 end;
 s:=s+π/4;
 return (0.564189583547/z^(1/4))*
 (sin(s)*new1-cos(s)*new2);
end;

local bismallz(z)
begin
  return (0.355028053887817*f(z)+
  0.258819403792807*g(z))*√3;
end;

local bibigz(z)
begin
 local s, k, old1, old2, new1, new2;
 if z ≥ 0 then
  s:=2/3*z^(3/2);
  new1:=1; 
  for k from 1 to 48 do
   old1:=new1;
   new1:=new1+c(k)/s^k;
   if new1 == old1 then break; end;
  end; 
  return 0.564189583548/z^(1/4)*exp(s)*new1;
 end;
 z:=−z;
 s:=2/3*z^(3/2);
 new1:=1;
 new2:=c(1)/s;
 for k from 1 to 24 do
  old1:=new1;
  old2:=new2;
  new1:=new1+(−1)^k*c(2*k)*s^(−2*k);
  new2:=new2+(−1)^k*c(2*k+1)*s^(−2*k-1);
  if (new1 == old1) and (new2 == old2) then
   break;
  end;
 end;
 s:=s+π/4;
 return (0.564189583548/z^(1/4))*
 (cos(s)*new1+sin(s)*new2);
end;

export Airy_Ai(z)
begin
 HAngle:=0;
 if ABS(z) < 7 then
  return aismallz(z);
 end;
 return aibigz(z);
end;

export Airy_Bi(z)
begin
 HAngle:=0;
 if ABS(z) < 7 then
  return bismallz(z);
 end;
 return bibigz(z);
end;

Road