function Cv = DebyeCv(t)
%DebyeCV: returns CV, heat capacity, per particle, as  function of t=theta_D/T
%in the Debye model.
%Cv should be mulitplied by 3*k_B to get the heat capacity in 
%physical units of erg/Kelvin

%formula for Cv from Pathria, page 196, section 7.3
    if (t <300)
        fun=@(x)((x.^4.*exp(x))./(exp(x)-1).^2);
        q=integral(fun,eps,t);
        Cv=(3/t.^3).*q;
    else
        Cv=((4*pi^4)/5)*t.^(-3);    %switch to low temp approximation
    end
    
end