%para graficar algunas cosas del problema 2 de la guia de Bose,
%%2C 2014/5/6/7

%defino algunas funciones 

%logaritmo de la funcion Gran Canonica
logZGC=@(z,T,V)(-log(1-z) + 2*V.*T.^(3/2).*log(2./(2-z)));

%de la derivacion analitica <N> = z d logZGC / dz
Nm=@(z,T,V)(z./(1-z) + 2*V.*T.^(3/2).*z./(2-z));
%separo explicitamente N0 y Nex

N0=@(z)(z./(1-z));
Nex=@(z,T,V)(2*V.*T.^(3/2).*z./(2-z));


%analiticamente encuentro que v T^(3/2) = 1/2 o Tc=1/(2v)^(2/3) me define la
%condicion critica, cuando el numero de particulas es igual al de
%excitados.
Tc=@(v)(1./(2*v).^(2./3.));



%La solucion analitica de <N> = N me determina z(N,T,V), esta es la solucion de la cuadratica

z=@(N,T,V)((2+3*N+2*T.^(3/2).*V-sqrt((N+2).^2-4.*(N-2).*T.^(3/2).*V+4.*T.^3.*V.^2))./(2*(1+N+2.*T.^(3/2).*V)));



%Tomo kB=1
%Presion, directamente del Gran Canonico
Pv=@(z,T,V)(T.*logZGC(z,T,V)./Nm(z,T,V));


%Energia Interna, derivando analiticamente logZGC)
U=@(z,T,V)(3*T.^(5/2).*V.*log(2./(2 - z)));



%EJEMPLOS:
%plots de la fraccion condensada, N0/N

v=1; % analizo solo este valor de v=V/N
T=linspace(0,3*Tc(v),100); % preparo un rango de temperaturas

figure(1)

clf;
for logN=0:8;
NN=10^logN;
subplot(2,1,1)
        plot(T./Tc(v),N0(z(NN,T,v*NN))./NN,'b-',T./Tc(v),Nex(z(NN,T,v*NN),T,v*NN)./NN,'r-');
hold on
subplot(2,1,2)
        plot(T./Tc(v),z(NN,T,v*NN));
hold on
end 
subplot(2,1,1)
xlabel('T/Tc');ylabel('Ni/N');
legend('N0/N','N_{exc}/N')
title('Fraccion condensada, excitaciones y fugacidad para N=1,10,100,..10^8')
subplot(2,1,2)
xlabel('T/Tc');ylabel('z');


%NB: aca estoy recalculando en cada plot z, pero deberia
%calcularlo al principio y re-usarlo, para hacerlo mas
%eficiente

%repito para la presion
figure(2)
hold off;
for logN=0:8;
NN=10^logN;
plot(T./Tc(v),Pv(z(NN,T,v*NN),T,v*NN));
hold on
end 
xlabel('T/Tc');ylabel('Pv');
title('Pv para  N=1,10,100,..10^8')


%Ahora la energia interna por particula
figure(3)
hold off;
for logN=0:8;
NN=10^logN;
plot(T./Tc(v),U(z(NN,T,v*NN),T,v*NN)./NN);
hold on
end 
xlabel('T/Tc');ylabel('U/N');
title('U/N para  N=1,10,100,..10^8')

%Ahora comparo energia interna por particula con 3/2 Pv, son
%iguales?
figure(4)
for logN=0:2:7;
NN=10^logN;
subplot(4,1,logN/2+1);
plot(T./Tc(v),U(z(NN,T,v*NN),T,v*NN)./NN,'r-',T./Tc(v),3/2*Pv(z(NN,T,v*NN),T,v*NN),'bo');
title(['N = ' num2str(NN)]);
legend('U/N','3/2Pv')
ylabel('U/N y 3/2 Pv');
end 
xlabel('T/Tc');



%Ahora el Cv (por particula) 
figure(5)
hold off;
for logN=0:8;
NN=10^logN;
%hago el calculo numericamente, la expresion analitica es un poco larga...
%NB: Cv es a N (no z) constante.
UiN=U(z(NN,T,v*NN),T,v*NN)./NN;
Cvi=gradient(UiN,T);
plot(T./Tc(v),Cvi);
hold on
end 
xlabel('T/Tc');ylabel('Cv');
title('Cv para N=1,10,100,..10^8')

%NB: todas las derivadas podrian hacerlas numericamente.