1 /* Copyright 2011-2012 Space Applications Services
2 * Licensed to CS Communication & Systèmes (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.models.earth;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.orekit.bodies.GeodeticPoint;
22 import org.orekit.errors.OrekitException;
23 import org.orekit.errors.OrekitMessages;
24 import org.orekit.utils.Constants;
25
26 /** Used to calculate the geomagnetic field at a given geodetic point on earth.
27 * The calculation is estimated using spherical harmonic expansion of the
28 * geomagnetic potential with coefficients provided by an actual geomagnetic
29 * field model (e.g. IGRF, WMM).
30 * <p>
31 * Based on original software written by Manoj Nair from the National
32 * Geophysical Data Center, NOAA, as part of the WMM 2010 software release
33 * (WMM_SubLibrary.c)
34 * </p>
35 * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml">World Magnetic Model Overview</a>
36 * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/soft.shtml">WMM Software Downloads</a>
37 * @author Thomas Neidhart
38 */
39 public class GeoMagneticField {
40
41 /** Semi major-axis of WGS-84 ellipsoid in km. */
42 private static double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS / 1000d;
43
44 /** The first eccentricity squared. */
45 private static double epssq = 0.0066943799901413169961;
46
47 /** Mean radius of IAU-66 ellipsoid, in km. */
48 private static double ellipsoidRadius = 6371.2;
49
50 /** The model name. */
51 private String modelName;
52
53 /** Base time of magnetic field model epoch (yrs). */
54 private double epoch;
55
56 /** C - Gauss coefficients of main geomagnetic model (nT). */
57 private double[] g;
58
59 /** C - Gauss coefficients of main geomagnetic model (nT). */
60 private double[] h;
61
62 /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
63 private double[] dg;
64
65 /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
66 private double[] dh;
67
68 /** maximum degree of spherical harmonic model. */
69 private int maxN;
70
71 /** maximum degree of spherical harmonic secular variations. */
72 private int maxNSec;
73
74 /** the validity start of this magnetic field model. */
75 private double validityStart;
76 /** the validity end of this magnetic field model. */
77 private double validityEnd;
78
79 /** Pre-calculated ratio between gauss-normalized and schmidt quasi-normalized
80 * associated Legendre functions.
81 */
82 private double[] schmidtQuasiNorm;
83
84 /** Create a new geomagnetic field model with the given parameters. Internal
85 * structures are initialized according to the specified degrees of the main
86 * and secular variations.
87 * @param modelName the model name
88 * @param epoch the epoch of the model
89 * @param maxN the maximum degree of the main model
90 * @param maxNSec the maximum degree of the secular variations
91 * @param validityStart validity start of this model
92 * @param validityEnd validity end of this model
93 */
94 protected GeoMagneticField(final String modelName, final double epoch,
95 final int maxN, final int maxNSec,
96 final double validityStart, final double validityEnd) {
97
98 this.modelName = modelName;
99 this.epoch = epoch;
100 this.maxN = maxN;
101 this.maxNSec = maxNSec;
102
103 this.validityStart = validityStart;
104 this.validityEnd = validityEnd;
105
106 // initialize main and secular field coefficient arrays
107 final int maxMainFieldTerms = (maxN + 1) * (maxN + 2) / 2;
108 g = new double[maxMainFieldTerms];
109 h = new double[maxMainFieldTerms];
110
111 final int maxSecularFieldTerms = (maxNSec + 1) * (maxNSec + 2) / 2;
112 dg = new double[maxSecularFieldTerms];
113 dh = new double[maxSecularFieldTerms];
114
115 // pre-calculate the ratio between gauss-normalized and schmidt quasi-normalized
116 // associated Legendre functions as they depend only on the degree of the model.
117
118 schmidtQuasiNorm = new double[maxMainFieldTerms + 1];
119 schmidtQuasiNorm[0] = 1.0;
120
121 int index;
122 int index1;
123 for (int n = 1; n <= maxN; n++) {
124 index = n * (n + 1) / 2;
125 index1 = (n - 1) * n / 2;
126
127 // for m = 0
128 schmidtQuasiNorm[index] =
129 schmidtQuasiNorm[index1] * (double) (2 * n - 1) / (double) n;
130
131 for (int m = 1; m <= n; m++) {
132 index = n * (n + 1) / 2 + m;
133 index1 = n * (n + 1) / 2 + m - 1;
134 schmidtQuasiNorm[index] =
135 schmidtQuasiNorm[index1] *
136 FastMath.sqrt((double) ((n - m + 1) * (m == 1 ? 2 : 1)) / (double) (n + m));
137 }
138 }
139 }
140
141 /** Returns the epoch for this magnetic field model.
142 * @return the epoch
143 */
144 public double getEpoch() {
145 return epoch;
146 }
147
148 /** Returns the model name.
149 * @return the model name
150 */
151 public String getModelName() {
152 return modelName;
153 }
154
155 /** Returns the start of the validity period for this model.
156 * @return the validity start as decimal year
157 */
158 public double validFrom() {
159 return validityStart;
160 }
161
162 /** Returns the end of the validity period for this model.
163 * @return the validity end as decimal year
164 */
165 public double validTo() {
166 return validityEnd;
167 }
168
169 /** Indicates whether this model supports time transformation or not.
170 * @return <code>true</code> if this model can be transformed within its
171 * validity period, <code>false</code> otherwise
172 */
173 public boolean supportsTimeTransform() {
174 return maxNSec > 0;
175 }
176
177 /** Set the given main field coefficients.
178 * @param n the n index
179 * @param m the m index
180 * @param gnm the g coefficient at position n,m
181 * @param hnm the h coefficient at position n,m
182 */
183 protected void setMainFieldCoefficients(final int n, final int m,
184 final double gnm, final double hnm) {
185 final int index = n * (n + 1) / 2 + m;
186 g[index] = gnm;
187 h[index] = hnm;
188 }
189
190 /** Set the given secular variation coefficients.
191 * @param n the n index
192 * @param m the m index
193 * @param dgnm the dg coefficient at position n,m
194 * @param dhnm the dh coefficient at position n,m
195 */
196 protected void setSecularVariationCoefficients(final int n, final int m,
197 final double dgnm, final double dhnm) {
198 final int index = n * (n + 1) / 2 + m;
199 dg[index] = dgnm;
200 dh[index] = dhnm;
201 }
202
203 /** Calculate the magnetic field at the specified geodetic point identified
204 * by latitude, longitude and altitude.
205 * @param latitude the WGS84 latitude in decimal degrees
206 * @param longitude the WGS84 longitude in decimal degrees
207 * @param height the height above the WGS84 ellipsoid in kilometers
208 * @return the {@link GeoMagneticElements} at the given geodetic point
209 */
210 public GeoMagneticElements calculateField(final double latitude,
211 final double longitude,
212 final double height) {
213
214 final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(latitude),
215 FastMath.toRadians(longitude),
216 height * 1000d);
217
218 final SphericalCoordinates sph = transformToSpherical(gp);
219 final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
220 final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));
221
222 // sum up the magnetic field vector components
223 final Vector3D magFieldSph = summation(sph, vars, legendre);
224 // rotate the field to geodetic coordinates
225 final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
226 // return the magnetic elements
227 return new GeoMagneticElements(magFieldGeo);
228 }
229
230 /** Time transform the model coefficients from the base year of the model
231 * using secular variation coefficients.
232 * @param year the year to which the model shall be transformed
233 * @return a time-transformed magnetic field model
234 * @throws OrekitException if the specified year is outside the validity period of the
235 * model or the model does not support time transformations
236 * (i.e. no secular variations available)
237 */
238 public GeoMagneticField transformModel(final double year) throws OrekitException {
239
240 if (!supportsTimeTransform()) {
241 throw new OrekitException(OrekitMessages.UNSUPPORTED_TIME_TRANSFORM, modelName, String.valueOf(epoch));
242 }
243
244 // the model can only be transformed within its validity period
245 if (year < validityStart || year > validityEnd) {
246 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
247 modelName, String.valueOf(epoch), year, validityStart, validityEnd);
248 }
249
250 final double dt = year - epoch;
251 final int maxSecIndex = maxNSec * (maxNSec + 1) / 2 + maxNSec;
252
253 final GeoMagneticField transformed = new GeoMagneticField(modelName, year, maxN, maxNSec,
254 validityStart, validityEnd);
255
256 for (int n = 1; n <= maxN; n++) {
257 for (int m = 0; m <= n; m++) {
258 final int index = n * (n + 1) / 2 + m;
259 if (index <= maxSecIndex) {
260 transformed.h[index] = h[index] + dt * dh[index];
261 transformed.g[index] = g[index] + dt * dg[index];
262 // we need a copy of the secular var coef to calculate secular change
263 transformed.dh[index] = dh[index];
264 transformed.dg[index] = dg[index];
265 } else {
266 // just copy the parts that do not have corresponding secular variation coefficients
267 transformed.h[index] = h[index];
268 transformed.g[index] = g[index];
269 }
270 }
271 }
272
273 return transformed;
274 }
275
276 /** Time transform the model coefficients from the base year of the model
277 * using a linear interpolation with a second model. The second model is
278 * required to have an adjacent validity period.
279 * @param otherModel the other magnetic field model
280 * @param year the year to which the model shall be transformed
281 * @return a time-transformed magnetic field model
282 * @throws OrekitException if the specified year is outside the validity period of the
283 * model or the model does not support time transformations
284 * (i.e. no secular variations available)
285 */
286 public GeoMagneticField transformModel(final GeoMagneticField otherModel, final double year)
287 throws OrekitException {
288
289 // the model can only be transformed within its validity period
290 if (year < validityStart || year > validityEnd) {
291 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
292 modelName, String.valueOf(epoch), year, validityStart, validityEnd);
293 }
294
295 final double factor = (year - epoch) / (otherModel.epoch - epoch);
296 final int maxNCommon = FastMath.min(maxN, otherModel.maxN);
297 final int maxNCommonIndex = maxNCommon * (maxNCommon + 1) / 2 + maxNCommon;
298
299 final int newMaxN = FastMath.max(maxN, otherModel.maxN);
300
301 final GeoMagneticField transformed = new GeoMagneticField(modelName, year, newMaxN, 0,
302 validityStart, validityEnd);
303
304 for (int n = 1; n <= newMaxN; n++) {
305 for (int m = 0; m <= n; m++) {
306 final int index = n * (n + 1) / 2 + m;
307 if (index <= maxNCommonIndex) {
308 transformed.h[index] = h[index] + factor * (otherModel.h[index] - h[index]);
309 transformed.g[index] = g[index] + factor * (otherModel.g[index] - g[index]);
310 } else {
311 if (maxN < otherModel.maxN) {
312 transformed.h[index] = factor * otherModel.h[index];
313 transformed.g[index] = factor * otherModel.g[index];
314 } else {
315 transformed.h[index] = h[index] + factor * -h[index];
316 transformed.g[index] = g[index] + factor * -g[index];
317 }
318 }
319 }
320 }
321
322 return transformed;
323 }
324
325 /** Utility function to get a decimal year for a given day.
326 * @param day the day (1-31)
327 * @param month the month (1-12)
328 * @param year the year
329 * @return the decimal year represented by the given day
330 */
331 public static double getDecimalYear(final int day, final int month, final int year) {
332 final int[] days = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
333 final int leapYear = (((year % 4) == 0) && (((year % 100) != 0) || ((year % 400) == 0))) ? 1 : 0;
334
335 final double dayInYear = days[month - 1] + (day - 1) + (month > 2 ? leapYear : 0);
336 return (double) year + (dayInYear / (365.0d + leapYear));
337 }
338
339 /** Transform geodetic coordinates to spherical coordinates.
340 * @param gp the geodetic point
341 * @return the spherical coordinates wrt to the reference ellipsoid of the model
342 */
343 private SphericalCoordinates transformToSpherical(final GeodeticPoint gp) {
344
345 // Convert geodetic coordinates (defined by the WGS-84 reference ellipsoid)
346 // to Earth Centered Earth Fixed Cartesian coordinates, and then to spherical coordinates.
347
348 final double lat = gp.getLatitude();
349 final double heightAboveEllipsoid = gp.getAltitude() / 1000d;
350 final double sinLat = FastMath.sin(lat);
351
352 // compute the local radius of curvature on the reference ellipsoid
353 final double rc = a / FastMath.sqrt(1.0d - epssq * sinLat * sinLat);
354
355 // compute ECEF Cartesian coordinates of specified point (for longitude=0)
356 final double xp = (rc + heightAboveEllipsoid) * FastMath.cos(lat);
357 final double zp = (rc * (1.0d - epssq) + heightAboveEllipsoid) * sinLat;
358
359 // compute spherical radius and angle lambda and phi of specified point
360 final double r = FastMath.hypot(xp, zp);
361 return new SphericalCoordinates(r, gp.getLongitude(), FastMath.asin(zp / r));
362 }
363
364 /** Rotate the magnetic vectors to geodetic coordinates.
365 * @param sph the spherical coordinates
366 * @param gp the geodetic point
367 * @param field the magnetic field in spherical coordinates
368 * @return the magnetic field in geodetic coordinates
369 */
370 private Vector3D rotateMagneticVector(final SphericalCoordinates sph,
371 final GeodeticPoint gp,
372 final Vector3D field) {
373
374 // difference between the spherical and geodetic latitudes
375 final double psi = sph.phi - gp.getLatitude();
376
377 // rotate spherical field components to the geodetic system
378 final double Bz = field.getX() * FastMath.sin(psi) + field.getZ() * FastMath.cos(psi);
379 final double Bx = field.getX() * FastMath.cos(psi) - field.getZ() * FastMath.sin(psi);
380 final double By = field.getY();
381
382 return new Vector3D(Bx, By, Bz);
383 }
384
385 /** Computes Geomagnetic Field Elements X, Y and Z in spherical coordinate
386 * system using spherical harmonic summation.
387 * The vector Magnetic field is given by -grad V, where V is geomagnetic
388 * scalar potential. The gradient in spherical coordinates is given by:
389 * <pre>
390 * dV ^ 1 dV ^ 1 dV ^
391 * grad V = -- r + - -- t + -------- -- p
392 * dr r dt r sin(t) dp
393 * </pre>
394 * @param sph the spherical coordinates
395 * @param vars the spherical harmonic variables
396 * @param legendre the legendre function
397 * @return the magnetic field vector in spherical coordinates
398 */
399 private Vector3D summation(final SphericalCoordinates sph, final SphericalHarmonicVars vars,
400 final LegendreFunction legendre) {
401
402 int index;
403 double Bx = 0.0;
404 double By = 0.0;
405 double Bz = 0.0;
406
407 for (int n = 1; n <= maxN; n++) {
408 for (int m = 0; m <= n; m++) {
409 index = n * (n + 1) / 2 + m;
410
411 /**
412 * <pre>
413 * nMax (n+2) n m m m
414 * Bz = -SUM (n + 1) * (a/r) * SUM [g cos(m p) + h sin(m p)] P (sin(phi))
415 * n=1 m=0 n n n
416 * </pre>
417 * Equation 12 in the WMM Technical report. Derivative with respect to radius.
418 */
419 Bz -= vars.relativeRadiusPower[n] *
420 (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * (1d + n) * legendre.mP[index];
421
422 /**
423 * <pre>
424 * nMax (n+2) n m m m
425 * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
426 * n=1 m=0 n n n
427 * </pre>
428 * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
429 */
430 By += vars.relativeRadiusPower[n] *
431 (g[index] * vars.smLambda[m] - h[index] * vars.cmLambda[m]) * (double) m * legendre.mP[index];
432 /**
433 * <pre>
434 * nMax (n+2) n m m m
435 * Bx = - SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
436 * n=1 m=0 n n n
437 * </pre>
438 * Equation 10 in the WMM Technical report. Derivative with respect to latitude, divided by radius.
439 */
440 Bx -= vars.relativeRadiusPower[n] *
441 (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * legendre.mPDeriv[index];
442 }
443 }
444
445 final double cosPhi = FastMath.cos(sph.phi);
446 if (FastMath.abs(cosPhi) > 1.0e-10) {
447 By = By / cosPhi;
448 } else {
449 // special calculation for component - By - at geographic poles.
450 // To avoid using this function, make sure that the latitude is not
451 // exactly +/-90.
452 By = summationSpecial(sph, vars);
453 }
454
455 return new Vector3D(Bx, By, Bz);
456 }
457
458 /** Special calculation for the component By at geographic poles.
459 * @param sph the spherical coordinates
460 * @param vars the spherical harmonic variables
461 * @return the By component of the magnetic field
462 */
463 private double summationSpecial(final SphericalCoordinates sph, final SphericalHarmonicVars vars) {
464
465 double k;
466 final double sinPhi = FastMath.sin(sph.phi);
467 final double[] mPcupS = new double[maxN + 1];
468 mPcupS[0] = 1;
469 double By = 0.0;
470
471 for (int n = 1; n <= maxN; n++) {
472 final int index = n * (n + 1) / 2 + 1;
473 if (n == 1) {
474 mPcupS[n] = mPcupS[n - 1];
475 } else {
476 k = (double) (((n - 1) * (n - 1)) - 1) / (double) ((2 * n - 1) * (2 * n - 3));
477 mPcupS[n] = sinPhi * mPcupS[n - 1] - k * mPcupS[n - 2];
478 }
479
480 /**
481 * <pre>
482 * nMax (n+2) n m m m
483 * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
484 * n=1 m=0 n n n
485 * </pre>
486 * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
487 */
488 By += vars.relativeRadiusPower[n] *
489 (g[index] * vars.smLambda[1] - h[index] * vars.cmLambda[1]) * mPcupS[n] * schmidtQuasiNorm[index];
490 }
491
492 return By;
493 }
494
495 /** Utility class to hold spherical coordinates. */
496 private static class SphericalCoordinates {
497
498 /** the radius. */
499 private double r;
500
501 /** the azimuth angle. */
502 private double lambda;
503
504 /** the polar angle. */
505 private double phi;
506
507 /** Create a new spherical coordinate object.
508 * @param r the radius
509 * @param lambda the lambda angle
510 * @param phi the phi angle
511 */
512 private SphericalCoordinates(final double r, final double lambda, final double phi) {
513 this.r = r;
514 this.lambda = lambda;
515 this.phi = phi;
516 }
517 }
518
519 /** Utility class to compute certain variables for magnetic field summation. */
520 private class SphericalHarmonicVars {
521
522 /** (Radius of Earth / Spherical radius r)^(n+2). */
523 private double[] relativeRadiusPower;
524
525 /** cos(m*lambda). */
526 private double[] cmLambda;
527
528 /** sin(m*lambda). */
529 private double[] smLambda;
530
531 /** Calculates the spherical harmonic variables for a given spherical coordinate.
532 * @param sph the spherical coordinate
533 */
534 private SphericalHarmonicVars(final SphericalCoordinates sph) {
535
536 relativeRadiusPower = new double[maxN + 1];
537
538 // Compute a table of (EARTH_REFERENCE_RADIUS_KM / radius)^n for i in
539 // 0 .. maxN (this is much faster than calling FastMath.pow maxN+1 times).
540
541 final double p = ellipsoidRadius / sph.r;
542 relativeRadiusPower[0] = p * p;
543 for (int n = 1; n <= maxN; n++) {
544 relativeRadiusPower[n] = relativeRadiusPower[n - 1] * (ellipsoidRadius / sph.r);
545 }
546
547 // Compute tables of sin(lon * m) and cos(lon * m) for m = 0 .. maxN
548 // this is much faster than calling FastMath.sin and FastMath.cos maxN+1 times.
549
550 cmLambda = new double[maxN + 1];
551 smLambda = new double[maxN + 1];
552
553 cmLambda[0] = 1.0d;
554 smLambda[0] = 0.0d;
555
556 final double cosLambda = FastMath.cos(sph.lambda);
557 final double sinLambda = FastMath.sin(sph.lambda);
558 cmLambda[1] = cosLambda;
559 smLambda[1] = sinLambda;
560
561 for (int m = 2; m <= maxN; m++) {
562 cmLambda[m] = cmLambda[m - 1] * cosLambda - smLambda[m - 1] * sinLambda;
563 smLambda[m] = cmLambda[m - 1] * sinLambda + smLambda[m - 1] * cosLambda;
564 }
565 }
566 }
567
568 /** Utility class to compute a table of Schmidt-semi normalized associated Legendre functions. */
569 private class LegendreFunction {
570
571 /** the vector of all associated Legendre polynomials. */
572 private double[] mP;
573
574 /** the vector of derivatives of the Legendre polynomials wrt latitude. */
575 private double[] mPDeriv;
576
577 /** Calculate the Schmidt-semi normalized Legendre function.
578 * <p>
579 * <b>Note:</b> In geomagnetism, the derivatives of ALF are usually
580 * found with respect to the colatitudes. Here the derivatives are found
581 * with respect to the latitude. The difference is a sign reversal for
582 * the derivative of the Associated Legendre Functions.
583 * </p>
584 * @param x sinus of the spherical latitude (or cosinus of the spherical colatitude)
585 */
586 private LegendreFunction(final double x) {
587
588 final int numTerms = (maxN + 1) * (maxN + 2) / 2;
589
590 mP = new double[numTerms + 1];
591 mPDeriv = new double[numTerms + 1];
592
593 mP[0] = 1.0;
594 mPDeriv[0] = 0.0;
595
596 // sin (geocentric latitude) - sin_phi
597 final double z = FastMath.sqrt((1.0d - x) * (1.0d + x));
598
599 int index;
600 int index1;
601 int index2;
602
603 // First, compute the Gauss-normalized associated Legendre functions
604 for (int n = 1; n <= maxN; n++) {
605 for (int m = 0; m <= n; m++) {
606 index = n * (n + 1) / 2 + m;
607 if (n == m) {
608 index1 = (n - 1) * n / 2 + m - 1;
609 mP[index] = z * mP[index1];
610 mPDeriv[index] = z * mPDeriv[index1] + x * mP[index1];
611 } else if (n == 1 && m == 0) {
612 index1 = (n - 1) * n / 2 + m;
613 mP[index] = x * mP[index1];
614 mPDeriv[index] = x * mPDeriv[index1] - z * mP[index1];
615 } else if (n > 1 && n != m) {
616 index1 = (n - 2) * (n - 1) / 2 + m;
617 index2 = (n - 1) * n / 2 + m;
618 if (m > n - 2) {
619 mP[index] = x * mP[index2];
620 mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2];
621 } else {
622 final double k = (double) ((n - 1) * (n - 1) - (m * m)) /
623 (double) ((2 * n - 1) * (2 * n - 3));
624
625 mP[index] = x * mP[index2] - k * mP[index1];
626 mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2] - k * mPDeriv[index1];
627 }
628 }
629
630 }
631 }
632
633 // Converts the Gauss-normalized associated Legendre functions to the Schmidt quasi-normalized
634 // version using pre-computed relation stored in the variable schmidtQuasiNorm
635
636 for (int n = 1; n <= maxN; n++) {
637 for (int m = 0; m <= n; m++) {
638 index = n * (n + 1) / 2 + m;
639
640 mP[index] = mP[index] * schmidtQuasiNorm[index];
641 // The sign is changed since the new WMM routines use derivative with
642 // respect to latitude instead of co-latitude
643 mPDeriv[index] = -mPDeriv[index] * schmidtQuasiNorm[index];
644 }
645 }
646 }
647 }
648 }