simgear/magvar/coremag.cxx

   1 // coremag.cxx -- compute local magnetic variation given position,
   2 //                altitude, and date
   3 //
   4 // This is an implementation of the NIMA (formerly DMA) WMM2000
   5 //
   6 //    http://www.nima.mil/GandG/ngdc-wmm2000.html
   7 //
   8 // Copyright (C) 2000  Edward A Williams <Ed_Williams@compuserve.com>
   9 //
  10 // Adapted from Excel 3.0 version 3/27/94 EAW
  11 // Recoded in C++ by Starry Chan
  12 // WMM95 added and rearranged in ANSI-C EAW 7/9/95
  13 // Put shell around program and made Borland & GCC compatible EAW 11/22/95
  14 // IGRF95 added 2/96 EAW
  15 // WMM2000 IGR2000 added 2/00 EAW
  16 // Released under GPL 3/26/00 EAW
  17 // Adaptions and modifications for the SimGear project  3/27/2000 CLO
  18 //
  19 // Removed all pow() calls and made static roots[][] arrays to
  20 // save many sqrt() calls on subsequent invocations
  21 // left old code as SGMagVarOrig() for testing purposes
  22 // 3/28/2000  Norman Vine -- nhv@yahoo.com
  23 //
  24 // Put in some bullet-proofing to handle magnetic and geographic poles.
  25 // 3/28/2000 EAW
  26 //
  27 // Updated coefficient arrays to use the current wmm2005 model,
  28 // (valid between 2005.0 and 2010.0)
  29 // Also removed unused variables and corrected earth radii constants
  30 // to the values for WGS84 and WMM2005.
  31 // Reference:
  32 //     McLean, S., S. Macmillan, S. Maus, V. Lesur, A.
  33 //     Thomson, and D. Dater, December 2004, The
  34 //     US/UK World Magnetic Model for 2005-2010,
  35 //     NOAA Technical Report NESDIS/NGDC-1.
  36 //
  37 // 25/10/2006  Wim Van Hoydonck -- wim.van.hoydonck@gmail.com
  38 //
  39 //
  40 // Updated coefficient arrays to use the current WMM2015 model,
  41 // (valid between 2015.0 and 2020.0)
  42 // Also removed unused variables and corrected earth radii constants
  43 // to the values for WGS84 and WMM2015.
  44 // Reference:
  45 //     A. Chulliat , S. Macmillan, P. Alken, C. Beggan, M.
  46 //     Nair, B. Hamilton, A. Woods, V. Ridley,
  47 //     S Maus, and A Thomson, December 2014, The
  48 //     US/UK World Magnetic Model for 2015-2020,
  49 //     NOAA Technical Report WMM2015_Report.pdf
  50 //
  51 // 18/06/2015  Jean-Paul Anceaux -- j.p.r.anceaux@gmail.com
  52
  53
  54 //  The routine uses a spherical harmonic expansion of the magnetic
  55 // potential up to twelfth order, together with its time variation, as
  56 // described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
  57 // Academic Press (London 1987)". The program first converts geodetic
  58 // coordinates (lat/long on elliptic earth and altitude) to spherical
  59 // geocentric (spherical lat/long and radius) coordinates. Using this,
  60 // the spherical (B_r, B_theta, B_phi) magnetic field components are
  61 // computed from the model. These are finally referred to surface (X, Y,
  62 // Z) coordinates.
  63 //
  64 //   Fields are accurate to better than 200nT, variation and dip to
  65 // better than 0.5 degrees, with the exception of the declination near
  66 // the magnetic poles (where it is ill-defined) where the error may reach
  67 // 4 degrees or more.
  68 //
  69 //   Variation is undefined at both the geographic and
  70 // magnetic poles, even though the field itself is well-behaved. To
  71 // avoid the routine blowing up, latitude entries corresponding to
  72 // the geographic poles are slightly offset. At the magnetic poles,
  73 // the routine returns zero variation.
  74
  75
  76 //
  77 // This library is free software; you can redistribute it and/or
  78 // modify it under the terms of the GNU Library General Public
  79 // License as published by the Free Software Foundation; either
  80 // version 2 of the License, or (at your option) any later version.
  81 //
  82 // This library is distributed in the hope that it will be useful,
  83 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  84 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  85 // Library General Public License for more details.
  86 //
  87 // You should have received a copy of the GNU General Public License
  88 // along with this program; if not, write to the Free Software
  89 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  90 //
  91 // $Id$
  92
  93
  94 #include <stdio.h>
  95 #include <stdlib.h>
  96 #include <cmath>
  97
  98 #include <simgear/constants.h>
  99 #include <simgear/sg_inlines.h>
 100
 101 #include "coremag.hxx"
 102
 103 static const double a = 6378.137;       /* semi-major axis (equatorial radius) of WGS84 ellipsoid */
 104 static const double b = 6356.7523142;   /* semi-minor axis referenced to the WGS84 ellipsoid */
 105 static const double r_0 = 6371.2;       /* standard Earth magnetic reference radius  */
 106
 107 static double gnm_wmm2015[13][13] =
 108 {
 109     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 110     {-29438.5, -1501.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 111     {-2445.3, 3012.5, 1676.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 112     {1351.1, -2352.3, 1225.6, 581.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 113     {907.2, 813.7, 120.3, -335.0, 70.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 114     {-232.6, 360.1, 192.4, -141.0, -157.4, 4.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 115     {69.5, 67.4, 72.8, -129.8, -29.0, 13.2, -70.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 116     {81.6, -76.1, -6.8, 51.9, 15.0, 9.3, -2.8, 6.7, 0.0, 0.0, 0.0, 0.0, 0.0},
 117     {24.0, 8.6, -16.9, -3.2, -20.6, 13.3, 11.7, -16.0, -2.0, 0.0, 0.0, 0.0, 0.0},
 118     {5.4, 8.8, 3.1, -3.1, 0.6, -13.3, -0.1, 8.7, -9.1, -10.5, 0.0, 0.0, 0.0},
 119     {-1.9, -6.5, 0.2, 0.6, -0.6, 1.7, -0.7, 2.1, 2.3, -1.8, -3.6, 0.0, 0.0},
 120     {3.1, -1.5, -2.3, 2.1, -0.9, 0.6, -0.7, 0.2, 1.7, -0.2, 0.4, 3.5, 0.0},
 121     {-2.0, -0.3, 0.4, 1.3, -0.9, 0.9, 0.1, 0.5, -0.4, -0.4, 0.2, -0.9, 0.0},
 122 };
 123
 124 static double hnm_wmm2015[13][13]=
 125 {
 126     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 127     {0.0, 4796.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 128     {0.0, -2845.6, -642.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 129     {0.0, -115.3, 245.0, -538.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 130     {0.0, 283.4, -188.6, 180.9, -329.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 131     {0.0, 47.4, 196.9, -119.4, 16.1, 100.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 132     {0.0, -20.7, 33.2, 58.8, -66.5, 7.3, 62.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 133     {0.0, -54.1, -19.4, 5.6, 24.4, 3.3, -27.5, -2.3, 0.0, 0.0, 0.0, 0.0, 0.0},
 134     {0.0, 10.2, -18.1, 13.2, -14.6, 16.2, 5.7, -9.1, 2.2, 0.0, 0.0, 0.0, 0.0},
 135     {0.0, -21.6, 10.8, 11.7, -6.8, -6.9, 7.8, 1.0, -3.9, 8.5, 0.0, 0.0, 0.0},
 136     {0.0, 3.3, -0.3, 4.6, 4.4, -7.9, -0.6, -4.1, -2.8, -1.1, -8.7, 0.0, 0.0},
 137     {0.0, -0.1, 2.1, -0.7, -1.1, 0.7, -0.2, -2.1, -1.5, -2.5, -2.0, -2.3, 0.0},
 138     {0.0, -1.0, 0.5, 1.8, -2.2, 0.3, 0.7, -0.1, 0.3, 0.2, -0.9, -0.2, 0.7},
 139 };
 140
 141 static double gtnm_wmm2015[13][13]=
 142 {
 143     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 144     {10.7, 17.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 145     {-8.6, -3.3, 2.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 146     {3.1, -6.2, -0.4, -10.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 147     {-0.4, 0.8, -9.2, 4.0, -4.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 148     {-0.2, 0.1, -1.4, 0.0, 1.3, 3.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 149     {-0.5, -0.2, -0.6, 2.4, -1.1, 0.3, 1.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 150     {0.2, -0.2, -0.4, 1.3, 0.2, -0.4, -0.9, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0},
 151     {0.0, 0.1, -0.5, 0.5, -0.2, 0.4, 0.2, -0.4, 0.3, 0.0, 0.0, 0.0, 0.0},
 152     {0.0, -0.1, -0.1, 0.4, -0.5, -0.2, 0.1, 0.0, -0.2, -0.1, 0.0, 0.0, 0.0},
 153     {0.0, 0.0, -0.1, 0.3, -0.1, -0.1, -0.1, 0.0, -0.2, -0.1, -0.2, 0.0, 0.0},
 154     {0.0, 0.0, -0.1, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.1, -0.1, 0.0},
 155     {0.1, 0.0, 0.0, 0.1, -0.1, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 156 };
 157
 158 static double htnm_wmm2015[13][13]=
 159 {
 160     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 161     {0.0, -26.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 162     {0.0, -27.1, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 163     {0.0, 8.4, -0.4, 2.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 164     {0.0, -0.6, 5.3, 3.0, -5.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 165     {0.0, 0.4, 1.6, -1.1, 3.3, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 166     {0.0, 0.0, -2.2, -0.7, 0.1, 1.0, 1.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 167     {0.0, 0.7, 0.5, -0.2, -0.1, -0.7, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
 168     {0.0, -0.3, 0.3, 0.3, 0.6, -0.1, -0.2, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0},
 169     {0.0, -0.2, -0.1, -0.2, 0.1, 0.1, 0.0, -0.2, 0.4, 0.3, 0.0, 0.0, 0.0},
 170     {0.0, 0.1, -0.1, 0.0, 0.0, -0.2, 0.1, -0.1, -0.2, 0.1, -0.1, 0.0, 0.0},
 171     {0.0, 0.0, 0.1, 0.0, 0.1, 0.0, 0.0, 0.1, 0.0, -0.1, 0.0, -0.1, 0.0},
 172     {0.0, 0.0, 0.0, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 173 };
 174
 175 static const int nmax = 12;
 176
 177 static double P[13][13];
 178 static double DP[13][13];
 179 static double gnm[13][13];
 180 static double hnm[13][13];
 181 static double sm[13];
 182 static double cm[13];
 183
 184 static double root[13];
 185 static double roots[13][13][2];
 186
 187 /* Convert date to Julian day    1950-2049 */
 188 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
 189 {
 190     unsigned long jd;
 191
 192     yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
 193     jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
 194     jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
 195     jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
 196
 197     /* printf("julian date = %d\n", jd ); */
 198     return jd;
 199 }
 200
 201
 202 /*
 203  * return variation (in radians) given geodetic latitude (radians),
 204  * longitude(radians), height (km) and (Julian) date
 205  * N and E lat and long are positive, S and W negative
 206 */
 207
 208 double calc_magvar( double lat, double lon, double h, long dat, double* field )
 209 {
 210     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 211     int n,m;
 212     /* reference date for current model is 1 januari 2015 */
 213     long date0_wmm2015 = yymmdd_to_julian_days(15,1,1);
 214
 215     double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
 216     double sinpsi, cospsi, inv_s;
 217
 218     static int been_here = 0;
 219
 220     double sinlat = sin(lat);
 221     double coslat = cos(lat);
 222
 223     /* convert to geocentric coords: */
 224     // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
 225     sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
 226     /* sr is effective radius */
 227     theta = atan2(coslat * (h*sr + a*a),
 228                   sinlat * (h*sr + b*b));
 229     /* theta is geocentric co-latitude */
 230
 231     r = h*h + 2.0*h * sr +
 232         (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
 233         (a*a - (a*a - b*b) * sinlat*sinlat );
 234
 235     r = sqrt(r);
 236
 237     /* r is geocentric radial distance */
 238     c = cos(theta);
 239     s = sin(theta);
 240     /* protect against zero divide at geographic poles */
 241     inv_s =  1.0 / (s + (s == 0.)*1.0e-8);
 242
 243     /* zero out arrays */
 244     for ( n = 0; n <= nmax; n++ ) {
 245         for ( m = 0; m <= n; m++ ) {
 246             P[n][m] = 0;
 247             DP[n][m] = 0;
 248         }
 249     }
 250
 251     /* diagonal elements */
 252     P[0][0] = 1;
 253     P[1][1] = s;
 254     DP[0][0] = 0;
 255     DP[1][1] = c;
 256     P[1][0] = c ;
 257     DP[1][0] = -s;
 258
 259     // these values will not change for subsequent function calls
 260     if( !been_here ) {
 261         for ( n = 2; n <= nmax; n++ ) {
 262             root[n] = sqrt((2.0*n-1) / (2.0*n));
 263         }
 264
 265         for ( m = 0; m <= nmax; m++ ) {
 266             double mm = m*m;
 267             for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 268                 roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
 269                 roots[m][n][1] = 1.0 / sqrt( n*n - mm);
 270             }
 271         }
 272         been_here = 1;
 273     }
 274
 275     for ( n=2; n <= nmax; n++ ) {
 276         // double root = sqrt((2.0*n-1) / (2.0*n));
 277         P[n][n] = P[n-1][n-1] * s * root[n];
 278         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 279             root[n];
 280     }
 281
 282     /* lower triangle */
 283     for ( m = 0; m <= nmax; m++ ) {
 284         // double mm = m*m;
 285         for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 286             // double root1 = sqrt((n-1)*(n-1) - mm);
 287             // double root2 = 1.0 / sqrt( n*n - mm);
 288             P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
 289                        P[n-2][m] * roots[m][n][0]) *
 290                 roots[m][n][1];
 291
 292             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 293                         (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) *
 294                 roots[m][n][1];
 295         }
 296     }
 297
 298     /* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
 299        achieved by adjusting the coefficients at time t0 for linear secular variation */
 300     /* WMM2015 */
 301     yearfrac = (dat - date0_wmm2015) / 365.25;
 302     for ( n = 1; n <= nmax; n++ ) {
 303         for ( m = 0; m <= nmax; m++ ) {
 304             gnm[n][m] = gnm_wmm2015[n][m] + yearfrac * gtnm_wmm2015[n][m];
 305             hnm[n][m] = hnm_wmm2015[n][m] + yearfrac * htnm_wmm2015[n][m];
 306         }
 307     }
 308
 309     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 310     for ( m = 0; m <= nmax; m++ ) {
 311         sm[m] = sin(m * lon);
 312         cm[m] = cos(m * lon);
 313     }
 314
 315     /* compute B fields */
 316     B_r = 0.0;
 317     B_theta = 0.0;
 318     B_phi = 0.0;
 319     fn_0 = r_0/r;
 320     fn = fn_0 * fn_0;
 321
 322     for ( n = 1; n <= nmax; n++ ) {
 323         double c1_n=0;
 324         double c2_n=0;
 325         double c3_n=0;
 326         for ( m = 0; m <= n; m++ ) {
 327             double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
 328             c1_n=c1_n + tmp * P[n][m];
 329             c2_n=c2_n + tmp * DP[n][m];
 330             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 331         }
 332         // fn=pow(r_0/r,n+2.0);
 333         fn *= fn_0;
 334         B_r = B_r + (n + 1) * c1_n * fn;
 335         B_theta = B_theta - c2_n * fn;
 336         B_phi = B_phi + c3_n * fn * inv_s;
 337     }
 338
 339     /* Find geodetic field components: */
 340     psi = theta - ((M_PI / 2.0) - lat);
 341     sinpsi = sin(psi);
 342     cospsi = cos(psi);
 343     X = -B_theta * cospsi - B_r * sinpsi;
 344     Y = B_phi;
 345     Z = B_theta * sinpsi - B_r * cospsi;
 346
 347     field[0]=B_r;
 348     field[1]=B_theta;
 349     field[2]=B_phi;
 350     field[3]=X;
 351     field[4]=Y;
 352     field[5]=Z;   /* output fields */
 353
 354     /* find variation in radians */
 355     /* return zero variation at magnetic pole X=Y=0. */
 356     /* E is positive */
 357     return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
 358 }
 359
 360
 361 #ifdef TEST_NHV_HACKS
 362 double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
 363 {
 364     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 365     int n,m;
 366     /* reference dates */
 367     long date0_wmm2015 = yymmdd_to_julian_days(5,1,1);
 368
 369     double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
 370
 371     /* convert to geocentric coords: */
 372     sr = sqrt(pow(a*cos(lat),2.0)+pow(b*sin(lat),2.0));
 373     /* sr is effective radius */
 374     theta = atan2(cos(lat) * (h * sr + a * a),
 375                   sin(lat) * (h * sr + b * b));
 376     /* theta is geocentric co-latitude */
 377
 378     r = h * h + 2.0*h * sr +
 379         (pow(a,4.0) - (pow(a,4.0) - pow(b,4.0)) * pow(sin(lat),2.0)) /
 380         (a * a - (a * a - b * b) * pow(sin(lat),2.0));
 381
 382     r = sqrt(r);
 383
 384     /* r is geocentric radial distance */
 385     c = cos(theta);
 386     s = sin(theta);
 387
 388     /* zero out arrays */
 389     for ( n = 0; n <= nmax; n++ ) {
 390         for ( m = 0; m <= n; m++ ) {
 391             P[n][m] = 0;
 392             DP[n][m] = 0;
 393         }
 394     }
 395
 396     /* diagonal elements */
 397     P[0][0] = 1;
 398     P[1][1] = s;
 399     DP[0][0] = 0;
 400     DP[1][1] = c;
 401     P[1][0] = c ;
 402     DP[1][0] = -s;
 403
 404     for ( n = 2; n <= nmax; n++ ) {
 405         P[n][n] = P[n-1][n-1] * s * sqrt((2.0*n-1) / (2.0*n));
 406         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 407             sqrt((2.0*n-1) / (2.0*n));
 408     }
 409
 410     /* lower triangle */
 411     for ( m = 0; m <= nmax; m++ ) {
 412         for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 413             P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
 414                        sqrt(1.0*(n-1)*(n-1) - m * m)) /
 415                 sqrt(1.0* n * n - m * m);
 416             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 417                         (2.0*n-1) - DP[n-2][m] *
 418                         sqrt(1.0*(n-1) * (n-1) - m * m)) /
 419                 sqrt(1.0* n * n - m * m);
 420         }
 421     }
 422
 423     /* compute gnm, hnm at dat */
 424     /* WMM2015 */
 425     yearfrac = (dat - date0_wmm2015) / 365.25;
 426     for ( n = 1; n <= nmax; n++ ) {
 427         for ( m = 0; m <= nmax; m++ ) {
 428             gnm[n][m] = gnm_wmm2015[n][m] + yearfrac * gtnm_wmm2015[n][m];
 429             hnm[n][m] = hnm_wmm2015[n][m] + yearfrac * htnm_wmm2015[n][m];
 430         }
 431     }
 432
 433     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 434     for ( m = 0; m <= nmax; m++ ) {
 435         sm[m] = sin(m * lon);
 436         cm[m] = cos(m * lon);
 437     }
 438
 439     /* compute B fields */
 440     B_r = 0.0;
 441     B_theta = 0.0;
 442     B_phi = 0.0;
 443
 444     for ( n = 1; n <= nmax; n++ ) {
 445         double c1_n=0;
 446         double c2_n=0;
 447         double c3_n=0;
 448         for ( m = 0; m <= n; m++ ) {
 449             c1_n=c1_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * P[n][m];
 450             c2_n=c2_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * DP[n][m];
 451             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 452         }
 453         fn=pow(r_0/r,n+2.0);
 454         B_r = B_r + (n + 1) * c1_n * fn;
 455         B_theta = B_theta - c2_n * fn;
 456         B_phi = B_phi + c3_n * fn / s;
 457     }
 458
 459     /* Find geodetic field components: */
 460     psi = theta - (pi / 2.0 - lat);
 461     X = -B_theta * cos(psi) - B_r * sin(psi);
 462     Y = B_phi;
 463     Z = B_theta * sin(psi) - B_r * cos(psi);
 464
 465     field[0]=B_r;
 466     field[1]=B_theta;
 467     field[2]=B_phi;
 468     field[3]=X;
 469     field[4]=Y;
 470     field[5]=Z;   /* output fields */
 471
 472     /* find variation, leave in radians! */
 473     return atan2(Y, X);  /* E is positive */
 474 }
 475 #endif // TEST_NHV_HACKS