function Y_lm=Y_lm(l,m,theta,varphi)
%Y_lm=Y_lm(l,m,theta,varphi)
%(c) Moo K. Chung 
% mkchung@wisc.edu
% Update history 2005; July 2007; Mar 2008; 
%-----------------------------------------


Y=(-1)^m*legendre(l,cos(theta));	
Y=Y(abs(m)+1,:);
Ylm=Y;

if m>0 
	Ylm=Ylm.*cos(abs(m)*varphi);
end;
if m<0
	Ylm=Ylm.*sin(abs(m)*varphi);
end;
if m==0
    Ylm=Ylm.*cos(abs(m)*varphi)/sqrt(2);
end;

% old line 
%clm=sqrt((2*l+1)*factorial(l-abs(m))/(2*pi*factorial(l+abs(m))));
% factorial routine is very time consuming

for k = 0:(2*l)
    fact(k+1) = factorial(k);
end
clm = sqrt(((2*l+1)/(2*pi))*(fact(l-abs(m)+1)./fact(l+abs(m)+1)));
Y_lm=clm*Ylm;
	