double plgndr(double l, double m, double x)
{
	int ll,i;
	double pmm,somx2,fact,temp,pmmp1,pll;
    if( m < 0 || m > l || fabs(x) > 1.0)
		cout<<"bad arguments";
    pmm = 1.0;
    if (m > 0)
	{
        somx2 = sqrt((1.0 - x) * (1.0 + x));
        fact = 1.0;
        for( i = 1; i<=m; i++)
		{
            pmm = -pmm * fact * somx2;
            fact = fact + 2.0;
        }
    }
    if (l == m)
	{
        temp = pmm;
	}
    else
	{
        pmmp1 = x * (2 * m + 1) * pmm;
        if (l == m + 1)
		{
            temp = pmmp1;
		}
        else
		{
            for (ll = m + 2; ll<=l; ll++)
			{
                pll = x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm;
                pll = pll / (ll - m);
                pmm = pmmp1;
                pmmp1 = pll;
            }
            temp = pll;
        }
    }
	return temp;
}