%% IDEAL BOSE GAS IN A 3D BOX
% mu=mu/Tc as a function of t=T/Tc
clc
clear all

% preliminary check for integrals
z=0.8;
num = quadgk(@(x) int_num(x,z),0,inf)

z=1;
den = quadgk(@(x) int_num(x,z),0,inf)

num/den

%%
% checks for find zero
N=100; t=0.2;

musol=fzero(@(mu) fun_zero(N, mu, t), [-100,-0.00001]);
musol

%%
% plot the chemical potential
Npts=100;
t=linspace(0.1,2,Npts);

N=100;
for i=1:Npts
    musol_100(i)=fzero(@(mu) fun_zero(N, mu, t(i)), [-100,-0.00001]);
end

N=1000;
for i=1:Npts
    musol_1000(i)=fzero(@(mu) fun_zero(N, mu, t(i)), [-100,-0.00001]);
end

N=4000;
for i=1:Npts
    musol_4000(i)=fzero(@(mu) fun_zero(N, mu, t(i)), [-100,-0.00001]);
end

figure(1)
hold on
plot(t,musol_100,'r-','LineWidth',2)
plot(t,musol_1000,'b-','LineWidth',2)
plot(t,musol_4000,'g-','LineWidth',2)
legend('N=100','N=1000', 'N=4000','Location','South')
xlabel('T/T_c','fontsize',16),ylabel('mu/T_c','fontsize',16)
hold off

%%
% plot the condesate fraction

n0_100=(1/100)*1./(exp(-musol_100./t)-1);
n0_1000=(1/1000)*1./(exp(-musol_1000./t)-1);
n0_4000=(1/4000)*1./(exp(-musol_4000./t)-1);

figure(2)
hold on
plot(t,n0_100,'r-','LineWidth',2)
plot(t,n0_1000,'b-','LineWidth',2)
plot(t,n0_4000,'g-','LineWidth',2)
legend('N=100','N=1000', 'N=4000','Location','North')
xlabel('T/T_c','fontsize',16),ylabel('n_0','fontsize',16)
hold off
