function Pnm = legendreFunctions(theta, maxDegree)

global legendreFactor1;
global legendreFactor2;

% init legendre factors at the first time
if (size(legendreFactor1,1) < maxDegree+1)
  legendreFactor1 = zeros(maxDegree+1, maxDegree+1);
  legendreFactor2 = zeros(maxDegree+1, maxDegree+1);
  legendreFactor1(2,2) = sqrt(3);
  for n=2:maxDegree
    legendreFactor1(n+1,n+1) = sqrt((2*n+1)/(2*n));
  end
  for m=0:maxDegree-1
    for n=m+1:maxDegree
      f=(2*n+1)/((n+m)*(n-m));
      legendreFactor1(n+1,m+1) = sqrt(f*(2*n-1));
      legendreFactor2(n+1,m+1) = sqrt(f*(n-m-1)*(n+m-1)/(2*n-3));
    end
  end
end

cosTheta = cos(theta);
sinTheta = sin(theta);

Pnm = zeros(maxDegree+1, maxDegree+1);
Pnm(1,1)=1;
for n=1:maxDegree
  Pnm(n+1,n+1) = legendreFactor1(n+1,n+1)*sinTheta*Pnm(n,n);
end
for m=0:maxDegree-1
  Pnm(m+2,m+1) = legendreFactor1(m+2,m+1)*cosTheta*Pnm(m+1,m+1);
end
for m=0:maxDegree-1
  for n=m+2:maxDegree
    Pnm(n+1,m+1) = legendreFactor1(n+1,m+1)*cosTheta*Pnm(n,m+1)-legendreFactor2(n+1,m+1)*Pnm(n-1,m+1);
  end
end
