function [f,fs]=BilliardFunction(x)
% BILLIARDFUNCTION evaluates the billiard function and its derivative
%   [f,fs]=BilliardFunction(x) evaluates the billiard function and
%   its derivative at x using algorithmic differentiation

global px py a
cs=-sin(x);   c=cos(x);
ss=cos(x);    s=sin(x);
xp1s=-cs;     xp1=px-c;
xp2s=-ss;     xp2=py-s;
xq1s=-cs;     xq1=a-c;  
xq2s=-ss;     xq2=-s;
hs=(xp1*xp1s+xp2*xp2s)/sqrt(xp1^2+xp2^2); h=sqrt(xp1^2+xp2^2);
ep1s=(h*xp1s-xp1*hs)/h^2;                 ep1=xp1/h;
ep2s=(h*xp2s-xp2*hs)/h^2;                 ep2=xp2/h;
hs=(xq1*xq1s+xq2*xq2s)/sqrt(xq1^2+xq2^2); h=sqrt(xq1^2 + xq2^2);
eq1s=(h*xq1s-xq1*hs)/h^2;                 eq1=xq1/h;
eq2s=(h*xq2s-xq2*hs)/h^2;                 eq2=xq2/h;
fs=(ep1s+eq1s)*s+(ep1+eq1)*ss-(ep2s+eq2s)*c-(ep2+eq2)*cs;
f=(ep1+eq1)*s-(ep2+eq2)*c;