Possible bug in Legendre Transform at fixed m

[Former user]

Hello, I am just starting out learning your library and really looking forward to using it. In trying to learn the calling sequences, I coded up a simple example based on the Fortran 2003 example from the webpage.

I simply want to do a Legendre transform at fixed m, including m values that are nonzero. I start with a pure l=2 m=1 function in physical space, then call the "spat_to_sh_ml" routine. The resulting spectral coefficients seem to be wrong, with power being placed in the l=1 m=0 mode and nowhere else. I have attached my short example Fortran code.

The routine works to machine epsilon if the physical space function is m=0, but seems to fail when m is nonzero. Is this an issue with SHTns or am I just using it incorrectly?

Thank you!

[Nathanaël Schaeffer]

Hello, this is actually the intended behavior for these “_ml” functions: the input or output arrays contain only the given m.
That is, the spatial data is only of size shtns%nlat (as you have setup in your program), and the spectral array is of size shtns%lmax+1-m, containing only the coefficients for the given m.

If you want to put the spherical harmonic coefficient at the right place of an array containing all m, you could do something like this (replacing lines 59 to 61 in your program):

    lmstart = shtns_lmidx(shtns_c,m,m)
    lmstop = shtns_lmidx(shtns_c,lmax,m)
    call spat_to_sh_ml(shtns_c, m, phys, spec(lmstart:lmstop), lmax)

Another mistake is found line 67, where you should put instead:

    do i=lmstart,lmstop

I acknowledge the documentation is totally lacking to describe this, sorry !

Note that the functions shtns_lmidx, shtns_lm2l, and shtns_lm2m all deal with indices for a full spectral array (containing all m requested by shtns_create. Instead, within a “subarray” of given m, the index (in Fortran convection) is simply l-m+1
I will try to fix the documentation, hopefully soon !

[Ryan Orvedahl]

Fantastic, thank you so much. I knew it had to be an issue on my end.

For the documentation, when ever you get around to it, I would suggest simply including an example of transforming something easy, like l=2, m=1 in the same Fortran 2003 example program that you already have. Based on that existing example, it probably only requires a few lines of new code and a useful comment. That way you can avoid fully writing up new documentation.