real function sph_harm(l,m,theta)
c Returns part of the spherical harmonic function 
c Y(l,m,theta) where l is the harmonic degree, 0=<m=<l,
c and theta is the colatitude angle. Unlike the calculation 
c for the associated Legendre function P(l,m,x), this routine 
c provides a normalization during the recursion calculations 
c that keeps the values at computationally viable magnitudes for 
c very large l.

	integer l,m,i,ll
	double precision x,fact,pll,pmm,pmmp1,somx2,fourpi,flm1,flm2
	real phi,theta

	x=cos(1.0*theta)
	fourpi=16.0*atan(1.0)
	pmm=1.0
	if(m.gt.0) then
	  somx2=sqrt((1.0-x)*(1.0+x))
	  fact=1.0
	  do i=1,m
	    pmm=-pmm*fact*somx2/sqrt(1.0*(i*(m+i)))
	    fact=fact+2.0
	  enddo
	endif
	if(l.eq.m) then
	  sph_harm=1.0*pmm
	else
	  pmmp1=x*(2*m+1)*pmm/sqrt(1.0*((m+m+1)))
	  if(l.eq.m+1) then
	    sph_harm=1.0*pmmp1
	  elseif(m.eq.0) then
	    do ll=m+2,l
	      pll=(x*(2*ll-1)*pmmp1-(ll-1)*pmm)/ll
	      pmm=pmmp1
	      pmmp1=pll
	    enddo
	    sph_harm=pll
	  else
	    do ll=m+2,l
	      flm1=sqrt(1.0*(ll-m))/sqrt(1.0*(ll+m))
	      flm2=flm1*sqrt(1.0*(ll-m-1))/sqrt(1.0*(ll+m-1))
	      pll=(x*(2*ll-1)*pmmp1*flm1-(ll+m-1)*pmm*flm2)/(ll-m)
	      pmm=pmmp1
	      pmmp1=pll
	    enddo
	    sph_harm=1.0*pll
	  endif
	endif
	sph_harm=sph_harm*sqrt(1.0*(2*l+1)/fourpi)

	return
	end