function [fermi]=fermi(beta,energy,mu1)
%	fermi is the Fermi-Dirac function
%   beta=1/kBT and mu1 is chemical potential
	x=(energy-mu1)*beta;
   if(x > 180.0)
      x=180.;
      end;
    if(x < -180.);
       x=-180.;
       end;
	fermi=1./((exp(x))+1.);
	return;