%Numerical calulation of persistence probability of Airy1 process based on
%Bornemann's method. 


%The following function is numerical calulation of (2.3) 

function [result] = persistenceProbability(L,m,c)
    [wx x]=ClenshawCurtis(-40,40,m);

    IdentityMatrix = eye(m);
    persistenceMatrix = zeros(m);
    
    for i=1:m
        for j=1:m
            persistenceMatrix(i,j)=wx(i)^(1/2)*wx(j)^(1/2)*exp(L*(x(j)-x(i)))*persistenceKernel1(L,x(i),x(j)+2*c);
        end
    end

    result = det(IdentityMatrix-persistenceMatrix);
end

%Kernel defined in Proposition 2.2
function [A] = persistenceKernel1(L,x,y)

    temp1=P_0_B_0(x,y);
    temp2=tildeK(x,y,L);
    temp3=hatK(x,y,L);
    A=temp1+temp2+temp3;
end

% This is the function 1_{x>=0}Airy(x,y), where Airy(*) is the usual Airy function. 

function [sum] = P_0_B_0(x,y)
sum=0;
if x>0
    sum = airy(x+y);
end
end
%Kernenl \tilde K_{L,c} defined in (2.4) 
function [sum] = tildeK(x,y,L)
sum = 0;
if x<0
[w,z]=ClenshawCurtis(0,inf,512);

for i =1:length(z)
    sum = sum+ exp(-(x-z(i))^2/(4*L)-2*L^3/3-L*(y+z(i)))/sqrt(4*pi*L)*airy(z(i)+y+L^2)*w(i);
end
end

end
%Kernel \hat K_{L,c} defined in (2.4)

function [sum] = hatK(x,y,L)
sum = 0;
if x<0
[w,z]=ClenshawCurtis(-inf,0,512);

for i =1:length(z)
    sum = sum+ exp(-(x+z(i))^2/(4*L)-2*L^3/3-L*(y+z(i)))/sqrt(4*pi*L)*airy(z(i)+y+L^2)*w(i);
end
end

end