IGRF -13 Geomagnetic field, Calculating Partial Derivatives of Legendre Polynomials

I have been trying to implement the IGRF-13 geomagnetic model on my own in C code. I have followed the equations from the website below and I am getting correct results for the B$$_{\phi}$$ and Br magnetic field strength components.

But right now I am completely stuck on the B$$_{\theta}$$ component which requires calculating the Partial derivatives of the Schmidt Normalized Associated Legendre Polynomials with respect to theta according to equation 3b.

Now, how do I do that? The recursive equations 19 a, b, and c are not producing the correct results, and there is no information on the internet about this.

I would really appreciate it if someone could provide a sample code that calculates the "Partial derivatives of the Associated Legendre Polynomials" or maybe provide the equations that will allow me to calculate them.

• Might be better off if you migrate this question to Scicomp SE. I have seen quite a few Legendre polynomial question over there Jun 10 at 7:48
• Does playing around with this on Wolfram Alpha help? This is a table of the first number of derivatives you're looking for (I think!), perhaps you can reconstruct the recurrence relationships from there? wolframalpha.com/input/…
– Erik
Jun 11 at 6:52
• And if the question is instead about (3a)--(3c), you may want to have a look at the "Wertz" reference to check for any typos in the notes you have.
– Erik
Jun 11 at 6:58
• If you want to migrate this to Scicomp, please flag for moderator attention and one of us will migrate it.
– gerrit
Jun 11 at 9:00

I don't have a specific reference that defines the recursion in the right form, but I've transcribed it from code I have that I've verified the output of. I think you want these, but please correct me if this doesn't give the result you expect:

Where P(n,m) is the Schmidt normalized associated Legendre polynomial of degree n and order m, and dP(n,m) is the derivative with respect to co-latitude theta.

dP(0, 0) = 0
dP(1, 1) = cos(theta)
dP(n, n) = sqrt(1 - 1/(2n))*(sin(theta)*dP(n-1, n-1) + cos(theta)*P(n-1, n-1))
dP(n, m) = (2n - 1)/sqrt(n^2 - m^2)*(cos(theta)*dP(n-1, m) - sin(theta)*P(n-1, m)) - sqrt(((n-1)^2 - m^2)/(n^2 - m^2))*dP(n-2, m)


I'm sure there is an IGRF implementation in C somewhere, but you can find C code for a different model that contains the required Legendre functionality here.

Also check the implementation of the Legendre polynomials you are using and whether they include the Condon-Shortley phase factor of (-1)^m, if they do, the normalisations also need to include it, so that it ultimately cancels back out. It is standard in geomagnetism to not apply this factor.