##  A Maple package to compute Airy-type assymptotic expansions,
##     supplementing the paper
##     "Symbolic Evaluation of Coefficients
##      in Airy-type Asymptotic Expansions"
##     by Raimundas Vidunas (University of Amsterdam)
##     and Nico Temme (CWI, Amsterdam).
##  The considered integrals are of form (1) in this paper,
##     with f(t) given by the power series (9)
##     at the two saddle points (b and -b).
##  The position of saddle points is given by the variable AiryB.
##  To give the input coefficients in (9) the user should (re)define
##     functions AiryPw(k) and AiryPm(k).
##  The output coefficients are returned by functions
##     AiryAlpha(k) and AiryBeta(k).
##  The intermediate coefficients are given by functions
##     AiryFo, AiryFe, AiryFoe, AiryC, AiryD, AiryGamma, AiryDelta.
##   
##  (c) Raimundas Vidunas
##  Date: 17 Jan 2001
##  Version: 3.0
##

`This package computes Airy-type asymptotic expansion of the integral`;
1/(2*'pi'*'i'),int(exp('z*(w^3/3-eta*w)')*'f(w)', 'w'), `.`;
`For the input one has to redefine Maple functions`;
'AiryPw(k)', ` and `, 'AiryPm(k)';
`which specify the k-th coefficient of `, 'f(w)',
 ` and `, 'f(-w)', ` at the saddle point`, 'w=AiryB', `.`;
`Here `, 'AiryB=sqrt(eta)', `, is a global variable`;
`For example, for general `, 'f(w)',`one may define:`;
`AiryPw:= k->p[k];     AiryPm:= k->q[k];`;
if (eval(AiryPw)='AiryPw') or (eval(AiryPm)='AiryPm') then
  `And for `, 'f(w)=1/(1-w)', ` we have:`;
  'AiryPw(k)=(-1)^k/(1+AiryB)^k', ` `, 'AiryPm(k)=1/(1-AiryB)^k'
fi;
`For the output, call functions `, 'AiryAlpha(k) and AiryBeta(k)';
`to get the coefficients in the assymptotic expansion`;
'z^`-1/3`*Ai(eta*z^`2/3`)*sum((-1)^k*alpha[k]/z^k,k=0..infinity)
-z^`-2/3`*D(Ai)(eta*z^`2/3`)*sum((-1)^k*beta[k]/z^k,k=0..infinity)';
`These coeffcients depend on `*'eta', `(or equivalently, on AiryB)`;
`It is convenient to rename AiryB, for example `, 'alias(b=`AiryB `)';
`(However, either define AipyPw, AiryPm using AiryB, or first define alias)`;
if (AiryPw='ParCyPw') or (AiryPm='ParCyPm') then
  print(`You already have chosen to work with Weber parabolic cylinder function`);
  print('AiryPw'=AiryPw);   print('AiryPm'=AiryPm);
  print(`For more information please call `*'ParCyInfo()')
else
  `To work with Weber parabolic cylinder function, call`*'ParCyInfo()';
  if (eval(AiryPw)='AiryPw') or (eval(AiryPm)='AiryPm') then
    print(`Please define functions`, 'AiryPw(k) and AiryPm(k)');
  else 
    print(`At the moment you have defined:`);
    print('AiryPw'=eval(AiryPw));
    print('AiryPm'=eval(AiryPm));
  fi;
fi;


AiryAlpha:= proc(n)  normal(AiryGamma(n,0))  end:
AiryBeta:=  proc(n)  normal(AiryDelta(n,0))  end:

AiryGamma:= proc( n, k )
  if n=0 then  AiryC(k)
  else factor(
   (2*k+1)*AiryDelta(n-1,k+1)+2*AiryB^2*(k+1)*AiryDelta(n-1,k+2) )
  fi
end:

AiryDelta:= proc(n, k)
  if n=0 then  AiryD(k)
  else 2*(k+1)*AiryGamma(n-1,k+2)
  fi
end:

AiryC:= proc( k )
local j;
  if k=0 then  AiryFe(0)
  else factor( sum( '(-1)^(k-j)*j/(2*k-j)*
    binomial(2*k-j,k)/(2*AiryB)^(2*k-j)*AiryFe(j)', 'j'=1..k ) )
  fi
end:

AiryD:= proc( k )
local j;
  if k=0 then  AiryFoe(0)
  else factor( sum( '(-1)^(k-j)*j/(2*k-j)*
    binomial(2*k-j,k)/(2*AiryB)^(2*k-j)*AiryFoe(j)', 'j'=1..k ) )
  fi
end:

AiryFoe:= proc( k )
  if k<0 then  0
    else expand( (AiryFo(k)-AiryFoe(k-1))/AiryB );
  fi
end:

AiryFe:= proc(k)  (AiryPw(k)+AiryPm(k))/2  end:
AiryFo:= proc(k)  (AiryPw(k)-AiryPm(k))/2  end:

###
### Weber Parabolic Cylinder function
###
# The variable ParCySU is the power series expansion of u in b.
# The variable ParCySXi is the power series expansion of xi in b.
# Other used names: ParCySw, ParCySm, ParCySqrtS.
#

ParCyInfo:= proc()
if (AiryPw='ParCyPw') or (AiryPm='ParCyPm') then
  print(`You already have chosen to work with Weber parabolic cylinder function`);
  print('AiryPw'=AiryPw);
  print('AiryPm'=AiryPm)
else
  print(`To work with Weber parabolic cylinder function, you have to assign:`);
  print(`AiryPw:=ParCyPw;   AiryPm:=ParCyPm;`)
fi;    
print(`Parabolic cylinder functions are solutions of differential equation`);
print('diff(y(x),x$2)=(x^2/4+a)*y(x)');
print(`Following paper [VT], we take the exponentially decreasing solution`);
print('U(a,x)',` = `,'exp(x^2/4)/i/sqrt(2*pi)',
  'int(exp(-x*s+s^2/2)*s^(-a-1/2), s)'); 
print(`This integral representation can be analytically transformed to`);
print(`a standard Airy-type integral. The obtained Airy-type asymptotic expansion`);
print(`is in the (non-positive) powers of`, 'z=sqrt(-a)');
print(`The coefficients depend on mutally related parameters`);
print('2/3*AiryB^3=t*sqrt(t^2-1)-log(t+sqrt(t^2-1))', ` where `, 't=x/2/sqrt(-a)');
#print('ParCyU=sqrt(2*AiryB)*((x-2*sqrt(-a))/(x+2*sqrt(-a)))^`1/4`');
print('ParCyU=sqrt(2*AiryB)*((t-1)/(t+1))^`1/4`', 'ParCyXi=t*AiryB/sqrt(t^2-1)');
#'ParCyXi=(ParCyU^4+4*ParCyB^2)/(4*ParCyB^2)')
print(`Coefficients in `, 'AiryB and ParCyU', ` are returned by `,
  'AiryAlpha(k), AiryBeta(k)');
print(`Coefficients in `, 'AiryB and ParCyXi', ` are returned by `,
  'ParCyAlpha(k), ParCyBeta(k)');
print(`It is convenient to rename AiryB, ParCyU and ParCyXi, say`);
print('alias(b=`AiryB `, u=`ParCyU `, xi=`ParCyXi `)')
end:

alias( ParCyUa=RootOf( z^4-4*ParCyXi*z^2+4*AiryB^2, z)):
# Algebraic relation between ParCyU and ParCyXi

ParCySXi:=1+2/5*AiryB^2-9/175*AiryB^4+92/7875*AiryB^6-1894/606375*AiryB^8
+19758/21896875*AiryB^10-2418092/8868234375*AiryB^12
+89543912/1055319890625*AiryB^14-1834815291/68074647265625*AiryB^16
+9059400557356/1042154774989453125*AiryB^18
-401773649760188/141638308055384765625*AiryB^20
+116469314266892/124670207136083984375*AiryB^22
-643655907483649612/2076363119819611669921875*AiryB^24+O(AiryB^26):

ParCySU:=AiryB-3/40*AiryB^3+69/4480*AiryB^5-1261/322560*AiryB^7
+7221551/6623232000*AiryB^9-1841161947/5740134400000*AiryB^11
+25424826618691/260372496384000000*AiryB^13
-51337733788603/1686221881344000000*AiryB^15
+30418747093245561/3148690285920256000000*AiryB^17+O(AiryB^18):

ParCyAlpha:=  proc(k)  factor( evala(subs(ParCyU=ParCyUa,AiryAlpha(k))) ) end:
ParCyBeta:= proc(k)  factor( evala(subs(ParCyU=ParCyUa,AiryBeta(k))) ) end:

ParCySw:= proc( k )
option remember;
local T, s, a, w;
if k=0 then  (2*AiryB+ParCyU^2)/(2*AiryB-ParCyU^2)
elif k=1 then  (2*AiryB+ParCyU^2)/2/ParCyU
else
  s:= sum('a[i]*w^i', 'i'=0..k);
  T:= coeff( expand(
   (s^2+2*(ParCyU^4+4*AiryB^2)/(ParCyU^4-4*AiryB^2)*s+1)*diff(s,w)
   -w*(w+2*AiryB)*s ), w, k);
  sort( factor( solve(
    subs( seq(a[i]=ParCySw(i),i=0..k-1), T), a[k])), ParCyU)
fi
end:

ParCySm:= proc( k )
  subs( ParCyU=-ParCyU, AiryB=-AiryB, ParCySw(k) )
end:

ParCySqrtS:= proc( k )
option remember;
local j;
  if k=0 then  4*ParCyU/(2*AiryB-ParCyU^2)
  else
    factor(  sum('(3/2*j-k)*ParCySw(j)*ParCySqrtS(k-j)','j'=1..k)
        /ParCySw(0)/k )
  fi
end:

ParCyPw:= proc(k)  (k+1)*ParCySqrtS(k+1)  end:

ParCyPm:= proc( k )
  (-1)^k*subs(AiryB=-AiryB, ParCyU=-ParCyU, ParCyPw(k))
end: