IGRF calculation

[leouieda]

Description of the desired feature:

I'd like to implement calculation of the IGRF field. This would require the following components:

1. A way to calculate the associated Legendre functions $P_n^m$ without the Condon-Shortly phase and with Schmidt quasi-normalization. Equations for this can be found in https://www.gnu.org/software/gsl/tr/tr001.pdf
2. Reading the model coefficients from a file.
3. Interpolating the model coefficients to the specified date. Not hard to do with Python's datetime module.

Since the model is only degree 13, we don't have to worry too much about the Legendre functions blowing up at larger orders. It would be nice to have the Holmes and Featherstone (2002) method implemented later but it's not required for this job.

Design-wise, I have a few questions:

1. Should this be a function or a class? An igrf function would be very straight forward: calculate the 3-components at the given time and location. The location can be a float or array. But this won't produce an xarray grid. A class could have IGRF.predict and IGRF.grid like our equivalent-sources. But it won't have a fit since the model is already fit. In that case, IGRF could be an instance of a SphericalHarmonics class but I'm not sure I like that.
2. The calculations are done in geocentric spherical coordinates. But usually we'll have geodetic coordinates as inputs and the vector components output should be rotated to a local geodetic frame as well. At least, this is what web calculators and the NOAA and BGS programs calculate. Doing the conversion requires knowing the ellipsoid and so Boule comes in. So: should the function/class require an bould.Ellipsoid to be passed so it can do the conversions? Should it be optional (with ellipsoid=None assuming values in geocentric coordinates)? Or do we do everything in spherical and delegate the conversions entirely to Boule outside the IGRF function (this would require adding a vector rotation function to boule.Ellipsoid)?
3. Do we package the coefficient file or download it from NOAA? I'd be hesitant of NOAA changing their download location or simply changing the file when a new IGRF comes out. So maybe packaging it would be best.

I might break this down into multiple PRs adding different part as possibly private functions (Legendre and coefficient reading ones).

Are you willing to help implement and maintain this feature?

Yes!