%--------------------------------------------------------------------------------------
% Error curves for various Ellipse Perimeter Approximations
% and the respective approximations of the elliptic function E(x).
%--------------------------------------------------------------------------------------
% This program was written in November-December 2005 by Stan Sykora while writing
% the article on Ellipse Perimeter Approximations, published on the web URL
% www.ebyte.it\library\docs\math05a\EllipsePerimeterApprox05.html.
% The program is an attachment to the said article.
%--------------------------------------------------------------------------------------
%
% 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, you should cite the Web article. If you don't, nothing
% bad will happen to you immediately but, in a long run, bad conscience will
% play havock with your body chemistry and you will get covered by smelly and
% itching sores.
%
% The AUTOR DISCLAIMS any LIABILITIES which any misguided soul might claim,
% due to whatever he/she has done with the program or, worse still,
% thinks to have done.

clear;clc;clf;echo on;
y = logspace(-4,0,10000);
pi = 3.14159265;
x = 1-y.*y;
h  = ((1-y)./(1+y)).^2;

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

%Enable/Disable the approximation type you prefer (0 value diables the whole class)
%--------------------------------------------------------------------------------------
typeK = 1;
typeO = 1;
typeE = 1;
typeC = 1;


%--------------------------------------------------------------------------------------
% Type K: Keplerian approximations
%--------------------------------------------------------------------------------------
if typeK

%Kepler; quadratic type or zero-th exponent Hoelder mean
Kepler = (pi/2)*sqrt(y);
plot(y,abs((Kepler-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Sipos 1792
Sipos = (pi/2)*(1+y).*(1+y)./((1+sqrt(y)).^2);
plot(y,((Sipos-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Naive
Naive = (pi/4)*(1+y);
plot(y,-((Naive-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Peano 1889
Peano = (pi/8)*(1+y).*(3-sqrt(1-h));
plot(y,abs((Peano-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Euler 1773
Euler = (pi/4)*sqrt(2+2*y.*y);
plot(y,abs((Euler-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Almkvist 1978
Almkvist = (pi/2)*(2*(1+y).^2-(1-sqrt(y)).^4)./((1+sqrt(y)).^2+2*sqrt(2)*sqrt(1+y).*sqrt(sqrt(y)));
plot(y,abs((Almkvist-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Quadratic
Anonymous = (pi/4)*sqrt(2*(1+y.*y)-0.5*(1-y).^2);
plot(y,abs((Anonymous-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Muir 1883; Hoelder mean type
s = 1.5;
Muir = (pi/2)*((1+y.^s)/2).^(1/s);
plot(y,abs((Muir-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%E.S.Selmer II 1975
Selmer2 = (pi/16)*((1+y).*(6+h/2)-sqrt(2+12*y+2*y.*y));
plot(y,abs((Selmer2-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Ramanujan 1997-1920
Ramanujan1 = (pi/4)*(3*(1+y)-sqrt((1+3*y).*(3+y)));
plot(y,abs((Ramanujan1-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

%Ramanujan 1897-1920
Ramanujan2 = (pi/4)*(1+y).*(1+3*h./(10+sqrt(4-3*h)));
plot(y,abs((Ramanujan2-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

end
%--------------------------------------------------------------------------------------

%--------------------------------------------------------------------------------------
if (typeK | typeE)

%Lindner 1904-1920
Lindner = (pi/4)*(1+y).*(1+h/8).^2;

%E.S.Selmer 1975
%It is a Padè11
Selmer = (pi/4)*(1+y).*((16+3*h)./(16-h));

%G.P.Michon 2000+
%It is a Padè12
Michon = (pi/4)*(1+y).*((64+16*h)./(64-h.*h));

%In "A Manual of Mathematics" by Hudson-Lipka 1917, real author unknown
%It is a Padè21
%Hudson = (pi/16)*(1+y).*(3*(1+h/4)+1./(1-h/4));
Hudson = (pi/4)*(1+y).*((64-3*h.*h)./(64-16*h));

%Jacobsen-Waadeland 1985
%It is a Padè22
Jacobsen = (pi/4)*(1+y).*((256-48*h-21*h.*h)./(256-112*h+3*h.*h));

%Padè32
Pade32 = (pi/4)*(1+y).*((3072-1280*h-252*h.*h + 33*h.*h.*h)./(3072-2048*h+212*h.*h));

%Padè33
Pade33 = (pi/4)*(1+y).*((135168-85760*h-5568*h.*h + 3867*h.*h.*h)./(135168-119552*h+22208*h.*h-345*h.*h.*h));

end
%--------------------------------------------------------------------------------------


%--------------------------------------------------------------------------------------
if typeK

plot(y,abs((Lindner-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;
    
plot(y,abs((Selmer-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

plot(y,abs((Michon-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

plot(y,abs((Hudson-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

plot(y,abs((Jacobsen-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

plot(y,abs((Pade32-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

plot(y,abs((Pade33-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0.5]);
hold on;

end
%--------------------------------------------------------------------------------------



%--------------------------------------------------------------------------------------
%Type O: Optimized Keplerian
%--------------------------------------------------------------------------------------
if typeO

%Arithmetic-Geometric
p = 1.32;
AriGeom = (pi/2)*(p*(1+y)/2+(1-p)*sqrt(y));
plot(y,abs((AriGeom-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

%Optimized Ramanujan mean
%Optimum coincides with arithmetic mean

%Optimized quadratic
Quadratic = (pi/4)*sqrt(2*(1+y.*y)-((1-y).^2)/2.458338);
plot(y,abs((Quadratic-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Ramanujan I optimized for p
p = 3.0273;
w = 3;
Ramanujan1OptP = (pi/2)*(p*(1+y)/2+(1-p)*sqrt((1+w*y).*(w+y))/(1+w));
plot(y,abs((Ramanujan1OptP-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Ramanujan I optimized for w
p = 3;
w = 3.02575;
Ramanujan1OptW = (pi/2)*(p*(1+y)/2+(1-p)*sqrt((1+w*y).*(w+y))/(1+w));
plot(y,abs((Ramanujan1OptW-E)./E),'-k','LineWidth',2,'Color',[0.8 0 0]);
hold on;

%Ramanujan I optimized for P,W does not lead anywhere (marginal gains)

end
%--------------------------------------------------------------------------------------



%--------------------------------------------------------------------------------------
%Type E: Exact for y=0 and y=1, with no crossings
%--------------------------------------------------------------------------------------
if typeE
    
   
%Takakazu Seki Kowa (1642-1708); quadratic type
Takakazu = 0.5*sqrt(pi*pi*y+4*(1-y).^2);
plot(y,abs((Takakazu-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

%Bartolomeu-Michon 2004
BartolomeuMichon = (pi/4)*(1-y)./atan((1-y)./(1+y));
plot(y,abs((BartolomeuMichon-E)./E),'-k','LineWidth',2,'Color',[0.8 0.4 0]);
hold on;

%Lockwood 1932
Lockwood = y.*y.*atan(1./y)+(1./y).*atan(y);
plot(y,abs((Lockwood-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%Cantrell 2004
k=74;
p=0.2410117;
f = p*(1+y)+((1-2*p)/(k+1))*sqrt((1+k*y).*(k+y));
Cantrell4 = (1+y)-(2-pi/2)*y./f;
%plot(y,abs((Cantrell4-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
%hold on;

%Sykora-Rivera 2005; renamed Cantrell's Particularly Fruitful
SykoraRivera = (pi*y+(1-y).*(1-y))./(1+y);
plot(y,abs(SykoraRivera-E)./E,'-k','LineWidth',2,'Color',[1 0 0]);
hold on;

%Maertens Roger's YNOT formula; a type of Holder mean
s = log(2)/log(pi/2);
ynot = (1+y.^s).^(1/s);
plot(y,abs((ynot-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;


% Padè combinations
%--------------------------------------------------------------------------------------
% To run these, keep enabled the Padè formuals in type K

%Selmer-Lindner, Padè11 with Padè20
d1 = (pi/4)*81/64-1;
d2 = (pi/4)*19/15-1;
p = d1/(d1-d2) %=7.30996434244961
comb = p*Selmer+(1-p)*Lindner;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

%Michon-Selmer, Padè12 with Padè11
d1 = (pi/4)*19/15-1;
d2 = (pi/4)*80/63-1;
p = d1/(d1-d2); %=2.07045704986602
comb = p*Michon+(1-p)*Selmer;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

%Hudson-Michon, Padè21 with Padè12
d1 = (pi/4)*80/63-1;
d2 = (pi/4)*61/48-1;
p = d1/(d1-d2); %=3.42546255957109
comb = p*Hudson+(1-p)*Michon;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

%JaconsenWaadeland-Hudson, Padè22 with Padè21
d1 = (pi/4)*61/48-1;
d2 = (pi/4)*187/147-1;
p = d1/(d1-d2); %=1.88647087966682
comb = p*Jacobsen+(1-p)*Hudson;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

%Padè32 with JaconsenWaadeland (Padè22)
d1 = (pi/4)*187/147-1;
d2 = (pi/4)*1573/1236-1;
p = d1/(d1-d2); %=2.07514774103763
comb = p*Pade32+(1-p)*Jacobsen;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

%Padè33 with Padè32
d1 = (pi/4)*1573/1236-1;
d2 = (pi/4)*47707/37479-1;
p = d1/(d1-d2); %=2.38364165550747
comb = p*Pade33+(1-p)*Pade32;
plot(y,abs((comb-E)./E),'-k','LineWidth',2,'Color',[0 0.5 0]);
hold on;

end
%--------------------------------------------------------------------------------------


%--------------------------------------------------------------------------------------
% Type C: Exact for y=0 and y=1, with crossings
%--------------------------------------------------------------------------------------
if typeC
    
%R.Bartolomeu 2004, 0.98%
t = (pi/4)*(1-y);
Bartolomeu2 = (pi/4)*sqrt(2*(1+y.*y)).*sin(t)./t;
plot(y,abs((Bartolomeu2-E)./E),'-k','LineWidth',2,'Color',[0 0 0.8]);
hold on;

%D.F.Rivera Asymmetric (y<=0.5), 0.45%
Rivera1 = 1+((pi-2)/2)*y.^1.456;
plot(y,abs((Rivera1-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0]);
hold on;

%D.F.Rivera 2004, 103 ppm
Rivera2 = (pi*y+(1-y).*(1-y))./(1+y)-(89/584)*((sqrt(y)-y)./(1+y)).^2;
plot(y,abs((Rivera2-E)./E),'-k','LineWidth',2,'Color',[0 0.7 0]);
hold on;

%Cantrell David 2001
%overall optimal, 83 ppm
s = 0.825056176207;
Cantrell1 = (1+y)-(2-pi/2)*y./(((1+y.^s)/2).^(1/s));
plot(y,abs((Cantrell1-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;
%round optimal
%s = (3*pi-8)/(8-2*pi);
%Cantrell2 = (1+y)-(2-pi/2)*y./(((1+y.^s)/2).^(1/s));
%plot(y,abs((Cantrell2-E)./E),'-k','LineWidth',2,'Color',[0 0.8 0]);
%hold on;
%elongated optimal
%s = log(2)/log(2/(4-pi));
%Cantrell3 = (1+y)-(2-pi/2)*y./(((1+y.^s)/2).^(1/s));
%plot(y,abs((Cantrell3-E)./E),'-k','LineWidth',2,'Color',[0 0.8 0]);
%hold on;

%Sykora 2005, 78.5 ppm
Sykora = (pi*y+(1-y).*(1-y))./(1+y)-y.*((1-y).^2)./(8*(1+y).*((1+y).^2+3.14159265*y));
plot(y,abs((Sykora-E)./E),'-k','LineWidth',2,'Color','red');
hold on;

%Cantrell-Ramanujan 2004, 14.5 ppm
CantrellRamanujan = (pi/4)*(1+y).*(1+3*h./(10+sqrt(4-3*h))+(4/pi-14/11)*h.^12);
plot(y,abs((CantrellRamanujan-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

%Cantrell2 2006, 5.2 ppm
p = 3.982901;
q = 66.71674;
s = 18.31287;
t = 23.39728;
r = 56.2007;
Cantrell2 = (1+y)-(y./(4*(1+y))).*((p*(1+y).^2 + q*y + r*(y./(1+y)).^2)./((1+y).^2 + s*y + t*(y./(1+y)).^2));
plot(y,abs((Cantrell2-E)./E),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;

end
%--------------------------------------------------------------------------------------

set(gcf,'Color','white');
set(gca,'FontSize',16);
set(gca,'LineWidth',2);
set(gca,'TickLength',[.03 .03]);
set(gca,'yScale','log');
set(gca,'xScale','log'); %Comment this line for linear x-scale
axis([0.0001 1 1e-7 1]); %Below 1e-8 the precision of the program is insufficient
%set(gca,'yTick',[-0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0]);
%set(gca,'xTick',[.2,.4,.6,.8,1]);