%--------------------------------------------------------------------------------------
% Error curves for the Algebraic Ellipse Perimeter Approximations
% and the respective approximations of the elliptic function E(x).
%--------------------------------------------------------------------------------------
%--------------------------------------------------------------------------------------
% Globally Best Algebraic Approximations of Ellipse Perimeters
%--------------------------------------------------------------------------------------
% This program was written in may 2006 by Stan Sykora while writing
% the article on 'Algebraic Approximations of Ellipse Perimeters',
% published online at www.ebyte.it/library/docs/math05a/EllipseAlgApprox06.html.
% The program is an attachment to the said article and was used to produce Fig.1.
%--------------------------------------------------------------------------------------
% The AUTHOR hereby authorizes EVERYBODY to copy, distribute, cut, or modify
% the program in whichever form and by whatever means. If you use it and obtain
% interesting results, however, you should cite the Web article.
%
% The AUTOR DISCLAIMS any LIABILITIES which any misguided soul might claim,
% due to whatever he/she has done with the program or just thinks to have done.

clear;clc;clf;echo off;
logscale = 1;                   % set 0 for linear scale

y = logspace(-4,0,10000);
y1 = 1+y;

[K,E] = ellipke(1-y.*y,1e-15);   %Compute the elliptic function E to the highest precision
                                 %Actually just 8 digits can be trusted

  
%Cantrell Particularly Fruitful = Rational1, 6313 ppm
t1 = (pi-2)/2; % 0.57079632679490
Rational1 = ((1+y.^2)+2*t1*y)./(1+y);
plot(y,(Rational1-E)./E,'-k','LineWidth',2,'Color',[0.5 0.3 0]);
hold on;
err1 = max(abs((Rational1-E)./E))

%Rational2, 557.3 ppm
p = [1.22694921875000];
s1 = p(1);
t1 = (pi/2)*(1+s1)-1; % 2.49808365277126
Rational2 = ((1+y.^3)+t1*(1+y).*y)./((1+y.^2)+2*s1*y);
plot(y,(Rational2-E)./E,'-k','LineWidth',2,'Color',[1 0 0]);
hold on;
plot(y,-(Rational2-E)./E,'-k','LineWidth',1,'Color',[1 0 0]);
hold on;
err2 = max(abs((Rational2-E)./E))

%Sykora 2005 = Rational3, 75.96 ppm
p = [6.15241658169936   6.16881239339582];
s1 = p(1);
t1 = p(2);
t2 = (pi/2)*(1+s1)-1-t1; % 4.06617730084445 
Rational3 = ((1+y.^4)+t1*y.*(1+y.^2)+2*t2*y.^2)./((1+y.^3)+s1*y.*(1+y));
plot(y,(Rational3-E)./E,'-k','LineWidth',2,'Color',[1 0.7 0]);
hold on;
plot(y,-(Rational3-E)./E,'-k','LineWidth',1,'Color',[1 0.7 0]);
hold on;
err3 = max(abs((Rational3-E)./E))

%Rational4, 15.29 ppm
p = [13.01750519704827  13.09922140579137  13.02487942169925];
s1 = p(1);
s2 = p(2);
t1 = p(3);
t2 = (pi/2)*(1+s1+s2)-1-t1; % 28.56997512074272
Rational4 = ((1+y.^5)+t1*y.*(1+y.^3)+t2*(y.^2).*(1+y))./((1+y.^4)+s1*y.*(1+y.^2)+2*s2*y.^2);
plot(y,(Rational4-E)./E,'-k','LineWidth',2,'Color',[0 0.6 0]);
hold on;
plot(y,-(Rational4-E)./E,'-k','LineWidth',1,'Color',[0 0.6 0]);
hold on;
err4 = max(abs((Rational4-E)./E))

%Rational5 = Cantrell 2006, 3.65 ppm
p = 100*[0.27297854333670   1.10551872985869   0.27301243680755   1.13302483206429];
s1 = p(1);
s2 = p(2);
t1 = p(3);
t2 = p(4);
t3 = (pi/2)*(1+s1+s2)-1-t2-t1; % 76.50091476282086
Rational5 = ((1+y.^6)+t1*y.*(1+y.^4)+t2*(y.^2).*(1+y.^2)+2*t3*(y.^3))./((1+y.^5)+s1*y.*(1+y.^3)+s2*(y.^2).*(1+y));
plot(y,(Rational5-E)./E,'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;
plot(y,-(Rational5-E)./E,'-k','LineWidth',1,'Color',[0 0 0.8]);
hold on;
err5 = max(abs((Rational5-E)./E))


%Rational6, 988 ppb
p = 100*[0.51445996674310   3.82256974433855   3.47212007887276   0.51447782789130   3.85327854851892];
s1 = p(1);
s2 = p(2);
s3 = p(3);
t1 = p(4);
t2 = p(5);
t3 = (pi/2)*(1+s1+s2+s3)-1-t2-t1; % 7.904535392309255e+002
Rational6 = ((1+y.^7)+t1*y.*(1+y.^5)+t2*(y.^2).*(1+y.^3)+t3*(y.^3).*(1+y))./((1+y.^6)+s1*y.*(1+y.^4)+s2*(y.^2).*(1+y.^2)+2*s3*(y.^3));
plot(y,(Rational6-E)./E,'-k','LineWidth',2,'Color',[0 0 0]);
hold on;
plot(y,-(Rational6-E)./E,'-k','LineWidth',1,'Color',[0 0 0]);
hold on;
err7 = max(abs((Rational6-E)./E))

%Rational7, 282 ppb
p = 1000*[0.09349135794687   1.25936473022183   4.09396201082922   0.09349235523473   1.26273239571330   4.29645229646421];
s1 = p(1);
s2 = p(2);
s3 = p(3);
t1 = p(4);
t2 = p(5);
t3 = p(6);
t4 = (pi/2)*(1+s1+s2+s3)-1-t2-t3-t1; % 2.903735611540449e+003
Rational7 = ((1+y.^8)+t1*y.*(1+y.^6)+t2*(y.^2).*(1+y.^4)+t3*(y.^3).*(1+y.^2)+2*t4*(y.^4))./((1+y.^7)+s1*y.*(1+y.^5)+s2*(y.^2).*(1+y.^3)+s3*(y.^3).*(1+y));
plot(y,(Rational7-E)./E,'-k','LineWidth',2,'Color',[1 0 0]);
hold on;
plot(y,-(Rational7-E)./E,'-k','LineWidth',1,'Color',[1 0 0]);
hold on;
err7 = max(abs((Rational7-E)./E))

%--------------------------------------------------------------------------------------

set(gcf,'Color','white');
set(gca,'FontSize',16);
set(gca,'LineWidth',2);
set(gca,'TickLength',[.03 .03]);
set(gca,'yScale','log');
axis([0.0001 1 1e-7 1e-2]); %Below 1e-8 the precision of the program may be insufficient
if logscale
  set(gca,'xScale','log'); %Comment this line for linear x-scale
else
  set(gca,'xTick',[.2,.4,.6,.8,1]);
end