%--------------------------------------------------------------------------------------
% Gamma Function Approximations
%--------------------------------------------------------------------------------------
% Used to compare Stirling, Laplace, Ramanujan, and Nemes formulae
clear;clc;clf;echo off;

% Specify what you want to plot:
plotcrude = 1; % 0-th order Stirling
plotstirl = 1; % Stirling 1,3,5
plotraman = 1; % Ramanujan 1,2,3
plotlapla = 1; % Laplace 1,2,3
plotnemes = 1; % Nemes 2,4,6
plotclosed= 1; % Nemes, closed formula

e = exp(1);

n = 100;
z = linspace(0.5,0.5*n,n);
lngammaz = 0*z;
lngammaz(1) = log(sqrt(pi));
lngammaz(2) = 0;
for k=3:n
  lngammaz(k) = lngammaz(k-2)+log(k/2.0-1);  
end

if plotcrude
% Crude
crd = z.*log(z/e)+0.5*log(2*pi./z);
plot(z,-(crd-lngammaz),'-k','LineWidth',2,'Color',[0 0 0]);
hold on;
end

if plotstirl
% Stirling
stirling1 = crd + (1/12)./z;
plot(z,+(stirling1-lngammaz),'-k','LineWidth',4,'Color',[0 0 1]);
hold on;
stirling3 = stirling1 - (1/360)./(z.^3);
plot(z,-(stirling3-lngammaz),'-k','LineWidth',2,'Color',[0 0 1]);
hold on;
stirling5 = stirling3 + (1/1260)./(z.^5);
plot(z,(stirling5-lngammaz),'-k','LineWidth',4,'Color',[0 0 1]);
hold on;
end

if plotraman
% Ramanujan
ramanujan1 = crd + (1/6)*log(1+(1/2)./z);
plot(z,-(ramanujan1-lngammaz),'-k','LineWidth',2,'Color',[0.8 0.5 0]);
hold on;
ramanujan2 = crd + (1/6)*log(1+(1/2)./z+(1/8)./(z.^2));
plot(z,-(ramanujan2-lngammaz),'-k','LineWidth',2,'Color',[0.8 0.5 0]);
hold on;
ramanujan3 = crd + (1/6)*log(1+(1/2)./z+(1/8)./(z.^2)+(1/240)./(z.^3));
plot(z,(ramanujan3-lngammaz),'-k','LineWidth',4,'Color',[0.8 0.5 0]);
hold on;
end

if plotlapla
% Laplace
laplace1 = crd + log(1+(1/12)./z);
plot(z,-(laplace1-lngammaz),'-k','LineWidth',2,'Color',[0 0.8 0]);
hold on;
laplace2 = crd + log(1+(1/12)./z+(1/288)./(z.^2));
plot(z,(laplace2-lngammaz),'-k','LineWidth',4,'Color',[0 0.8 0]);
hold on;
laplace3 = crd + log(1+(1/12)./z+(1/288)./(z.^2)-(139/51840)./(z.^3));
plot(z,(laplace3-lngammaz),'-k','LineWidth',4,'Color',[0 0.8 0]);
hold on;
plot(z,-(laplace3-lngammaz),'-k','LineWidth',2,'Color',[0 0.8 0]);
hold on;
end

if plotnemes
% Nemes
nemes2 = crd + z.*log((1+(1/12)./(z.^2)));
plot(z,-(nemes2-lngammaz),'-k','LineWidth',2,'Color',[1 0 0]);
hold on;
nemes4 = crd + z.*log((1+(1/12)./z./z)+(1/1440)./(z.^4));
plot(z,-(nemes4-lngammaz),'-k','LineWidth',2,'Color',[1 0 0]);
hold on;
nemes6 = crd + z.*log((1+(1/12)./z./z)+(1/1440)./(z.^4)+(239/362880)./(z.^6));
plot(z,(nemes6-lngammaz),'-k','LineWidth',4,'Color',[1 0 0]);
hold on;
end

if plotclosed
% Nemes Closed
nemesc = crd + (5/4)*z.*log((1+(1/15)./(z.^2)));
plot(z,-(nemesc-lngammaz),'-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',[.01 .01]);
set(gca,'yScale','log');
axis([0 n/2 1e-12 1]);