clear all;

hbar = 1.0546e-34;

clight = 3e8;

echarge = 1.6e-19;

kb = 1.38e-23/1.6e-19; %kb in eV

me = 9.11e-31; %electron mass in kg

Te = logspace(2,3,100);



%****************** Finding the allowable energy levels ****************

% Actually, I'm not finding the levels themselves - (l^2+m^2+n^2)

% is only proportional to the energy levels, but thats good enough 

% to determine degeneracy. It is not necessary to determine degeneracy 

% to solvethis problem.

index = 1;

max = 25;

for l = 1:max

    for m = 1:max

        for n = 1:max

            E(index) = (l^2+m^2+n^2);

            index = index+1;

        end

    end

end



E = sort(E');



energy = unique(E);



%****************** Finding the degenaracy of each level **********

index = 1;

degen = zeros(size(energy,1),1);

for l = 1:size(E,1)-1

    degen(index) = degen(index)+1;

    if (E(l) ~= E(l+1))

        index = index+1;

    end

end

  

%************** quantum dot specifications ********

L = 20e-10; % quantum dot of 20 A 

%L = 100e-10;

Nelects = 2e28*L^3;

energy = (1/(2*me))*(pi*hbar/L)^2*energy; %energy in joules

energy = energy/echarge; %energy in eV

%mu = [0:0.0005:5];



for Tindex = 1:size(Te, 2)

    T = Te(Tindex);

    % Initially, I used this crude method (graphical) to determine mu. mu should be 

    % calculated faiirly accurately to determine energy density and specific heat.

    % Hence I stopped using this method and used a bisection method to

    % determine mu.

    %     for index = 1:size(mu,2)

    %         Ne(index) = 2*sum(degen./(exp((energy-mu(index))/(kb*T))+1));

    %         if (Ne(index) > Nelects), break, end;

    %     end

    %     Ne(index:size(mu,2)) = 0;

    %     [best, bestindex] = min(abs(Ne-Nelects));

    %     chempot(Tindex) = mu(bestindex);

    

    % This is the bisection method to determine mu

    mulow = 0.0;

    muhigh = 50.0;

    f = 10;

    while (f > 1e-10)

        mumid = (mulow+muhigh)/2;

        flow = (Nelects-2*sum(degen./(exp((energy-mulow)/(kb*T))+1)));

        fhigh = (Nelects-2*sum(degen./(exp((energy-muhigh)/(kb*T))+1)));

        fmid = (Nelects-2*sum(degen./(exp((energy-mumid)/(kb*T))+1)));

        f = abs(fmid);

        if (fmid*flow > 0)

            mulow = mumid;

        else

            muhigh = mumid;

        end

    end

    chempot(Tindex) = mumid; % chemical potentia as a function of temperature

    error(Tindex) = abs(fmid/Nelects); % a measure of error

    dmudT(Tindex) = -(1/T)*...

        sum(degen.*(energy-mumid)./cosh((energy-mumid)/(2*kb*T)).^2)./...

        sum(degen./cosh((energy-mumid)/(2*kb*T)).^2); % derivative of mu w.r.t temperature

    

    totenergy(Tindex) = sum(degen.*energy./...

        (exp((energy-mumid)/(kb*T))+1)); %total energy as a function of temperature

    spheat(Tindex) = (1/(kb*T))*sum(degen.*energy./...

        (cosh((energy-mumid)/(2*kb*T)).^2).*((energy-mumid)/T+dmudT(Tindex)));

    %specific heat as a function of temperature

        

end



color = 'b';



subplot(2,2,1)

plot(Te, chempot, color)

ylabel('\mu (eV)')

hold on



subplot(2,2,2)

plot(Te, dmudT, color)

ylabel('d\mu/dT (eV K^{-1})')

hold on



subplot(2,2,3)

plot(Te, totenergy, color)

xlabel('Temperature (K)')

ylabel('eV')

hold on



subplot(2,2,4)

loglog(Te, spheat*echarge/L^3, color)

xlabel('Temperature (K)')

ylabel('JK^{-1})m^{-3}')

hold on

%plot(E)