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;

  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;

  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 {

  41.     /** Semi major-axis of WGS-84 ellipsoid in km. */
  42.     private static double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS / 1000d;

  44.     /** The first eccentricity squared. */
  45.     private static double epssq = 0.0066943799901413169961;

  47.     /** Mean radius of IAU-66 ellipsoid, in km. */
  48.     private static double ellipsoidRadius = 6371.2;

  50.     /** The model name. */
  51.     private String modelName;

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

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

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

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

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

  68.     /** maximum degree of spherical harmonic model. */
  69.     private int maxN;

  71.     /** maximum degree of spherical harmonic secular variations. */
  72.     private int maxNSec;

  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;

  79.     /** Pre-calculated ratio between gauss-normalized and schmidt quasi-normalized
  80.      * associated Legendre functions.
  81.      */
  82.     private double[] schmidtQuasiNorm;

  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) {

  98.         this.modelName = modelName;
  99.         this.epoch = epoch;
  100.         this.maxN = maxN;
  101.         this.maxNSec = maxNSec;

  103.         this.validityStart = validityStart;
  104.         this.validityEnd = validityEnd;

  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];

  111.         final int maxSecularFieldTerms = (maxNSec + 1) * (maxNSec + 2) / 2;
  112.         dg = new double[maxSecularFieldTerms];
  113.         dh = new double[maxSecularFieldTerms];

  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.

  118.         schmidtQuasiNorm = new double[maxMainFieldTerms + 1];
  119.         schmidtQuasiNorm[0] = 1.0;

  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;

  127.             // for m = 0
  128.             schmidtQuasiNorm[index] =
  129.                 schmidtQuasiNorm[index1] * (double) (2 * n - 1) / (double) n;

  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.     }

  141.     /** Returns the epoch for this magnetic field model.
  142.      * @return the epoch
  143.      */
  144.     public double getEpoch() {
  145.         return epoch;
  146.     }

  148.     /** Returns the model name.
  149.      * @return the model name
  150.      */
  151.     public String getModelName() {
  152.         return modelName;
  153.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  203.     /** Calculate the magnetic field at the specified geodetic point identified
  204.      * by latitude, longitude and altitude.
  205.      * @param latitude the WGS84 latitude in decimal degrees
  206.      * @param longitude the WGS84 longitude in decimal degrees
  207.      * @param height the height above the WGS84 ellipsoid in kilometers
  208.      * @return the {@link GeoMagneticElements} at the given geodetic point
  209.      */
  210.     public GeoMagneticElements calculateField(final double latitude,
  211.                                               final double longitude,
  212.                                               final double height) {

  214.         final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(latitude),
  215.                                                    FastMath.toRadians(longitude),
  216.                                                    height * 1000d);

  218.         final SphericalCoordinates sph = transformToSpherical(gp);
  219.         final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
  220.         final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));

  222.         // sum up the magnetic field vector components
  223.         final Vector3D magFieldSph = summation(sph, vars, legendre);
  224.         // rotate the field to geodetic coordinates
  225.         final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
  226.         // return the magnetic elements
  227.         return new GeoMagneticElements(magFieldGeo);
  228.     }

  230.     /** Time transform the model coefficients from the base year of the model
  231.      * using secular variation coefficients.
  232.      * @param year the year to which the model shall be transformed
  233.      * @return a time-transformed magnetic field model
  234.      */
  235.     public GeoMagneticField transformModel(final double year) {

  237.         if (!supportsTimeTransform()) {
  238.             throw new OrekitException(OrekitMessages.UNSUPPORTED_TIME_TRANSFORM, modelName, String.valueOf(epoch));
  239.         }

  241.         // the model can only be transformed within its validity period
  242.         if (year < validityStart || year > validityEnd) {
  243.             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
  244.                                       modelName, String.valueOf(epoch), year, validityStart, validityEnd);
  245.         }

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

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

  253.         for (int n = 1; n <= maxN; n++) {
  254.             for (int m = 0; m <= n; m++) {
  255.                 final int index = n * (n + 1) / 2 + m;
  256.                 if (index <= maxSecIndex) {
  257.                     transformed.h[index] = h[index] + dt * dh[index];
  258.                     transformed.g[index] = g[index] + dt * dg[index];
  259.                     // we need a copy of the secular var coef to calculate secular change
  260.                     transformed.dh[index] = dh[index];
  261.                     transformed.dg[index] = dg[index];
  262.                 } else {
  263.                     // just copy the parts that do not have corresponding secular variation coefficients
  264.                     transformed.h[index] = h[index];
  265.                     transformed.g[index] = g[index];
  266.                 }
  267.             }
  268.         }

  270.         return transformed;
  271.     }

  273.     /** Time transform the model coefficients from the base year of the model
  274.      * using a linear interpolation with a second model. The second model is
  275.      * required to have an adjacent validity period.
  276.      * @param otherModel the other magnetic field model
  277.      * @param year the year to which the model shall be transformed
  278.      * @return a time-transformed magnetic field model
  279.      */
  280.     public GeoMagneticField transformModel(final GeoMagneticField otherModel, final double year) {

  282.         // the model can only be transformed within its validity period
  283.         if (year < validityStart || year > validityEnd) {
  284.             throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
  285.                                       modelName, String.valueOf(epoch), year, validityStart, validityEnd);
  286.         }

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

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

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

  297.         for (int n = 1; n <= newMaxN; n++) {
  298.             for (int m = 0; m <= n; m++) {
  299.                 final int index = n * (n + 1) / 2 + m;
  300.                 if (index <= maxNCommonIndex) {
  301.                     transformed.h[index] = h[index] + factor * (otherModel.h[index] - h[index]);
  302.                     transformed.g[index] = g[index] + factor * (otherModel.g[index] - g[index]);
  303.                 } else {
  304.                     if (maxN < otherModel.maxN) {
  305.                         transformed.h[index] = factor * otherModel.h[index];
  306.                         transformed.g[index] = factor * otherModel.g[index];
  307.                     } else {
  308.                         transformed.h[index] = h[index] + factor * -h[index];
  309.                         transformed.g[index] = g[index] + factor * -g[index];
  310.                     }
  311.                 }
  312.             }
  313.         }

  315.         return transformed;
  316.     }

  318.     /** Utility function to get a decimal year for a given day.
  319.      * @param day the day (1-31)
  320.      * @param month the month (1-12)
  321.      * @param year the year
  322.      * @return the decimal year represented by the given day
  323.      */
  324.     public static double getDecimalYear(final int day, final int month, final int year) {
  325.         final int[] days = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
  326.         final int leapYear = (((year % 4) == 0) && (((year % 100) != 0) || ((year % 400) == 0))) ? 1 : 0;

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

  332.     /** Transform geodetic coordinates to spherical coordinates.
  333.      * @param gp the geodetic point
  334.      * @return the spherical coordinates wrt to the reference ellipsoid of the model
  335.      */
  336.     private SphericalCoordinates transformToSpherical(final GeodeticPoint gp) {

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

  341.         final double lat = gp.getLatitude();
  342.         final double heightAboveEllipsoid = gp.getAltitude() / 1000d;
  343.         final double sinLat = FastMath.sin(lat);

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

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

  352.         // compute spherical radius and angle lambda and phi of specified point
  353.         final double r = FastMath.hypot(xp, zp);
  354.         return new SphericalCoordinates(r, gp.getLongitude(), FastMath.asin(zp / r));
  355.     }

  357.     /** Rotate the magnetic vectors to geodetic coordinates.
  358.      * @param sph the spherical coordinates
  359.      * @param gp the geodetic point
  360.      * @param field the magnetic field in spherical coordinates
  361.      * @return the magnetic field in geodetic coordinates
  362.      */
  363.     private Vector3D rotateMagneticVector(final SphericalCoordinates sph,
  364.                                           final GeodeticPoint gp,
  365.                                           final Vector3D field) {

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

  370.         // rotate spherical field components to the geodetic system
  371.         final double Bz = field.getX() * FastMath.sin(psi) + field.getZ() * FastMath.cos(psi);
  372.         final double Bx = field.getX() * FastMath.cos(psi) - field.getZ() * FastMath.sin(psi);
  373.         final double By = field.getY();

  375.         return new Vector3D(Bx, By, Bz);
  376.     }

  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) {

  395.         int index;
  396.         double Bx = 0.0;
  397.         double By = 0.0;
  398.         double Bz = 0.0;

  400.         for (int n = 1; n <= maxN; n++) {
  401.             for (int m = 0; m <= n; m++) {
  402.                 index = n * (n + 1) / 2 + m;

  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];

  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.         }

  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 +/-90.
  445.             By = summationSpecial(sph, vars);
  446.         }

  448.         return new Vector3D(Bx, By, Bz);
  449.     }

  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) {

  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;

  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.             }

  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.         }

  485.         return By;
  486.     }

  488.     /** Utility class to hold spherical coordinates. */
  489.     private static class SphericalCoordinates {

  491.         /** the radius. */
  492.         private double r;

  494.         /** the azimuth angle. */
  495.         private double lambda;

  497.         /** the polar angle. */
  498.         private double phi;

  500.         /** Create a new spherical coordinate object.
  501.          * @param r the radius
  502.          * @param lambda the lambda angle
  503.          * @param phi the phi angle
  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.     }

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

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

  518.         /** cos(m*lambda). */
  519.         private double[] cmLambda;

  521.         /** sin(m*lambda). */
  522.         private double[] smLambda;

  524.         /** Calculates the spherical harmonic variables for a given spherical coordinate.
  525.          * @param sph the spherical coordinate
  526.          */
  527.         private SphericalHarmonicVars(final SphericalCoordinates sph) {

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

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

  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.             }

  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.

  543.             cmLambda = new double[maxN + 1];
  544.             smLambda = new double[maxN + 1];

  546.             cmLambda[0] = 1.0d;
  547.             smLambda[0] = 0.0d;

  549.             final double cosLambda = FastMath.cos(sph.lambda);
  550.             final double sinLambda = FastMath.sin(sph.lambda);
  551.             cmLambda[1] = cosLambda;
  552.             smLambda[1] = sinLambda;

  554.             for (int m = 2; m <= maxN; m++) {
  555.                 cmLambda[m] = cmLambda[m - 1] * cosLambda - smLambda[m - 1] * sinLambda;
  556.                 smLambda[m] = cmLambda[m - 1] * sinLambda + smLambda[m - 1] * cosLambda;
  557.             }
  558.         }
  559.     }

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

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

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

  570.         /** Calculate the Schmidt-semi normalized Legendre function.
  571.          * <p>
  572.          * <b>Note:</b> In geomagnetism, the derivatives of ALF are usually
  573.          * found with respect to the colatitudes. Here the derivatives are found
  574.          * with respect to the latitude. The difference is a sign reversal for
  575.          * the derivative of the Associated Legendre Functions.
  576.          * </p>
  577.          * @param x sinus of the spherical latitude (or cosinus of the spherical colatitude)
  578.          */
  579.         private LegendreFunction(final double x) {

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

  583.             mP = new double[numTerms + 1];
  584.             mPDeriv = new double[numTerms + 1];

  586.             mP[0] = 1.0;
  587.             mPDeriv[0] = 0.0;

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

  592.             int index;
  593.             int index1;
  594.             int index2;

  596.             // First, compute the Gauss-normalized associated Legendre functions
  597.             for (int n = 1; n <= maxN; n++) {
  598.                 for (int m = 0; m <= n; m++) {
  599.                     index = n * (n + 1) / 2 + m;
  600.                     if (n == m) {
  601.                         index1 = (n - 1) * n / 2 + m - 1;
  602.                         mP[index] = z * mP[index1];
  603.                         mPDeriv[index] = z * mPDeriv[index1] + x * mP[index1];
  604.                     } else if (n == 1 && m == 0) {
  605.                         index1 = (n - 1) * n / 2 + m;
  606.                         mP[index] = x * mP[index1];
  607.                         mPDeriv[index] = x * mPDeriv[index1] - z * mP[index1];
  608.                     } else if (n > 1 && n != m) {
  609.                         index1 = (n - 2) * (n - 1) / 2 + m;
  610.                         index2 = (n - 1) * n / 2 + m;
  611.                         if (m > n - 2) {
  612.                             mP[index] = x * mP[index2];
  613.                             mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2];
  614.                         } else {
  615.                             final double k = (double) ((n - 1) * (n - 1) - (m * m)) /
  616.                                              (double) ((2 * n - 1) * (2 * n - 3));

  618.                             mP[index] = x * mP[index2] - k * mP[index1];
  619.                             mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2] - k * mPDeriv[index1];
  620.                         }
  621.                     }

  623.                 }
  624.             }

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

  629.             for (int n = 1; n <= maxN; n++) {
  630.                 for (int m = 0; m <= n; m++) {
  631.                     index = n * (n + 1) / 2 + m;

  633.                     mP[index] = mP[index] * schmidtQuasiNorm[index];
  634.                     // The sign is changed since the new WMM routines use derivative with
  635.                     // respect to latitude instead of co-latitude
  636.                     mPDeriv[index] = -mPDeriv[index] * schmidtQuasiNorm[index];
  637.                 }
  638.             }
  639.         }
  640.     }
  641. }
Created with JaCoCo 0.8.2.201808211720