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