function [LPOL] = calc_legendre(l, m, x)


%   calc_legendre.m
%
%   This will compute the associated Legendre polynomial
% 
%          P  m (x)
%             l 
%
%   m and l are integers satisfying :  0 <= m <= l, while x lies in the range -1 <= x <= 1
%
%   x == cos(phi)

%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.   


% Check that input variables are correct
if ((m < 0) | (m > l) | (abs(x) > 1.0))
	error('Bad arguments in routine calc_legendre')
end


pmm = 1.0;
if (m > 0)
	somx2 = sqrt((1.0 - x) * (1.0 + x));
	fact = 1.0;
	for i = 1: m
		pmm = pmm * (-1.0 * fact * somx2);
		fact = fact + 2.0;
	end
end

if (l == m)
	LPOL = pmm;
	return;
else
	pmmp1 = x * (2 * m + 1) * pmm;
	if (l == (m + 1))
		LPOL = pmmp1;
		return;
	else
		for ll = (m + 2):l,
			pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
			pmm = pmmp1;
			pmmp1 = pll;
		end
		LPOL = pll;
		return;
	end
end
LPOL=LPOL
return;

