function [mu]=mu(emass,ncarrier,kelvin,rerr)
% mu uses function fermi.m
% carrier density n(m-3), temperature kelvin(K)
% returns chemical potential mu in eV

	echarge=1.60219e-19;				%Electron charge C
	hbar=1.0545928e-34;				%Planck's constant J s
   kB=8.617e-5;						%Boltzmann constant eV K-1
   
   kF1=(3*(pi^2)*ncarrier)^(1/3);							%Fermi wave vector m-1
   eF=((hbar*kF1)^2)/(2*emass*echarge);		%Fermi energy eV
   kBT=kelvin*kB;										%Thermal energy eV
   beta=1./kBT;
   
	mumax=eF;
	x=(ncarrier/2.)*(((2*pi*beta*(hbar^2))/(echarge*emass))^1.5);
	mumin=(+1./beta)*log(x);
   
   emax=eF+(15./beta);
   de=emax/1000.;
	const=((2*echarge*emass)^.5)*echarge*emass/((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=(ncarrier-nprime)/ncarrier;
   	if((abs(delta)) < rerr)
      	break;
     			elseif(delta < 0.) 
         	mumax=mu1;
			else 
			mumin=mu1;
   	end;
end;
			if j >= 25
   			'check convergence!'
			end;
mu=mu1;
return;