GeoMagneticField.java
- /* Copyright 2011-2012 Space Applications Services
- * Licensed to CS Communication & Systèmes (CS) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * CS licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- package org.orekit.models.earth;
- import org.hipparchus.geometry.euclidean.threed.Vector3D;
- import org.hipparchus.util.FastMath;
- import org.hipparchus.util.SinCos;
- import org.orekit.bodies.GeodeticPoint;
- import org.orekit.errors.OrekitException;
- import org.orekit.errors.OrekitMessages;
- import org.orekit.utils.Constants;
- /** Used to calculate the geomagnetic field at a given geodetic point on earth.
- * The calculation is estimated using spherical harmonic expansion of the
- * geomagnetic potential with coefficients provided by an actual geomagnetic
- * field model (e.g. IGRF, WMM).
- * <p>
- * Based on original software written by Manoj Nair from the National
- * Geophysical Data Center, NOAA, as part of the WMM 2010 software release
- * (WMM_SubLibrary.c)
- * </p>
- * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml">World Magnetic Model Overview</a>
- * @see <a href="http://www.ngdc.noaa.gov/geomag/WMM/soft.shtml">WMM Software Downloads</a>
- * @author Thomas Neidhart
- */
- public class GeoMagneticField {
- /** Semi major-axis of WGS-84 ellipsoid in m. */
- private static double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
- /** The first eccentricity squared. */
- private static double epssq = 0.0066943799901413169961;
- /** Mean radius of IAU-66 ellipsoid, in m. */
- private static double ellipsoidRadius = 6371200.0;
- /** The model name. */
- private String modelName;
- /** Base time of magnetic field model epoch (yrs). */
- private double epoch;
- /** C - Gauss coefficients of main geomagnetic model (nT). */
- private double[] g;
- /** C - Gauss coefficients of main geomagnetic model (nT). */
- private double[] h;
- /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
- private double[] dg;
- /** CD - Gauss coefficients of secular geomagnetic model (nT/yr). */
- private double[] dh;
- /** maximum degree of spherical harmonic model. */
- private int maxN;
- /** maximum degree of spherical harmonic secular variations. */
- private int maxNSec;
- /** the validity start of this magnetic field model. */
- private double validityStart;
- /** the validity end of this magnetic field model. */
- private double validityEnd;
- /** Pre-calculated ratio between gauss-normalized and schmidt quasi-normalized
- * associated Legendre functions.
- */
- private double[] schmidtQuasiNorm;
- /** Create a new geomagnetic field model with the given parameters. Internal
- * structures are initialized according to the specified degrees of the main
- * and secular variations.
- * @param modelName the model name
- * @param epoch the epoch of the model
- * @param maxN the maximum degree of the main model
- * @param maxNSec the maximum degree of the secular variations
- * @param validityStart validity start of this model
- * @param validityEnd validity end of this model
- */
- protected GeoMagneticField(final String modelName, final double epoch,
- final int maxN, final int maxNSec,
- final double validityStart, final double validityEnd) {
- this.modelName = modelName;
- this.epoch = epoch;
- this.maxN = maxN;
- this.maxNSec = maxNSec;
- this.validityStart = validityStart;
- this.validityEnd = validityEnd;
- // initialize main and secular field coefficient arrays
- final int maxMainFieldTerms = (maxN + 1) * (maxN + 2) / 2;
- g = new double[maxMainFieldTerms];
- h = new double[maxMainFieldTerms];
- final int maxSecularFieldTerms = (maxNSec + 1) * (maxNSec + 2) / 2;
- dg = new double[maxSecularFieldTerms];
- dh = new double[maxSecularFieldTerms];
- // pre-calculate the ratio between gauss-normalized and schmidt quasi-normalized
- // associated Legendre functions as they depend only on the degree of the model.
- schmidtQuasiNorm = new double[maxMainFieldTerms + 1];
- schmidtQuasiNorm[0] = 1.0;
- int index;
- int index1;
- for (int n = 1; n <= maxN; n++) {
- index = n * (n + 1) / 2;
- index1 = (n - 1) * n / 2;
- // for m = 0
- schmidtQuasiNorm[index] =
- schmidtQuasiNorm[index1] * (double) (2 * n - 1) / (double) n;
- for (int m = 1; m <= n; m++) {
- index = n * (n + 1) / 2 + m;
- index1 = n * (n + 1) / 2 + m - 1;
- schmidtQuasiNorm[index] =
- schmidtQuasiNorm[index1] *
- FastMath.sqrt((double) ((n - m + 1) * (m == 1 ? 2 : 1)) / (double) (n + m));
- }
- }
- }
- /** Returns the epoch for this magnetic field model.
- * @return the epoch
- */
- public double getEpoch() {
- return epoch;
- }
- /** Returns the model name.
- * @return the model name
- */
- public String getModelName() {
- return modelName;
- }
- /** Returns the start of the validity period for this model.
- * @return the validity start as decimal year
- */
- public double validFrom() {
- return validityStart;
- }
- /** Returns the end of the validity period for this model.
- * @return the validity end as decimal year
- */
- public double validTo() {
- return validityEnd;
- }
- /** Indicates whether this model supports time transformation or not.
- * @return <code>true</code> if this model can be transformed within its
- * validity period, <code>false</code> otherwise
- */
- public boolean supportsTimeTransform() {
- return maxNSec > 0;
- }
- /** Set the given main field coefficients.
- * @param n the n index
- * @param m the m index
- * @param gnm the g coefficient at position n,m
- * @param hnm the h coefficient at position n,m
- */
- protected void setMainFieldCoefficients(final int n, final int m,
- final double gnm, final double hnm) {
- final int index = n * (n + 1) / 2 + m;
- g[index] = gnm;
- h[index] = hnm;
- }
- /** Set the given secular variation coefficients.
- * @param n the n index
- * @param m the m index
- * @param dgnm the dg coefficient at position n,m
- * @param dhnm the dh coefficient at position n,m
- */
- protected void setSecularVariationCoefficients(final int n, final int m,
- final double dgnm, final double dhnm) {
- final int index = n * (n + 1) / 2 + m;
- dg[index] = dgnm;
- dh[index] = dhnm;
- }
- /** Calculate the magnetic field at the specified geodetic point identified
- * by latitude, longitude and altitude.
- * @param latitude the WGS84 latitude in radians
- * @param longitude the WGS84 longitude in radians
- * @param height the height above the WGS84 ellipsoid in meters
- * @return the {@link GeoMagneticElements} at the given geodetic point
- */
- public GeoMagneticElements calculateField(final double latitude,
- final double longitude,
- final double height) {
- final GeodeticPoint gp = new GeodeticPoint(latitude, longitude, height);
- final SphericalCoordinates sph = transformToSpherical(gp);
- final SphericalHarmonicVars vars = new SphericalHarmonicVars(sph);
- final LegendreFunction legendre = new LegendreFunction(FastMath.sin(sph.phi));
- // sum up the magnetic field vector components
- final Vector3D magFieldSph = summation(sph, vars, legendre);
- // rotate the field to geodetic coordinates
- final Vector3D magFieldGeo = rotateMagneticVector(sph, gp, magFieldSph);
- // return the magnetic elements
- return new GeoMagneticElements(magFieldGeo);
- }
- /** Time transform the model coefficients from the base year of the model
- * using secular variation coefficients.
- * @param year the year to which the model shall be transformed
- * @return a time-transformed magnetic field model
- */
- public GeoMagneticField transformModel(final double year) {
- if (!supportsTimeTransform()) {
- throw new OrekitException(OrekitMessages.UNSUPPORTED_TIME_TRANSFORM, modelName, String.valueOf(epoch));
- }
- // the model can only be transformed within its validity period
- if (year < validityStart || year > validityEnd) {
- throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
- modelName, String.valueOf(epoch), year, validityStart, validityEnd);
- }
- final double dt = year - epoch;
- final int maxSecIndex = maxNSec * (maxNSec + 1) / 2 + maxNSec;
- final GeoMagneticField transformed = new GeoMagneticField(modelName, year, maxN, maxNSec,
- validityStart, validityEnd);
- for (int n = 1; n <= maxN; n++) {
- for (int m = 0; m <= n; m++) {
- final int index = n * (n + 1) / 2 + m;
- if (index <= maxSecIndex) {
- transformed.h[index] = h[index] + dt * dh[index];
- transformed.g[index] = g[index] + dt * dg[index];
- // we need a copy of the secular var coef to calculate secular change
- transformed.dh[index] = dh[index];
- transformed.dg[index] = dg[index];
- } else {
- // just copy the parts that do not have corresponding secular variation coefficients
- transformed.h[index] = h[index];
- transformed.g[index] = g[index];
- }
- }
- }
- return transformed;
- }
- /** Time transform the model coefficients from the base year of the model
- * using a linear interpolation with a second model. The second model is
- * required to have an adjacent validity period.
- * @param otherModel the other magnetic field model
- * @param year the year to which the model shall be transformed
- * @return a time-transformed magnetic field model
- */
- public GeoMagneticField transformModel(final GeoMagneticField otherModel, final double year) {
- // the model can only be transformed within its validity period
- if (year < validityStart || year > validityEnd) {
- throw new OrekitException(OrekitMessages.OUT_OF_RANGE_TIME_TRANSFORM,
- modelName, String.valueOf(epoch), year, validityStart, validityEnd);
- }
- final double factor = (year - epoch) / (otherModel.epoch - epoch);
- final int maxNCommon = FastMath.min(maxN, otherModel.maxN);
- final int maxNCommonIndex = maxNCommon * (maxNCommon + 1) / 2 + maxNCommon;
- final int newMaxN = FastMath.max(maxN, otherModel.maxN);
- final GeoMagneticField transformed = new GeoMagneticField(modelName, year, newMaxN, 0,
- validityStart, validityEnd);
- for (int n = 1; n <= newMaxN; n++) {
- for (int m = 0; m <= n; m++) {
- final int index = n * (n + 1) / 2 + m;
- if (index <= maxNCommonIndex) {
- transformed.h[index] = h[index] + factor * (otherModel.h[index] - h[index]);
- transformed.g[index] = g[index] + factor * (otherModel.g[index] - g[index]);
- } else {
- if (maxN < otherModel.maxN) {
- transformed.h[index] = factor * otherModel.h[index];
- transformed.g[index] = factor * otherModel.g[index];
- } else {
- transformed.h[index] = h[index] + factor * -h[index];
- transformed.g[index] = g[index] + factor * -g[index];
- }
- }
- }
- }
- return transformed;
- }
- /** Utility function to get a decimal year for a given day.
- * @param day the day (1-31)
- * @param month the month (1-12)
- * @param year the year
- * @return the decimal year represented by the given day
- */
- public static double getDecimalYear(final int day, final int month, final int year) {
- final int[] days = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
- final int leapYear = ((year % 4) == 0 && ((year % 100) != 0 || (year % 400) == 0)) ? 1 : 0;
- final double dayInYear = days[month - 1] + (day - 1) + (month > 2 ? leapYear : 0);
- return (double) year + (dayInYear / (365.0d + leapYear));
- }
- /** Transform geodetic coordinates to spherical coordinates.
- * @param gp the geodetic point
- * @return the spherical coordinates wrt to the reference ellipsoid of the model
- */
- private SphericalCoordinates transformToSpherical(final GeodeticPoint gp) {
- // Convert geodetic coordinates (defined by the WGS-84 reference ellipsoid)
- // to Earth Centered Earth Fixed Cartesian coordinates, and then to spherical coordinates.
- final double lat = gp.getLatitude();
- final SinCos sc = FastMath.sinCos(lat);
- final double heightAboveEllipsoid = gp.getAltitude();
- // compute the local radius of curvature on the reference ellipsoid
- final double rc = a / FastMath.sqrt(1.0d - epssq * sc.sin() * sc.sin());
- // compute ECEF Cartesian coordinates of specified point (for longitude=0)
- final double xp = (rc + heightAboveEllipsoid) * sc.cos();
- final double zp = (rc * (1.0d - epssq) + heightAboveEllipsoid) * sc.sin();
- // compute spherical radius and angle lambda and phi of specified point
- final double r = FastMath.hypot(xp, zp);
- return new SphericalCoordinates(r, gp.getLongitude(), FastMath.asin(zp / r));
- }
- /** Rotate the magnetic vectors to geodetic coordinates.
- * @param sph the spherical coordinates
- * @param gp the geodetic point
- * @param field the magnetic field in spherical coordinates
- * @return the magnetic field in geodetic coordinates
- */
- private Vector3D rotateMagneticVector(final SphericalCoordinates sph,
- final GeodeticPoint gp,
- final Vector3D field) {
- // difference between the spherical and geodetic latitudes
- final double psi = sph.phi - gp.getLatitude();
- final SinCos sc = FastMath.sinCos(psi);
- // rotate spherical field components to the geodetic system
- final double Bz = field.getX() * sc.sin() + field.getZ() * sc.cos();
- final double Bx = field.getX() * sc.cos() - field.getZ() * sc.sin();
- final double By = field.getY();
- return new Vector3D(Bx, By, Bz);
- }
- /** Computes Geomagnetic Field Elements X, Y and Z in spherical coordinate
- * system using spherical harmonic summation.
- * The vector Magnetic field is given by -grad V, where V is geomagnetic
- * scalar potential. The gradient in spherical coordinates is given by:
- * <pre>
- * dV ^ 1 dV ^ 1 dV ^
- * grad V = -- r + - -- t + -------- -- p
- * dr r dt r sin(t) dp
- * </pre>
- * @param sph the spherical coordinates
- * @param vars the spherical harmonic variables
- * @param legendre the legendre function
- * @return the magnetic field vector in spherical coordinates
- */
- private Vector3D summation(final SphericalCoordinates sph, final SphericalHarmonicVars vars,
- final LegendreFunction legendre) {
- int index;
- double Bx = 0.0;
- double By = 0.0;
- double Bz = 0.0;
- for (int n = 1; n <= maxN; n++) {
- for (int m = 0; m <= n; m++) {
- index = n * (n + 1) / 2 + m;
- /**
- * <pre>
- * nMax (n+2) n m m m
- * Bz = -SUM (n + 1) * (a/r) * SUM [g cos(m p) + h sin(m p)] P (sin(phi))
- * n=1 m=0 n n n
- * </pre>
- * Equation 12 in the WMM Technical report. Derivative with respect to radius.
- */
- Bz -= vars.relativeRadiusPower[n] *
- (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * (1d + n) * legendre.mP[index];
- /**
- * <pre>
- * nMax (n+2) n m m m
- * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
- * n=1 m=0 n n n
- * </pre>
- * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
- */
- By += vars.relativeRadiusPower[n] *
- (g[index] * vars.smLambda[m] - h[index] * vars.cmLambda[m]) * (double) m * legendre.mP[index];
- /**
- * <pre>
- * nMax (n+2) n m m m
- * Bx = - SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
- * n=1 m=0 n n n
- * </pre>
- * Equation 10 in the WMM Technical report. Derivative with respect to latitude, divided by radius.
- */
- Bx -= vars.relativeRadiusPower[n] *
- (g[index] * vars.cmLambda[m] + h[index] * vars.smLambda[m]) * legendre.mPDeriv[index];
- }
- }
- final double cosPhi = FastMath.cos(sph.phi);
- if (FastMath.abs(cosPhi) > 1.0e-10) {
- By = By / cosPhi;
- } else {
- // special calculation for component - By - at geographic poles.
- // To avoid using this function, make sure that the latitude is not
- // exactly +/- π/2.
- By = summationSpecial(sph, vars);
- }
- return new Vector3D(Bx, By, Bz);
- }
- /** Special calculation for the component By at geographic poles.
- * @param sph the spherical coordinates
- * @param vars the spherical harmonic variables
- * @return the By component of the magnetic field
- */
- private double summationSpecial(final SphericalCoordinates sph, final SphericalHarmonicVars vars) {
- double k;
- final double sinPhi = FastMath.sin(sph.phi);
- final double[] mPcupS = new double[maxN + 1];
- mPcupS[0] = 1;
- double By = 0.0;
- for (int n = 1; n <= maxN; n++) {
- final int index = n * (n + 1) / 2 + 1;
- if (n == 1) {
- mPcupS[n] = mPcupS[n - 1];
- } else {
- k = (double) (((n - 1) * (n - 1)) - 1) / (double) ((2 * n - 1) * (2 * n - 3));
- mPcupS[n] = sinPhi * mPcupS[n - 1] - k * mPcupS[n - 2];
- }
- /**
- * <pre>
- * nMax (n+2) n m m m
- * By = SUM (a/r) * SUM [g cos(m p) + h sin(m p)] dP (sin(phi))
- * n=1 m=0 n n n
- * </pre>
- * Equation 11 in the WMM Technical report. Derivative with respect to longitude, divided by radius.
- */
- By += vars.relativeRadiusPower[n] *
- (g[index] * vars.smLambda[1] - h[index] * vars.cmLambda[1]) * mPcupS[n] * schmidtQuasiNorm[index];
- }
- return By;
- }
- /** Utility class to hold spherical coordinates. */
- private static class SphericalCoordinates {
- /** the radius (m). */
- private double r;
- /** the azimuth angle (radians). */
- private double lambda;
- /** the polar angle (radians). */
- private double phi;
- /** Create a new spherical coordinate object.
- * @param r the radius, meters
- * @param lambda the lambda angle, radians
- * @param phi the phi angle, radians
- */
- private SphericalCoordinates(final double r, final double lambda, final double phi) {
- this.r = r;
- this.lambda = lambda;
- this.phi = phi;
- }
- }
- /** Utility class to compute certain variables for magnetic field summation. */
- private class SphericalHarmonicVars {
- /** (Radius of Earth / Spherical radius r)^(n+2). */
- private double[] relativeRadiusPower;
- /** cos(m*lambda). */
- private double[] cmLambda;
- /** sin(m*lambda). */
- private double[] smLambda;
- /** Calculates the spherical harmonic variables for a given spherical coordinate.
- * @param sph the spherical coordinate
- */
- private SphericalHarmonicVars(final SphericalCoordinates sph) {
- relativeRadiusPower = new double[maxN + 1];
- // Compute a table of (EARTH_REFERENCE_RADIUS / radius)^n for i in
- // 0 .. maxN (this is much faster than calling FastMath.pow maxN+1 times).
- final double p = ellipsoidRadius / sph.r;
- relativeRadiusPower[0] = p * p;
- for (int n = 1; n <= maxN; n++) {
- relativeRadiusPower[n] = relativeRadiusPower[n - 1] * (ellipsoidRadius / sph.r);
- }
- // Compute tables of sin(lon * m) and cos(lon * m) for m = 0 .. maxN
- // this is much faster than calling FastMath.sin and FastMath.cos maxN+1 times.
- cmLambda = new double[maxN + 1];
- smLambda = new double[maxN + 1];
- cmLambda[0] = 1.0d;
- smLambda[0] = 0.0d;
- final SinCos scLambda = FastMath.sinCos(sph.lambda);
- cmLambda[1] = scLambda.cos();
- smLambda[1] = scLambda.sin();
- for (int m = 2; m <= maxN; m++) {
- cmLambda[m] = cmLambda[m - 1] * scLambda.cos() - smLambda[m - 1] * scLambda.sin();
- smLambda[m] = cmLambda[m - 1] * scLambda.sin() + smLambda[m - 1] * scLambda.cos();
- }
- }
- }
- /** Utility class to compute a table of Schmidt-semi normalized associated Legendre functions. */
- private class LegendreFunction {
- /** the vector of all associated Legendre polynomials. */
- private double[] mP;
- /** the vector of derivatives of the Legendre polynomials wrt latitude. */
- private double[] mPDeriv;
- /** Calculate the Schmidt-semi normalized Legendre function.
- * <p>
- * <b>Note:</b> In geomagnetism, the derivatives of ALF are usually
- * found with respect to the colatitudes. Here the derivatives are found
- * with respect to the latitude. The difference is a sign reversal for
- * the derivative of the Associated Legendre Functions.
- * </p>
- * @param x sinus of the spherical latitude (or cosinus of the spherical colatitude)
- */
- private LegendreFunction(final double x) {
- final int numTerms = (maxN + 1) * (maxN + 2) / 2;
- mP = new double[numTerms + 1];
- mPDeriv = new double[numTerms + 1];
- mP[0] = 1.0;
- mPDeriv[0] = 0.0;
- // sin (geocentric latitude) - sin_phi
- final double z = FastMath.sqrt((1.0d - x) * (1.0d + x));
- int index;
- int index1;
- int index2;
- // First, compute the Gauss-normalized associated Legendre functions
- for (int n = 1; n <= maxN; n++) {
- for (int m = 0; m <= n; m++) {
- index = n * (n + 1) / 2 + m;
- if (n == m) {
- index1 = (n - 1) * n / 2 + m - 1;
- mP[index] = z * mP[index1];
- mPDeriv[index] = z * mPDeriv[index1] + x * mP[index1];
- } else if (n == 1 && m == 0) {
- index1 = (n - 1) * n / 2 + m;
- mP[index] = x * mP[index1];
- mPDeriv[index] = x * mPDeriv[index1] - z * mP[index1];
- } else if (n > 1) {
- index1 = (n - 2) * (n - 1) / 2 + m;
- index2 = (n - 1) * n / 2 + m;
- if (m > n - 2) {
- mP[index] = x * mP[index2];
- mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2];
- } else {
- final double k = (double) ((n - 1) * (n - 1) - (m * m)) /
- (double) ((2 * n - 1) * (2 * n - 3));
- mP[index] = x * mP[index2] - k * mP[index1];
- mPDeriv[index] = x * mPDeriv[index2] - z * mP[index2] - k * mPDeriv[index1];
- }
- }
- }
- }
- // Converts the Gauss-normalized associated Legendre functions to the Schmidt quasi-normalized
- // version using pre-computed relation stored in the variable schmidtQuasiNorm
- for (int n = 1; n <= maxN; n++) {
- for (int m = 0; m <= n; m++) {
- index = n * (n + 1) / 2 + m;
- mP[index] = mP[index] * schmidtQuasiNorm[index];
- // The sign is changed since the new WMM routines use derivative with
- // respect to latitude instead of co-latitude
- mPDeriv[index] = -mPDeriv[index] * schmidtQuasiNorm[index];
- }
- }
- }
- }
- }