function C = pset4problem2c(T,M)
%This function evaluates the heat capacities at a given temperature given
%the size of the basis set M to use
% Input:
%       T: vector of Temperatures at which the value of heat capacity is required
%       M: The size of the basis set used
%Output:
%       C: The vector of heat capacity at each temperature in the vector T

%For a given value of M we can calculate 2*M+1 values of energies
% we use problem2b to solve for the energies

I = 3e-45; %the moment of inertia Kg-m2
kb = 1.38e-23; %boltzmann constant

%calculate the energies of the system
E = calculateE(M,I);

%loop over all the temperatures and calculate the values of C for each
%temperature and store it in the vector C.
for k=1:length(T)
    E_avg = 0;
    E2_avg = 0;
    denominator = 0;

    for n=1:length(E)
        E_avg = E_avg + E(n)*exp(-E(n)/kb/T(k));
        E2_avg = E2_avg + E(n)^2*exp(-E(n)/kb/T(k));
        denominator = denominator + exp(-E(n)/kb/T(k));
    end
    C(k) = (E2_avg/denominator - (E_avg/denominator)^2)/kb/T(k)^2;
end

%convert the units of C from J/molecule/K to J/mol/K by multiplying by
%Avagadro number
C = C*6.023e23; 
return;

%**********************************************************************

function E = calculateE(M,I)
%calculates the energies of the system given the size of the basis function
% written by Sandeep Sharma on 26th Sep 2006
%Taken input 
%       M which is the size of the basis function to be used
%       I the angular moment of the hindered rotor
%outputs the vector of energies.

%Start by defining all the constants
h = 6.625e-34;  %J-s


%get the coefficients ys of the potential function
%y = [2.0667 -2.0333 1.2333 -1.7667 0.5]'*1e-19;
y = [2.067 -2.033 1.056 -1.767 0.678]'*1e-20;

%fill the elements due to the differential operators
for n=1:2*M+1
    H(n,n) = h^2*(n-M-1)^2/8/pi^2/I;
end

%fill in the elements for on the upper triagular side and diagonal of the matrix.

for n=1:2*M-3
    H(n,n) = H(n,n) + y(1);
    H(n,n+1) = y(2)/2;
    H(n,n+2) = y(3)/2;
    H(n,n+3) = y(4)/2;
    H(n,n+4) = y(5)/2;
end
H(2*M-2,2*M-2) = H(2*M-2,2*M-2) + y(1);
H(2*M-2,2*M-1) = y(2)/2;
H(2*M-2,2*M)  = y(3)/2;
H(2*M-2,2*M+1) = y(4)/2;

H(2*M-1,2*M-1) = H(2*M-1,2*M-1) + y(1);
H(2*M-1,2*M) = y(2)/2;
H(2*M-1,2*M+1) = y(3)/2;

H(2*M,2*M) = H(2*M,2*M) + y(1);
H(2*M,2*M+1) = y(2)/2;

H(2*M+1,2*M+1) = H(2*M+1,2*M+1) + y(1);

%fill in the elements on the lower triangular part of the matrix.

for n= 5:2*M+1
    H(n,n-1) = y(2)/2;
    H(n,n-2) = y(3)/2;
    H(n,n-3) = y(4)/2;
    H(n,n-4) = y(5)/2;
end
H(2,1) = y(2)/2;

H(3,2) = y(2)/2;
H(3,1) = y(3)/2;

H(4,3) = y(2)/2;
H(4,2) = y(3)/2;
H(4,1) = y(4)/2;

%The eigenvalues of matrix H are the allowed energy states of the system
E = eig(H);


return;