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.orekit.bodies.GeodeticPoint;
  21. import org.orekit.errors.OrekitException;
  22. import org.orekit.errors.OrekitMessages;
  23. import org.orekit.utils.Constants;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  85.         this.validityStart = validityStart;
  86.         this.validityEnd = validityEnd;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  225.         return transformed;
  226.     }

  227.     /** Time transform the model coefficients from the base year of the model
  228.      * using a linear interpolation with a second model. The second model is
  229.      * required to have an adjacent validity period.
  230.      * @param otherModel the other magnetic field model
  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 GeoMagneticField otherModel, final double year) {

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

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

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

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

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

  263.         return transformed;
  264.     }

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

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

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

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

  284.         final double lat = gp.getLatitude();
  285.         final double heightAboveEllipsoid = gp.getAltitude();
  286.         final double sinLat = FastMath.sin(lat);

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

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

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

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

  305.         // difference between the spherical and geodetic latitudes
  306.         final double psi = sph.phi - gp.getLatitude();

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

  311.         return new Vector3D(Bx, By, Bz);
  312.     }

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

  329.         int index;
  330.         double Bx = 0.0;
  331.         double By = 0.0;
  332.         double Bz = 0.0;

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

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

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

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

  377.         return new Vector3D(Bx, By, Bz);
  378.     }

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

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

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

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

  409.         return By;
  410.     }

  411.     /** Utility class to hold spherical coordinates. */
  412.     private static class SphericalCoordinates {

  413.         /** the radius (m). */
  414.         private double r;

  415.         /** the azimuth angle (radians). */
  416.         private double lambda;

  417.         /** the polar angle (radians). */
  418.         private double phi;

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

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

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

  434.         /** cos(m*lambda). */
  435.         private double[] cmLambda;

  436.         /** sin(m*lambda). */
  437.         private double[] smLambda;

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

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

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

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

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

  452.             cmLambda = new double[maxN + 1];
  453.             smLambda = new double[maxN + 1];

  454.             cmLambda[0] = 1.0d;
  455.             smLambda[0] = 0.0d;

  456.             final double cosLambda = FastMath.cos(sph.lambda);
  457.             final double sinLambda = FastMath.sin(sph.lambda);
  458.             cmLambda[1] = cosLambda;
  459.             smLambda[1] = sinLambda;

  460.             for (int m = 2; m <= maxN; m++) {
  461.                 cmLambda[m] = cmLambda[m - 1] * cosLambda - smLambda[m - 1] * sinLambda;
  462.                 smLambda[m] = cmLambda[m - 1] * sinLambda + smLambda[m - 1] * cosLambda;
  463.             }
  464.         }
  465.     }

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

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

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

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

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

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

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

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

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

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

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

  517.                 }
  518.             }

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

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

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