GeoMagneticField.java

  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. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  19. import org.hipparchus.util.FastMath;
  20. import org.hipparchus.util.SinCos;
  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. /** Used to calculate the geomagnetic field at a given geodetic point on earth.
  26.  * The calculation is estimated using spherical harmonic expansion of the
  27.  * geomagnetic potential with coefficients provided by an actual geomagnetic
  28.  * field model (e.g. IGRF, WMM).
  29.  * <p>
  30.  * Based on original software written by Manoj Nair from the National
  31.  * Geophysical Data Center, NOAA, as part of the WMM 2010 software release
  32.  * (WMM_SubLibrary.c)
  33.  * </p>
  34.  * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml">World Magnetic Model Overview</a>
  35.  * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/soft.shtml">WMM Software Downloads</a>
  36.  * @author Thomas Neidhart
  37.  */
  38. public class GeoMagneticField {

  39.     /** Semi major-axis of WGS-84 ellipsoid in m. */
  40.     private static double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;

  41.     /** The first eccentricity squared. */
  42.     private static double epssq = 0.0066943799901413169961;

  43.     /** Mean radius of IAU-66 ellipsoid, in m. */
  44.     private static double ellipsoidRadius = 6371200.0;

  45.     /** The model name. */
  46.     private String modelName;

  47.     /** Base time of magnetic field model epoch (yrs). */
  48.     private double epoch;

  49.     /** C - Gauss coefficients of main geomagnetic model (nT). */
  50.     private double[] g;

  51.     /** C - Gauss coefficients of main geomagnetic model (nT). */
  52.     private double[] h;

  53.     /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
  54.     private double[] dg;

  55.     /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
  56.     private double[] dh;

  57.     /** maximum degree of spherical harmonic model. */
  58.     private int maxN;

  59.     /** maximum degree of spherical harmonic secular variations. */
  60.     private int maxNSec;

  61.     /** the validity start of this magnetic field model. */
  62.     private double validityStart;
  63.     /** the validity end of this magnetic field model. */
  64.     private double validityEnd;

  65.     /** Pre-calculated ratio between gauss-normalized and schmidt quasi-normalized
  66.      * associated Legendre functions.
  67.      */
  68.     private double[] schmidtQuasiNorm;

  69.     /** Create a new geomagnetic field model with the given parameters. Internal
  70.      * structures are initialized according to the specified degrees of the main
  71.      * and secular variations.
  72.      * @param modelName the model name
  73.      * @param epoch the epoch of the model
  74.      * @param maxN the maximum degree of the main model
  75.      * @param maxNSec the maximum degree of the secular variations
  76.      * @param validityStart validity start of this model
  77.      * @param validityEnd validity end of this model
  78.      */
  79.     protected GeoMagneticField(final String modelName, final double epoch,
  80.                                final int maxN, final int maxNSec,
  81.                                final double validityStart, final double validityEnd) {

  82.         this.modelName = modelName;
  83.         this.epoch = epoch;
  84.         this.maxN = maxN;
  85.         this.maxNSec = maxNSec;

  86.         this.validityStart = validityStart;
  87.         this.validityEnd = validityEnd;

  88.         // initialize main and secular field coefficient arrays
  89.         final int maxMainFieldTerms = (maxN + 1) * (maxN + 2) / 2;
  90.         g = new double[maxMainFieldTerms];
  91.         h = new double[maxMainFieldTerms];

  92.         final int maxSecularFieldTerms = (maxNSec + 1) * (maxNSec + 2) / 2;
  93.         dg = new double[maxSecularFieldTerms];
  94.         dh = new double[maxSecularFieldTerms];

  95.         // pre-calculate the ratio between gauss-normalized and schmidt quasi-normalized
  96.         // associated Legendre functions as they depend only on the degree of the model.

  97.         schmidtQuasiNorm = new double[maxMainFieldTerms + 1];
  98.         schmidtQuasiNorm[0] = 1.0;

  99.         int index;
  100.         int index1;
  101.         for (int n = 1; n <= maxN; n++) {
  102.             index = n * (n + 1) / 2;
  103.             index1 = (n - 1) * n / 2;

  104.             // for m = 0
  105.             schmidtQuasiNorm[index] =
  106.                 schmidtQuasiNorm[index1] * (double) (2 * n - 1) / (double) n;

  107.             for (int m = 1; m <= n; m++) {
  108.                 index = n * (n + 1) / 2 + m;
  109.                 index1 = n * (n + 1) / 2 + m - 1;
  110.                 schmidtQuasiNorm[index] =
  111.                     schmidtQuasiNorm[index1] *
  112.                     FastMath.sqrt((double) ((n - m + 1) * (m == 1 ? 2 : 1)) / (double) (n + m));
  113.             }
  114.         }
  115.     }

  116.     /** Returns the epoch for this magnetic field model.
  117.      * @return the epoch
  118.      */
  119.     public double getEpoch() {
  120.         return epoch;
  121.     }

  122.     /** Returns the model name.
  123.      * @return the model name
  124.      */
  125.     public String getModelName() {
  126.         return modelName;
  127.     }

  128.     /** Returns the start of the validity period for this model.
  129.      * @return the validity start as decimal year
  130.      */
  131.     public double validFrom() {
  132.         return validityStart;
  133.     }

  134.     /** Returns the end of the validity period for this model.
  135.      * @return the validity end as decimal year
  136.      */
  137.     public double validTo() {
  138.         return validityEnd;
  139.     }

  140.     /** Indicates whether this model supports time transformation or not.
  141.      * @return <code>true</code> if this model can be transformed within its
  142.      *         validity period, <code>false</code> otherwise
  143.      */
  144.     public boolean supportsTimeTransform() {
  145.         return maxNSec > 0;
  146.     }

  147.     /** Set the given main field coefficients.
  148.      * @param n the n index
  149.      * @param m the m index
  150.      * @param gnm the g coefficient at position n,m
  151.      * @param hnm the h coefficient at position n,m
  152.      */
  153.     protected void setMainFieldCoefficients(final int n, final int m,
  154.                                             final double gnm, final double hnm) {
  155.         final int index = n * (n + 1) / 2 + m;
  156.         g[index] = gnm;
  157.         h[index] = hnm;
  158.     }

  159.     /** Set the given secular variation coefficients.
  160.      * @param n the n index
  161.      * @param m the m index
  162.      * @param dgnm the dg coefficient at position n,m
  163.      * @param dhnm the dh coefficient at position n,m
  164.      */
  165.     protected void setSecularVariationCoefficients(final int n, final int m,
  166.                                                    final double dgnm, final double dhnm) {
  167.         final int index = n * (n + 1) / 2 + m;
  168.         dg[index] = dgnm;
  169.         dh[index] = dhnm;
  170.     }

  171.     /** Calculate the magnetic field at the specified geodetic point identified
  172.      * by latitude, longitude and altitude.
  173.      * @param latitude the WGS84 latitude in radians
  174.      * @param longitude the WGS84 longitude in radians
  175.      * @param height the height above the WGS84 ellipsoid in meters
  176.      * @return the {@link GeoMagneticElements} at the given geodetic point
  177.      */
  178.     public GeoMagneticElements calculateField(final double latitude,
  179.                                               final double longitude,
  180.                                               final double height) {

  181.         final GeodeticPoint gp = new GeodeticPoint(latitude, longitude, height);

  182.         final SphericalCoordinates sph = transformToSpherical(gp);
  183.         final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
  184.         final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));

  185.         // sum up the magnetic field vector components
  186.         final Vector3D magFieldSph = summation(sph, vars, legendre);
  187.         // rotate the field to geodetic coordinates
  188.         final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
  189.         // return the magnetic elements
  190.         return new GeoMagneticElements(magFieldGeo);
  191.     }

  192.     /** Time transform the model coefficients from the base year of the model
  193.      * using secular variation coefficients.
  194.      * @param year the year to which the model shall be transformed
  195.      * @return a time-transformed magnetic field model
  196.      */
  197.     public GeoMagneticField transformModel(final double year) {

  198.         if (!supportsTimeTransform()) {
  199.             throw new OrekitException(OrekitMessages.UNSUPPORTED_TIME_TRANSFORM, modelName, String.valueOf(epoch));
  200.         }

  201.         // the model can only be transformed within its validity period
  202.         if (year < validityStart || year > validityEnd) {
  203.             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
  204.                                       modelName, String.valueOf(epoch), year, validityStart, validityEnd);
  205.         }

  206.         final double dt = year - epoch;
  207.         final int maxSecIndex = maxNSec * (maxNSec + 1) / 2 + maxNSec;

  208.         final GeoMagneticField transformed = new GeoMagneticField(modelName, year, maxN, maxNSec,
  209.                                                                   validityStart, validityEnd);

  210.         for (int n = 1; n <= maxN; n++) {
  211.             for (int m = 0; m <= n; m++) {
  212.                 final int index = n * (n + 1) / 2 + m;
  213.                 if (index <= maxSecIndex) {
  214.                     transformed.h[index] = h[index] + dt * dh[index];
  215.                     transformed.g[index] = g[index] + dt * dg[index];
  216.                     // we need a copy of the secular var coef to calculate secular change
  217.                     transformed.dh[index] = dh[index];
  218.                     transformed.dg[index] = dg[index];
  219.                 } else {
  220.                     // just copy the parts that do not have corresponding secular variation coefficients
  221.                     transformed.h[index] = h[index];
  222.                     transformed.g[index] = g[index];
  223.                 }
  224.             }
  225.         }

  226.         return transformed;
  227.     }

  228.     /** Time transform the model coefficients from the base year of the model
  229.      * using a linear interpolation with a second model. The second model is
  230.      * required to have an adjacent validity period.
  231.      * @param otherModel the other magnetic field model
  232.      * @param year the year to which the model shall be transformed
  233.      * @return a time-transformed magnetic field model
  234.      */
  235.     public GeoMagneticField transformModel(final GeoMagneticField otherModel, final double year) {

  236.         // the model can only be transformed within its validity period
  237.         if (year < validityStart || year > validityEnd) {
  238.             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
  239.                                       modelName, String.valueOf(epoch), year, validityStart, validityEnd);
  240.         }

  241.         final double factor = (year - epoch) / (otherModel.epoch - epoch);
  242.         final int maxNCommon = FastMath.min(maxN, otherModel.maxN);
  243.         final int maxNCommonIndex = maxNCommon * (maxNCommon + 1) / 2 + maxNCommon;

  244.         final int newMaxN = FastMath.max(maxN, otherModel.maxN);

  245.         final GeoMagneticField transformed = new GeoMagneticField(modelName, year, newMaxN, 0,
  246.                                                                   validityStart, validityEnd);

  247.         for (int n = 1; n <= newMaxN; n++) {
  248.             for (int m = 0; m <= n; m++) {
  249.                 final int index = n * (n + 1) / 2 + m;
  250.                 if (index <= maxNCommonIndex) {
  251.                     transformed.h[index] = h[index] + factor * (otherModel.h[index] - h[index]);
  252.                     transformed.g[index] = g[index] + factor * (otherModel.g[index] - g[index]);
  253.                 } else {
  254.                     if (maxN < otherModel.maxN) {
  255.                         transformed.h[index] = factor * otherModel.h[index];
  256.                         transformed.g[index] = factor * otherModel.g[index];
  257.                     } else {
  258.                         transformed.h[index] = h[index] + factor * -h[index];
  259.                         transformed.g[index] = g[index] + factor * -g[index];
  260.                     }
  261.                 }
  262.             }
  263.         }

  264.         return transformed;
  265.     }

  266.     /** Utility function to get a decimal year for a given day.
  267.      * @param day the day (1-31)
  268.      * @param month the month (1-12)
  269.      * @param year the year
  270.      * @return the decimal year represented by the given day
  271.      */
  272.     public static double getDecimalYear(final int day, final int month, final int year) {
  273.         final int[] days = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
  274.         final int leapYear = ((year % 4) == 0 && ((year % 100) != 0 || (year % 400) == 0)) ? 1 : 0;

  275.         final double dayInYear = days[month - 1] + (day - 1) + (month > 2 ? leapYear : 0);
  276.         return (double) year + (dayInYear / (365.0d + leapYear));
  277.     }

  278.     /** Transform geodetic coordinates to spherical coordinates.
  279.      * @param gp the geodetic point
  280.      * @return the spherical coordinates wrt to the reference ellipsoid of the model
  281.      */
  282.     private SphericalCoordinates transformToSpherical(final GeodeticPoint gp) {

  283.         // Convert geodetic coordinates (defined by the WGS-84 reference ellipsoid)
  284.         // to Earth Centered Earth Fixed Cartesian coordinates, and then to spherical coordinates.

  285.         final double lat = gp.getLatitude();
  286.         final SinCos sc  = FastMath.sinCos(lat);
  287.         final double heightAboveEllipsoid = gp.getAltitude();

  288.         // compute the local radius of curvature on the reference ellipsoid
  289.         final double rc = a / FastMath.sqrt(1.0d - epssq * sc.sin() * sc.sin());

  290.         // compute ECEF Cartesian coordinates of specified point (for longitude=0)
  291.         final double xp = (rc + heightAboveEllipsoid) * sc.cos();
  292.         final double zp = (rc * (1.0d - epssq) + heightAboveEllipsoid) * sc.sin();

  293.         // compute spherical radius and angle lambda and phi of specified point
  294.         final double r = FastMath.hypot(xp, zp);
  295.         return new SphericalCoordinates(r, gp.getLongitude(), FastMath.asin(zp / r));
  296.     }

  297.     /** Rotate the magnetic vectors to geodetic coordinates.
  298.      * @param sph the spherical coordinates
  299.      * @param gp the geodetic point
  300.      * @param field the magnetic field in spherical coordinates
  301.      * @return the magnetic field in geodetic coordinates
  302.      */
  303.     private Vector3D rotateMagneticVector(final SphericalCoordinates sph,
  304.                                           final GeodeticPoint gp,
  305.                                           final Vector3D field) {

  306.         // difference between the spherical and geodetic latitudes
  307.         final double psi = sph.phi - gp.getLatitude();
  308.         final SinCos sc  = FastMath.sinCos(psi);

  309.         // rotate spherical field components to the geodetic system
  310.         final double Bz = field.getX() * sc.sin() + field.getZ() * sc.cos();
  311.         final double Bx = field.getX() * sc.cos() - field.getZ() * sc.sin();
  312.         final double By = field.getY();

  313.         return new Vector3D(Bx, By, Bz);
  314.     }

  315.     /** Computes Geomagnetic Field Elements X, Y and Z in spherical coordinate
  316.      * system using spherical harmonic summation.
  317.      * The vector Magnetic field is given by -grad V, where V is geomagnetic
  318.      * scalar potential. The gradient in spherical coordinates is given by:
  319.      * <pre>
  320.      *          dV ^   1 dV ^       1    dV ^
  321.      * grad V = -- r + - -- t + -------- -- p
  322.      *          dr     r dt     r sin(t) dp
  323.      * </pre>
  324.      * @param sph the spherical coordinates
  325.      * @param vars the spherical harmonic variables
  326.      * @param legendre the legendre function
  327.      * @return the magnetic field vector in spherical coordinates
  328.      */
  329.     private Vector3D summation(final SphericalCoordinates sph, final SphericalHarmonicVars vars,
  330.                                final LegendreFunction legendre) {

  331.         int index;
  332.         double Bx = 0.0;
  333.         double By = 0.0;
  334.         double Bz = 0.0;

  335.         for (int n = 1; n <= maxN; n++) {
  336.             for (int m = 0; m <= n; m++) {
  337.                 index = n * (n + 1) / 2 + m;

  338.                 /**
  339.                  * <pre>
  340.                  *       nMax               (n+2)   n    m            m           m
  341.                  * Bz = -SUM (n + 1) * (a/r)     * SUM [g cos(m p) + h sin(m p)] P (sin(phi))
  342.                  *       n=1                       m=0   n            n           n
  343.                  * </pre>
  344.                  * Equation 12 in the WMM Technical report. Derivative with respect to radius.
  345.                  */
  346.                 Bz -= vars.relativeRadiusPower[n] *
  347.                       (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * (1d + n) * legendre.mP[index];

  348.                 /**
  349.                  * <pre>
  350.                  *      nMax     (n+2)   n    m            m            m
  351.                  * By = SUM (a/r)     * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
  352.                  *      n=1             m=0   n            n            n
  353.                  * </pre>
  354.                  * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
  355.                  */
  356.                 By += vars.relativeRadiusPower[n] *
  357.                       (g[index] * vars.smLambda[m] - h[index] * vars.cmLambda[m]) * (double) m * legendre.mP[index];
  358.                 /**
  359.                  * <pre>
  360.                  *        nMax     (n+2)   n    m            m            m
  361.                  * Bx = - SUM (a/r)     * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
  362.                  *        n=1             m=0   n            n            n
  363.                  * </pre>
  364.                  * Equation 10 in the WMM Technical report. Derivative with respect to latitude, divided by radius.
  365.                  */
  366.                 Bx -= vars.relativeRadiusPower[n] *
  367.                       (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * legendre.mPDeriv[index];
  368.             }
  369.         }

  370.         final double cosPhi = FastMath.cos(sph.phi);
  371.         if (FastMath.abs(cosPhi) > 1.0e-10) {
  372.             By = By / cosPhi;
  373.         } else {
  374.             // special calculation for component - By - at geographic poles.
  375.             // To avoid using this function, make sure that the latitude is not
  376.             // exactly +/- π/2.
  377.             By = summationSpecial(sph, vars);
  378.         }

  379.         return new Vector3D(Bx, By, Bz);
  380.     }

  381.     /** Special calculation for the component By at geographic poles.
  382.      * @param sph the spherical coordinates
  383.      * @param vars the spherical harmonic variables
  384.      * @return the By component of the magnetic field
  385.      */
  386.     private double summationSpecial(final SphericalCoordinates sph, final SphericalHarmonicVars vars) {

  387.         double k;
  388.         final double sinPhi = FastMath.sin(sph.phi);
  389.         final double[] mPcupS = new double[maxN + 1];
  390.         mPcupS[0] = 1;
  391.         double By = 0.0;

  392.         for (int n = 1; n <= maxN; n++) {
  393.             final int index = n * (n + 1) / 2 + 1;
  394.             if (n == 1) {
  395.                 mPcupS[n] = mPcupS[n - 1];
  396.             } else {
  397.                 k = (double) (((n - 1) * (n - 1)) - 1) / (double) ((2 * n - 1) * (2 * n - 3));
  398.                 mPcupS[n] = sinPhi * mPcupS[n - 1] - k * mPcupS[n - 2];
  399.             }

  400.             /**
  401.              * <pre>
  402.              *      nMax     (n+2)   n    m            m            m
  403.              * By = SUM (a/r)     * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
  404.              *      n=1             m=0   n            n            n
  405.              * </pre>
  406.              * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
  407.              */
  408.             By += vars.relativeRadiusPower[n] *
  409.                   (g[index] * vars.smLambda[1] - h[index] * vars.cmLambda[1]) * mPcupS[n] * schmidtQuasiNorm[index];
  410.         }

  411.         return By;
  412.     }

  413.     /** Utility class to hold spherical coordinates. */
  414.     private static class SphericalCoordinates {

  415.         /** the radius (m). */
  416.         private double r;

  417.         /** the azimuth angle (radians). */
  418.         private double lambda;

  419.         /** the polar angle (radians). */
  420.         private double phi;

  421.         /** Create a new spherical coordinate object.
  422.          * @param r the radius, meters
  423.          * @param lambda the lambda angle, radians
  424.          * @param phi the phi angle, radians
  425.          */
  426.         private SphericalCoordinates(final double r, final double lambda, final double phi) {
  427.             this.r = r;
  428.             this.lambda = lambda;
  429.             this.phi = phi;
  430.         }
  431.     }

  432.     /** Utility class to compute certain variables for magnetic field summation. */
  433.     private class SphericalHarmonicVars {

  434.         /** (Radius of Earth / Spherical radius r)^(n+2). */
  435.         private double[] relativeRadiusPower;

  436.         /** cos(m*lambda). */
  437.         private double[] cmLambda;

  438.         /** sin(m*lambda). */
  439.         private double[] smLambda;

  440.         /** Calculates the spherical harmonic variables for a given spherical coordinate.
  441.          * @param sph the spherical coordinate
  442.          */
  443.         private SphericalHarmonicVars(final SphericalCoordinates sph) {

  444.             relativeRadiusPower = new double[maxN + 1];

  445.             // Compute a table of (EARTH_REFERENCE_RADIUS / radius)^n for i in
  446.             // 0 .. maxN (this is much faster than calling FastMath.pow maxN+1 times).

  447.             final double p = ellipsoidRadius / sph.r;
  448.             relativeRadiusPower[0] = p * p;
  449.             for (int n = 1; n <= maxN; n++) {
  450.                 relativeRadiusPower[n] = relativeRadiusPower[n - 1] * (ellipsoidRadius / sph.r);
  451.             }

  452.             // Compute tables of sin(lon * m) and cos(lon * m) for m = 0 .. maxN
  453.             // this is much faster than calling FastMath.sin and FastMath.cos maxN+1 times.

  454.             cmLambda = new double[maxN + 1];
  455.             smLambda = new double[maxN + 1];

  456.             cmLambda[0] = 1.0d;
  457.             smLambda[0] = 0.0d;

  458.             final SinCos scLambda = FastMath.sinCos(sph.lambda);
  459.             cmLambda[1] = scLambda.cos();
  460.             smLambda[1] = scLambda.sin();

  461.             for (int m = 2; m <= maxN; m++) {
  462.                 cmLambda[m] = cmLambda[m - 1] * scLambda.cos() - smLambda[m - 1] * scLambda.sin();
  463.                 smLambda[m] = cmLambda[m - 1] * scLambda.sin() + smLambda[m - 1] * scLambda.cos();
  464.             }
  465.         }
  466.     }

  467.     /** Utility class to compute a table of Schmidt-semi normalized associated Legendre functions. */
  468.     private class LegendreFunction {

  469.         /** the vector of all associated Legendre polynomials. */
  470.         private double[] mP;

  471.         /** the vector of derivatives of the Legendre polynomials wrt latitude. */
  472.         private double[] mPDeriv;

  473.         /** Calculate the Schmidt-semi normalized Legendre function.
  474.          * <p>
  475.          * <b>Note:</b> In geomagnetism, the derivatives of ALF are usually
  476.          * found with respect to the colatitudes. Here the derivatives are found
  477.          * with respect to the latitude. The difference is a sign reversal for
  478.          * the derivative of the Associated Legendre Functions.
  479.          * </p>
  480.          * @param x sinus of the spherical latitude (or cosinus of the spherical colatitude)
  481.          */
  482.         private LegendreFunction(final double x) {

  483.             final int numTerms = (maxN + 1) * (maxN + 2) / 2;

  484.             mP = new double[numTerms + 1];
  485.             mPDeriv = new double[numTerms + 1];

  486.             mP[0] = 1.0;
  487.             mPDeriv[0] = 0.0;

  488.             // sin (geocentric latitude) - sin_phi
  489.             final double z = FastMath.sqrt((1.0d - x) * (1.0d + x));

  490.             int index;
  491.             int index1;
  492.             int index2;

  493.             // First, compute the Gauss-normalized associated Legendre functions
  494.             for (int n = 1; n <= maxN; n++) {
  495.                 for (int m = 0; m <= n; m++) {
  496.                     index = n * (n + 1) / 2 + m;
  497.                     if (n == m) {
  498.                         index1 = (n - 1) * n / 2 + m - 1;
  499.                         mP[index] = z * mP[index1];
  500.                         mPDeriv[index] = z * mPDeriv[index1] + x * mP[index1];
  501.                     } else if (n == 1 && m == 0) {
  502.                         index1 = (n - 1) * n / 2 + m;
  503.                         mP[index] = x * mP[index1];
  504.                         mPDeriv[index] = x * mPDeriv[index1] - z * mP[index1];
  505.                     } else if (n > 1) {
  506.                         index1 = (n - 2) * (n - 1) / 2 + m;
  507.                         index2 = (n - 1) * n / 2 + m;
  508.                         if (m > n - 2) {
  509.                             mP[index] = x * mP[index2];
  510.                             mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2];
  511.                         } else {
  512.                             final double k = (double) ((n - 1) * (n - 1) - (m * m)) /
  513.                                              (double) ((2 * n - 1) * (2 * n - 3));

  514.                             mP[index] = x * mP[index2] - k * mP[index1];
  515.                             mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2] - k * mPDeriv[index1];
  516.                         }
  517.                     }

  518.                 }
  519.             }

  520.             // Converts the Gauss-normalized associated Legendre functions to the Schmidt quasi-normalized
  521.             // version using pre-computed relation stored in the variable schmidtQuasiNorm

  522.             for (int n = 1; n <= maxN; n++) {
  523.                 for (int m = 0; m <= n; m++) {
  524.                     index = n * (n + 1) / 2 + m;

  525.                     mP[index] = mP[index] * schmidtQuasiNorm[index];
  526.                     // The sign is changed since the new WMM routines use derivative with
  527.                     // respect to latitude instead of co-latitude
  528.                     mPDeriv[index] = -mPDeriv[index] * schmidtQuasiNorm[index];
  529.                 }
  530.             }
  531.         }
  532.     }
  533. }