% Tarea de termodinamica
% Problemo seis punto treinta
% Oscar Zambrano y Benjamin Hatchett

epsilon = 0.0076; % eV
K = 8.617e-5; %ev/K
T = 1:.1:2000;
j = 0:2:20;
jodd = 1:2:21;
beta = 1./(K * T);

%Here is parahydrogen (Even integers for j, j=0,2,4,..)
z = NaN(length(T),length(j));
clear p;

for k = 1:length(T)
    clear Z;
    for i = 1:length(j);
        
        a = (2*j(i))+1;
        b = (j(i)*(j(i)+1));
        
        Z(i) = (a*exp(-((b*epsilon)*(beta(k)))));
        
    end %i--> j loop
    
    p(k) = sum(Z);
    z(k,:) = Z(:);
    
end % k--> temperature loop

% Now orthohydrogen (Odd integers for j, j=1,3,5,...)
clear po;
zo = NaN(length(T),length(jodd));
for k = 1:length(T)
    clear Z;
    for i = 1:length(jodd);
        
        a = (2*jodd(i))+1;
        b = (jodd(i)*(jodd(i)+1));
        % Triplet state so multiply by 3 to account for spin states
        Z(i) = 3*(a*exp(-((b*epsilon)*(beta(k)))));
    end %i--> j loop
    
    po(k) = sum(Z);
    zo(k,:) = Z(:);
    
end % k--> temperature loop

% calculate ln Z so that partial Z / partial E can be represented as
% partial ln Z / partial E since d/dx ln x = 1/x
ln_p = log(p);
ln_po = log(po); 

clear Eeven Ceven Eodd Codd;

% Using forward differencing like Mac Dre
for i = 1 : (length(T)-1)
    Eeven(i) = -(ln_p(i+1)-ln_p(i))/(beta(i+1)-beta(i));
    Eodd(i) = -(ln_po(i+1)-ln_po(i))/(beta(i+1)-beta(i));
end;

for i = 1 : (length(T)-2)
    Ceven(i) = (Eeven(i+1)-Eeven(i))/(T(i+1)-T(i));
    Codd(i) = (Eodd(i+1)-Eodd(i))/(T(i+1)-T(i));
end

KTeps = (K*T)/epsilon;

Cvkeven = Ceven(:)./K;
Cvkodd = Codd(:)./K;

%Finally we do the mixture of 1/4 parahydrogen and 3/4 orthohydrogen
Cvkmix = .75*(Cvkodd(:)) + .25*(Cvkeven(:));

figure;
plot(KTeps(1:6000),Cvkeven(1:6000),'b','linewidth',2)
hold on
plot(KTeps(1:6000),Cvkodd(1:6000),'k','linewidth',2)
plot(KTeps(1:6000),Cvkmix(1:6000),'r','linewidth',2);

title('Rotational Heat Capacities of the Hydrogen Atom','fontsize',14,'fontweight','bold');
ylabel('C_v / k','fontsize',14,'fontweight','bold');
xlabel('kT / \epsilon','fontsize',14,'fontweight','bold');
set(gca,'fontsize',12,'fontweight','bold','linewidth',1.5)
axis([0 7 0 1.5]);
legend('Parahydrogen','Orthohydrogen','3:1 Mixture','location','southeast');



% plotting partition functions
figure;
bar(j,Z);
hold on; 
plot(j,Z,'k','linewidth',5);
axis([-0.5 14 0 11]);

xlabel('j','fontsize',14,'fontweight','bold');
ylabel('Z(j)','fontsize',14,'fontweight','bold');
title('Orthohydrogen Partition Sums','fontsize',14,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold');