Overflow in Legendre polynomials

[astrojuanlu]

I am getting infinite overflows when computing Legendre polynomials of high order:

>> [data, dataD] = associated_legendre(360, 0.0, 'none');
>> data(150, 150)
ans =  -1.2553e+304
>> data(151, 151)
ans = Inf

and, as expected, NaN when computing the normalized polynomial:

>> [dataN, dataDN] = associated_legendre(360, 0.0, 'egm96');
>> dataN(150, 150)
ans = 0
>> dataN(151, 151)
ans = NaN

Does this affect the results displayed in https://people.sc.fsu.edu/~lb13f/projects/space_environment/egm96.php in some way? I haven't run the whole script yet because I am not sure whether it's going to work in Octave, but I stopped here and wanted to ask anyway.

[lukasbystricky]

Yes, this is actually because the method I used to compute the coefficients (though relatively simple) is numerically unstable for large n values. A better approach is described here: http://www.scielo.org.co/pdf/racefn/v37n145/v37n145a09.pdf

Additionally, this code has not been optimized at all so it is very slow to run anything higher than N_MAX=25. I would definitely consider vectorizing it before trying any high degree problems.

[astrojuanlu]

Thanks for your response. From the paper you linked:

The maximum possible n is 120. This limitation is due to the small value of P_{120}^{-120}, which cannot be accurately described in double precision.

So I guess the problem for representing the whole EGM-96 model still stand.

[lukasbystricky]

Sorry, you are quite right.

Here is an article dedicated to higher order polynomials for use in geopotential modeling:
http://www.dsr.inpe.br/sbsr2013/files/p0595.pdf

[astrojuanlu]

This one is interesting. Too bad the link to the source code is broken... http://www2.dem.inpe.br/hkk/software/high_geopot.htm We might have the chance to email the corresponding authors. Or implement the whole thing from scratch.

[lukasbystricky]

I won't have time to work on this in the near future (finishing up my PhD) but I can give you push access if you want to modify the code.

[astrojuanlu]

I appreciate your offer but in fact I'm a Python guy, so I plan to implement this as part of scipy/scipy#7778.

[lukasbystricky]

If you're interested, this link to the referenced C and Fortran code seems to work:
http://www.dem.inpe.br/~hkk/software/high_geopot.htm

I haven't tested the code.

[astrojuanlu]

Awesome, thanks!