To Spherical People,
In many calculations it is necessary to obtain spherical harmonic functions for a large range of l. For example, the Mars topography data set can include harmonic degrees up to 100. If you pass a given l and m to the function plgndr in Numerical Recipes the result is just fine for l up to about 16 or 17. Of course, this result is ordinarily normalized by a factor so that it is orthonormal over the sphere. The problem is that for l =18 and higher the numbers involved (for m close to l) in plgndr become unmanageably large for straightforward computation. Trying to normalize these numbers with the factor after applying the routine gives very inaccurate results. Another source of inaccuracy is in the recursion used to calculate P(l,m,x) for l>m>1, which requires taking the difference of very large numbers, leading to errors attributable to round off.
After trying several different approaches, I have solved this problem by modifying the recursions used in plgndr to use the normalization factors during the calculations themselves at all levels. This ensures that the calculations will only involve numbers of the same order of magnitude as the spherical harmonic function that is being calculated. As a result, this type of calculation can give accurate and stable results for values of l well into the hundreds. It is actually a fairly simple fix, just multiply the recursions in the routine by the appropriate factors for the l=m part and the l>m part.
I do not know the rules for posting code in these forums, so I will abstain from posting my subroutine here. If you want to know how this is done or have any questions on this topic, then feel free to e-mail me at hernlund@ess.ucla.edu.

Cheers!
John Hernlund

Hi John,

Yes, plgndr is on our list of routines to be upgraded. In the meantime, attached is a gif file of the renormalized functions and the recurrence relation, and below is a Fortran 77 implementation. (Sorry for the loss of tabs)

Saul Teukolsky

FUNCTION plgndr(l,m,x)
INTEGER l,m
REAL plgndr,x
INTEGER i,ll
REAL fact,oldfact,pll,pmm,pmmp1,omx2,PI
PARAMETER (PI=3.1415927)
if(m.lt.0.or.m.gt.l.or.abs(x).gt.1.)pause 'bad arguments in plgndr'
pmm=1.
if(m.gt.0) then
omx2=(1.-x)*(1.+x)
fact=1.
do i=1,m
pmm=pmm*omx2*fact/(fact+1.)
fact=fact+2.
enddo
endif
pmm=sqrt((2*m+1)*pmm/(4.*PI))
if(mod(m,2).eq.1)pmm=-pmm
if(l.eq.m) then
plgndr=pmm
else
pmmp1=x*sqrt(2.*m+3.)*pmm
if(l.eq.m+1) then
plgndr=pmmp1
else
oldfact=sqrt(2.*m+3.)
do ll=m+2,l
fact=sqrt((4.*ll**2-1.)/(ll**2-m**2))
pll=(x*pmmp1-pmm/oldfact)*fact
oldfact=fact
pmm=pmmp1
pmmp1=pll
enddo
plgndr=pll
endif
endif
return
END

I wondered if you have a similar function (to plgndr, with the spherical harmonic multipliers included) that returns the first derivative with respect to the colatitudinal argument of the spherical harmonic, usually the first angular argument.

My application is a small graphics program that needs the partials to compute the normal vector to a surface, defined as a spherical harmonic added to a radius.

Thanks in advance,

Don Ensley

Dear Don,

I don't have such a program. But it should be easy for you to do yourself using a recurrence relation to express the derivative in terms of 2 functions with l+1 and l-1.

Saul Teukolsky