[SciPy-user] Spherical Harmonics

[Pearu Peterson]

On Sat, 22 Jan 2005, Arne Keller wrote:

> I'm trying to use the sph_harm special function, but it has a very
> strange behavior, sometimes it fails with a segmentation fault.
>
> If someone can try to execute these lines of codes, to see if he obtains
> the same result:
>
> ##########################################################
> ## file test.py
> ##########################################################
> import scipy
> from scipy.special import sph_harm,lpmn,gammaln
> from scipy import *
>
> #n=1 #uncommenting this line produce a segmentation fault!
>
> #taken from special/basic.py:
> def sph_harmonic(m,n,theta,phi):
>    """inputs of (m,n,theta,phi) returns spherical harmonic of order
>    m,n (|m|<=n) and argument theta and phi:  Y^m_n(theta,phi)
>    """
>    x = cos(phi)
>    m,n = int(m), int(n)
>    Pmn,Pmnd = lpmn(m,n,x)
>    val = Pmn[m,n]
>    val *= sqrt((2*n+1)/4.0/pi)
>    val *= exp(0.5*(gammaln(n-m+1)-gammaln(n+m+1)))
>    val *= exp(1j*m*theta)
>    return val
>
> n=1
> m=1
> theta=pi/2.
> phi=0.
> print sph_harmonic(m,n,phi,theta)
>
> ################End test #################################
> it produces (-0.345494149471+0j)
>
>
> but When I uncomment the first line ('n=1'), then a segmentation fault
> occurs (if the line is uncommented but with 'n=2' instead of 'n=1' then
> the execution is normal).
>
> Please may anybody help?

This segfault is probably related to the following code in 
special/specfun/specfun.f starting at line #6927:

         DO 35 I=0,M
35         PM(I,I+1)=(2.0D0*I+1.0D0)*X*PM(I,I)

where PM has DIMENSION(0:M,0:N). So, if N==M then PM(M,M+1) will point out 
of the memory area located for PM.

At the moment I am not sure whether the bug is in specfun.f or there 
should be an requirement m < n, specfun.pyf uses requirement m<=n.

Pearu