function Y = real_sh(theta,phi,l,m)
P = legendre(l,cos(theta));
P = P(abs(m)+1,:);

[x,y]=size(theta);
if y==1
    P=P';
end


Y=0;
if m<0
    Y = sqrt(2)* sqrt((2*l+1)*factorial(l-abs(m))/(4*pi*factorial(l+abs(m)))) *P.*sin(abs(m)*phi);
end
if m==0
    Y = sqrt((2*l+1)/(4*pi)) *P;
end

if m>0
    Y = sqrt(2)* sqrt((2*l+1)*factorial(l-m)/(4*pi*factorial(l+m))) *P.*cos(m*phi);
end
