% CalulateChemicalPotential.m
% Input:  Doping, Temperature (m^-3), (K)
% Output: Chemical potential (eV)
%**************************************************************************
function [Mu] = CalculateChemicalPotential(Temp, NDoping, m, EFermi);

hbar = 1.05457159e-34;                       %Planck's constant (Js)
echarge = 1.6021764e-19;                     %Electron charge (C)
kBoltz = 1.3806505e-23;                      %Boltzmann constant (eV/K)
const = (2.0*m/(hbar*hbar)).^1.5/(2.0*pi*pi); 

if abs(Temp) > 1.0e-10
    beta = 1.0/(kBoltz*Temp);
else
    beta = 1.0/(kBoltz*1.0e-10);
end

MuMin = kBoltz/echarge*Temp*log(NDoping/2.0*...
    ((2.0*pi*hbar*hbar)/m*beta).^1.5);
MuMax = EFermi;

EMin = 0.0;
%Energy cutoff for temperature T (K)
EMax = EFermi + 15.0*kBoltz*Temp/echarge; 

nSteps = 10000; %Number of energy points in integral
EStep = (EMin+EMax)/nSteps;                 
rerr = 1.0;
nprime = 0.0;
E = EMin;
%**************************************************************************
% Applied Quantum Mechanics, AFJ Levi, p.337
% (1) Start with an initial guess of the chemical potential 
%     "Mu = MuMin + (MuMax-MuMin)/2.0;".
% (2) Calculate the corresponding number of occupied states "nprime"
%     by performing the energy integral over the Fermi function. 
% (3) Compare "nprime" with the given doping concentration "NDoping"
%     and choose a new "Mu"
% (4) Go back to (2), iterate until integral over Fermi function gives 
%     "NDoping" up to error "rerr"
%**************************************************************************
if Temp < 1.0e-10
    Mu = EFermi;
else    
    while rerr > 0.0001
        Mu = MuMin + (MuMax-MuMin)/2.0;
        nprime = 0.0;
        E = EMin;
        for j = 0 : nSteps
 nprime = nprime + sqrt(E*echarge)*1.0/(1.0+exp(-beta*echarge*(Mu-E)));        
 E = E + EStep;
        end
        nprime = const*nprime*EStep*echarge;
        if nprime > NDoping
            MuMax = Mu;
        else
            MuMin = Mu;
            rerr = abs(NDoping-nprime)/NDoping;
        end
    end
end