% chempot uses function fermi.m
% carrier density n(cm-3), temperature kelvin(K), relative error rerr
% returns chemical potential mu1 in meV
function [chempot]=chempot(kelvin,m1,n,rerr)
%   n=1.e18;
%   kelvin=300.0;
%   rerr=1.e-3;

   m0=9.10956;
%   m1=0.5;
	echarge=1.60219;
	hbar=1.0545928;
   kB=8.617e-5;
   
   kF1=(3.0*(pi^2.0)*n)^0.33333333;
   eF=1.e-11*((hbar*kF1)^2)/(2.0*m0*m1*echarge);
   kBT=1000.*kelvin*kB;
   beta=1./kBT;
   
	mumax=eF;
	x=((n*1.e6)/2.)*(((2.e-15*pi*beta*(hbar^2))/(echarge*m0*m1))^1.5);
	mumin=(+1./beta)*log(x);
   
   emax=eF+(15./beta);
   de=emax/1000.;
	const=1.e16*((20.*echarge*m1*m0)^.5)*echarge*m1*m0/((pi^2)*(hbar^3));
   
for j=1:25
   mu1=mumin+((mumax-mumin)/2.);
	energy=0.0;
	ainter=0.0;
   
%	calculate carrier density n'
	for i=1:1000;
		energy=energy+de;
		ainter=ainter+(((sqrt(energy))*de)*fermi(beta,energy,mu1));
   end;
   nprime=const*ainter;
%  delta is relative error in carrier density
   delta=(n-nprime)/n;
   	if((abs(delta)) < rerr)
      	break;
     			elseif(delta < 0.) 
         	mumax=mu1;
			else 
			mumin=mu1;
   	end;
end;
%	print output: chemical potential mu1(meV), Fermi energy eF(meV), 
%   temperature kelvin(K) 
%	kBT (meV), carrier density nprime(cm-3), number of iterations j
%mu1
%eF
%kelvin
%kBT
%nprime
%j
if j >= 25
   'check convergence!'
end;
chempot=mu1;
return

%	compare with second-order mu2(meV) and fourth-order mu4(meV) expansion
%mu2=eF-(((pi*kBT)^2)/(12.*eF))
%mu4=eF-(((pi*kBT)^2)/(12.*eF))-(7*((pi*kBT)^4)/(960.*(eF^3)))

