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 m. */
42 private static double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
43
44 /** The first eccentricity squared. */
45 private static double epssq = 0.0066943799901413169961;
46
47 /** Mean radius of IAU-66 ellipsoid, in m. */
48 private static double ellipsoidRadius = 6371200.0;
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 radians
206 * @param longitude the WGS84 longitude in radians
207 * @param height the height above the WGS84 ellipsoid in meters
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.html#GeodeticPoint">GeodeticPoint gp = new GeodeticPoint(latitude, longitude, height);
215
216 final SphericalCoordinates sph = transformToSpherical(gp);
217 final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
218 final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));
219
220 // sum up the magnetic field vector components
221 final Vector3D magFieldSph = summation(sph, vars, legendre);
222 // rotate the field to geodetic coordinates
223 final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
224 // return the magnetic elements
225 return new GeoMagneticElements(magFieldGeo);
226 }
227
228 /** Time transform the model coefficients from the base year of the model
229 * using secular variation coefficients.
230 * @param year the year to which the model shall be transformed
231 * @return a time-transformed magnetic field model
232 */
233 public GeoMagneticField transformModel(final double year) {
234
235 if (!supportsTimeTransform()) {
236 throw new OrekitException(OrekitMessages.UNSUPPORTED_TIME_TRANSFORM, modelName, String.valueOf(epoch));
237 }
238
239 // the model can only be transformed within its validity period
240 if (year < validityStart || year > validityEnd) {
241 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
242 modelName, String.valueOf(epoch), year, validityStart, validityEnd);
243 }
244
245 final double dt = year - epoch;
246 final int maxSecIndex = maxNSec * (maxNSec + 1) / 2 + maxNSec;
247
248 final GeoMagneticFieldml#GeoMagneticField">GeoMagneticField transformed = new GeoMagneticField(modelName, year, maxN, maxNSec,
249 validityStart, validityEnd);
250
251 for (int n = 1; n <= maxN; n++) {
252 for (int m = 0; m <= n; m++) {
253 final int index = n * (n + 1) / 2 + m;
254 if (index <= maxSecIndex) {
255 transformed.h[index] = h[index] + dt * dh[index];
256 transformed.g[index] = g[index] + dt * dg[index];
257 // we need a copy of the secular var coef to calculate secular change
258 transformed.dh[index] = dh[index];
259 transformed.dg[index] = dg[index];
260 } else {
261 // just copy the parts that do not have corresponding secular variation coefficients
262 transformed.h[index] = h[index];
263 transformed.g[index] = g[index];
264 }
265 }
266 }
267
268 return transformed;
269 }
270
271 /** Time transform the model coefficients from the base year of the model
272 * using a linear interpolation with a second model. The second model is
273 * required to have an adjacent validity period.
274 * @param otherModel the other magnetic field model
275 * @param year the year to which the model shall be transformed
276 * @return a time-transformed magnetic field model
277 */
278 public GeoMagneticFieldGeoMagneticField">GeoMagneticField transformModel(final GeoMagneticField otherModel, final double year) {
279
280 // the model can only be transformed within its validity period
281 if (year < validityStart || year > validityEnd) {
282 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
283 modelName, String.valueOf(epoch), year, validityStart, validityEnd);
284 }
285
286 final double factor = (year - epoch) / (otherModel.epoch - epoch);
287 final int maxNCommon = FastMath.min(maxN, otherModel.maxN);
288 final int maxNCommonIndex = maxNCommon * (maxNCommon + 1) / 2 + maxNCommon;
289
290 final int newMaxN = FastMath.max(maxN, otherModel.maxN);
291
292 final GeoMagneticFieldml#GeoMagneticField">GeoMagneticField transformed = new GeoMagneticField(modelName, year, newMaxN, 0,
293 validityStart, validityEnd);
294
295 for (int n = 1; n <= newMaxN; n++) {
296 for (int m = 0; m <= n; m++) {
297 final int index = n * (n + 1) / 2 + m;
298 if (index <= maxNCommonIndex) {
299 transformed.h[index] = h[index] + factor * (otherModel.h[index] - h[index]);
300 transformed.g[index] = g[index] + factor * (otherModel.g[index] - g[index]);
301 } else {
302 if (maxN < otherModel.maxN) {
303 transformed.h[index] = factor * otherModel.h[index];
304 transformed.g[index] = factor * otherModel.g[index];
305 } else {
306 transformed.h[index] = h[index] + factor * -h[index];
307 transformed.g[index] = g[index] + factor * -g[index];
308 }
309 }
310 }
311 }
312
313 return transformed;
314 }
315
316 /** Utility function to get a decimal year for a given day.
317 * @param day the day (1-31)
318 * @param month the month (1-12)
319 * @param year the year
320 * @return the decimal year represented by the given day
321 */
322 public static double getDecimalYear(final int day, final int month, final int year) {
323 final int[] days = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
324 final int leapYear = (((year % 4) == 0) && (((year % 100) != 0) || ((year % 400) == 0))) ? 1 : 0;
325
326 final double dayInYear = days[month - 1] + (day - 1) + (month > 2 ? leapYear : 0);
327 return (double) year + (dayInYear / (365.0d + leapYear));
328 }
329
330 /** Transform geodetic coordinates to spherical coordinates.
331 * @param gp the geodetic point
332 * @return the spherical coordinates wrt to the reference ellipsoid of the model
333 */
334 private SphericalCoordinates transformToSpherical(final GeodeticPoint gp) {
335
336 // Convert geodetic coordinates (defined by the WGS-84 reference ellipsoid)
337 // to Earth Centered Earth Fixed Cartesian coordinates, and then to spherical coordinates.
338
339 final double lat = gp.getLatitude();
340 final double heightAboveEllipsoid = gp.getAltitude();
341 final double sinLat = FastMath.sin(lat);
342
343 // compute the local radius of curvature on the reference ellipsoid
344 final double rc = a / FastMath.sqrt(1.0d - epssq * sinLat * sinLat);
345
346 // compute ECEF Cartesian coordinates of specified point (for longitude=0)
347 final double xp = (rc + heightAboveEllipsoid) * FastMath.cos(lat);
348 final double zp = (rc * (1.0d - epssq) + heightAboveEllipsoid) * sinLat;
349
350 // compute spherical radius and angle lambda and phi of specified point
351 final double r = FastMath.hypot(xp, zp);
352 return new SphericalCoordinates(r, gp.getLongitude(), FastMath.asin(zp / r));
353 }
354
355 /** Rotate the magnetic vectors to geodetic coordinates.
356 * @param sph the spherical coordinates
357 * @param gp the geodetic point
358 * @param field the magnetic field in spherical coordinates
359 * @return the magnetic field in geodetic coordinates
360 */
361 private Vector3D rotateMagneticVector(final SphericalCoordinates sph,
362 final GeodeticPoint gp,
363 final Vector3D field) {
364
365 // difference between the spherical and geodetic latitudes
366 final double psi = sph.phi - gp.getLatitude();
367
368 // rotate spherical field components to the geodetic system
369 final double Bz = field.getX() * FastMath.sin(psi) + field.getZ() * FastMath.cos(psi);
370 final double Bx = field.getX() * FastMath.cos(psi) - field.getZ() * FastMath.sin(psi);
371 final double By = field.getY();
372
373 return new Vector3D(Bx, By, Bz);
374 }
375
376 /** Computes Geomagnetic Field Elements X, Y and Z in spherical coordinate
377 * system using spherical harmonic summation.
378 * The vector Magnetic field is given by -grad V, where V is geomagnetic
379 * scalar potential. The gradient in spherical coordinates is given by:
380 * <pre>
381 * dV ^ 1 dV ^ 1 dV ^
382 * grad V = -- r + - -- t + -------- -- p
383 * dr r dt r sin(t) dp
384 * </pre>
385 * @param sph the spherical coordinates
386 * @param vars the spherical harmonic variables
387 * @param legendre the legendre function
388 * @return the magnetic field vector in spherical coordinates
389 */
390 private Vector3D summation(final SphericalCoordinates sph, final SphericalHarmonicVars vars,
391 final LegendreFunction legendre) {
392
393 int index;
394 double Bx = 0.0;
395 double By = 0.0;
396 double Bz = 0.0;
397
398 for (int n = 1; n <= maxN; n++) {
399 for (int m = 0; m <= n; m++) {
400 index = n * (n + 1) / 2 + m;
401
402 /**
403 * <pre>
404 * nMax (n+2) n m m m
405 * Bz = -SUM (n + 1) * (a/r) * SUM [g cos(m p) + h sin(m p)] P (sin(phi))
406 * n=1 m=0 n n n
407 * </pre>
408 * Equation 12 in the WMM Technical report. Derivative with respect to radius.
409 */
410 Bz -= vars.relativeRadiusPower[n] *
411 (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * (1d + n) * legendre.mP[index];
412
413 /**
414 * <pre>
415 * nMax (n+2) n m m m
416 * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
417 * n=1 m=0 n n n
418 * </pre>
419 * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
420 */
421 By += vars.relativeRadiusPower[n] *
422 (g[index] * vars.smLambda[m] - h[index] * vars.cmLambda[m]) * (double) m * legendre.mP[index];
423 /**
424 * <pre>
425 * nMax (n+2) n m m m
426 * Bx = - SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
427 * n=1 m=0 n n n
428 * </pre>
429 * Equation 10 in the WMM Technical report. Derivative with respect to latitude, divided by radius.
430 */
431 Bx -= vars.relativeRadiusPower[n] *
432 (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * legendre.mPDeriv[index];
433 }
434 }
435
436 final double cosPhi = FastMath.cos(sph.phi);
437 if (FastMath.abs(cosPhi) > 1.0e-10) {
438 By = By / cosPhi;
439 } else {
440 // special calculation for component - By - at geographic poles.
441 // To avoid using this function, make sure that the latitude is not
442 // exactly +/- π/2.
443 By = summationSpecial(sph, vars);
444 }
445
446 return new Vector3D(Bx, By, Bz);
447 }
448
449 /** Special calculation for the component By at geographic poles.
450 * @param sph the spherical coordinates
451 * @param vars the spherical harmonic variables
452 * @return the By component of the magnetic field
453 */
454 private double summationSpecial(final SphericalCoordinates sph, final SphericalHarmonicVars vars) {
455
456 double k;
457 final double sinPhi = FastMath.sin(sph.phi);
458 final double[] mPcupS = new double[maxN + 1];
459 mPcupS[0] = 1;
460 double By = 0.0;
461
462 for (int n = 1; n <= maxN; n++) {
463 final int index = n * (n + 1) / 2 + 1;
464 if (n == 1) {
465 mPcupS[n] = mPcupS[n - 1];
466 } else {
467 k = (double) (((n - 1) * (n - 1)) - 1) / (double) ((2 * n - 1) * (2 * n - 3));
468 mPcupS[n] = sinPhi * mPcupS[n - 1] - k * mPcupS[n - 2];
469 }
470
471 /**
472 * <pre>
473 * nMax (n+2) n m m m
474 * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
475 * n=1 m=0 n n n
476 * </pre>
477 * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
478 */
479 By += vars.relativeRadiusPower[n] *
480 (g[index] * vars.smLambda[1] - h[index] * vars.cmLambda[1]) * mPcupS[n] * schmidtQuasiNorm[index];
481 }
482
483 return By;
484 }
485
486 /** Utility class to hold spherical coordinates. */
487 private static class SphericalCoordinates {
488
489 /** the radius (m). */
490 private double r;
491
492 /** the azimuth angle (radians). */
493 private double lambda;
494
495 /** the polar angle (radians). */
496 private double phi;
497
498 /** Create a new spherical coordinate object.
499 * @param r the radius, meters
500 * @param lambda the lambda angle, radians
501 * @param phi the phi angle, radians
502 */
503 private SphericalCoordinates(final double r, final double lambda, final double phi) {
504 this.r = r;
505 this.lambda = lambda;
506 this.phi = phi;
507 }
508 }
509
510 /** Utility class to compute certain variables for magnetic field summation. */
511 private class SphericalHarmonicVars {
512
513 /** (Radius of Earth / Spherical radius r)^(n+2). */
514 private double[] relativeRadiusPower;
515
516 /** cos(m*lambda). */
517 private double[] cmLambda;
518
519 /** sin(m*lambda). */
520 private double[] smLambda;
521
522 /** Calculates the spherical harmonic variables for a given spherical coordinate.
523 * @param sph the spherical coordinate
524 */
525 private SphericalHarmonicVars(final SphericalCoordinates sph) {
526
527 relativeRadiusPower = new double[maxN + 1];
528
529 // Compute a table of (EARTH_REFERENCE_RADIUS / radius)^n for i in
530 // 0 .. maxN (this is much faster than calling FastMath.pow maxN+1 times).
531
532 final double p = ellipsoidRadius / sph.r;
533 relativeRadiusPower[0] = p * p;
534 for (int n = 1; n <= maxN; n++) {
535 relativeRadiusPower[n] = relativeRadiusPower[n - 1] * (ellipsoidRadius / sph.r);
536 }
537
538 // Compute tables of sin(lon * m) and cos(lon * m) for m = 0 .. maxN
539 // this is much faster than calling FastMath.sin and FastMath.cos maxN+1 times.
540
541 cmLambda = new double[maxN + 1];
542 smLambda = new double[maxN + 1];
543
544 cmLambda[0] = 1.0d;
545 smLambda[0] = 0.0d;
546
547 final double cosLambda = FastMath.cos(sph.lambda);
548 final double sinLambda = FastMath.sin(sph.lambda);
549 cmLambda[1] = cosLambda;
550 smLambda[1] = sinLambda;
551
552 for (int m = 2; m <= maxN; m++) {
553 cmLambda[m] = cmLambda[m - 1] * cosLambda - smLambda[m - 1] * sinLambda;
554 smLambda[m] = cmLambda[m - 1] * sinLambda + smLambda[m - 1] * cosLambda;
555 }
556 }
557 }
558
559 /** Utility class to compute a table of Schmidt-semi normalized associated Legendre functions. */
560 private class LegendreFunction {
561
562 /** the vector of all associated Legendre polynomials. */
563 private double[] mP;
564
565 /** the vector of derivatives of the Legendre polynomials wrt latitude. */
566 private double[] mPDeriv;
567
568 /** Calculate the Schmidt-semi normalized Legendre function.
569 * <p>
570 * <b>Note:</b> In geomagnetism, the derivatives of ALF are usually
571 * found with respect to the colatitudes. Here the derivatives are found
572 * with respect to the latitude. The difference is a sign reversal for
573 * the derivative of the Associated Legendre Functions.
574 * </p>
575 * @param x sinus of the spherical latitude (or cosinus of the spherical colatitude)
576 */
577 private LegendreFunction(final double x) {
578
579 final int numTerms = (maxN + 1) * (maxN + 2) / 2;
580
581 mP = new double[numTerms + 1];
582 mPDeriv = new double[numTerms + 1];
583
584 mP[0] = 1.0;
585 mPDeriv[0] = 0.0;
586
587 // sin (geocentric latitude) - sin_phi
588 final double z = FastMath.sqrt((1.0d - x) * (1.0d + x));
589
590 int index;
591 int index1;
592 int index2;
593
594 // First, compute the Gauss-normalized associated Legendre functions
595 for (int n = 1; n <= maxN; n++) {
596 for (int m = 0; m <= n; m++) {
597 index = n * (n + 1) / 2 + m;
598 if (n == m) {
599 index1 = (n - 1) * n / 2 + m - 1;
600 mP[index] = z * mP[index1];
601 mPDeriv[index] = z * mPDeriv[index1] + x * mP[index1];
602 } else if (n == 1 && m == 0) {
603 index1 = (n - 1) * n / 2 + m;
604 mP[index] = x * mP[index1];
605 mPDeriv[index] = x * mPDeriv[index1] - z * mP[index1];
606 } else if (n > 1 && n != m) {
607 index1 = (n - 2) * (n - 1) / 2 + m;
608 index2 = (n - 1) * n / 2 + m;
609 if (m > n - 2) {
610 mP[index] = x * mP[index2];
611 mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2];
612 } else {
613 final double k = (double) ((n - 1) * (n - 1) - (m * m)) /
614 (double) ((2 * n - 1) * (2 * n - 3));
615
616 mP[index] = x * mP[index2] - k * mP[index1];
617 mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2] - k * mPDeriv[index1];
618 }
619 }
620
621 }
622 }
623
624 // Converts the Gauss-normalized associated Legendre functions to the Schmidt quasi-normalized
625 // version using pre-computed relation stored in the variable schmidtQuasiNorm
626
627 for (int n = 1; n <= maxN; n++) {
628 for (int m = 0; m <= n; m++) {
629 index = n * (n + 1) / 2 + m;
630
631 mP[index] = mP[index] * schmidtQuasiNorm[index];
632 // The sign is changed since the new WMM routines use derivative with
633 // respect to latitude instead of co-latitude
634 mPDeriv[index] = -mPDeriv[index] * schmidtQuasiNorm[index];
635 }
636 }
637 }
638 }
639 }