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 //  The routine uses a spherical harmonic expansion of the magnetic
  41 // potential up to twelfth order, together with its time variation, as
  42 // described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
  43 // Academic Press (London 1987)". The program first converts geodetic
  44 // coordinates (lat/long on elliptic earth and altitude) to spherical
  45 // geocentric (spherical lat/long and radius) coordinates. Using this,
  46 // the spherical (B_r, B_theta, B_phi) magnetic field components are
  47 // computed from the model. These are finally referred to surface (X, Y,
  48 // Z) coordinates.
  49 //
  50 //   Fields are accurate to better than 200nT, variation and dip to
  51 // better than 0.5 degrees, with the exception of the declination near
  52 // the magnetic poles (where it is ill-defined) where the error may reach
  53 // 4 degrees or more.
  54 //
  55 //   Variation is undefined at both the geographic and
  56 // magnetic poles, even though the field itself is well-behaved. To
  57 // avoid the routine blowing up, latitude entries corresponding to
  58 // the geographic poles are slightly offset. At the magnetic poles,
  59 // the routine returns zero variation.
  60
  61
  62 //
  63 // This library is free software; you can redistribute it and/or
  64 // modify it under the terms of the GNU Library General Public
  65 // License as published by the Free Software Foundation; either
  66 // version 2 of the License, or (at your option) any later version.
  67 //
  68 // This library is distributed in the hope that it will be useful,
  69 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  70 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  71 // Library General Public License for more details.
  72 //
  73 // You should have received a copy of the GNU General Public License
  74 // along with this program; if not, write to the Free Software
  75 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  76 //
  77 // $Id$
  78
  79
  80 #include <stdio.h>
  81 #include <stdlib.h>
  82 #include <math.h>
  83
  84 #include <simgear/constants.h>
  85 #include <simgear/sg_inlines.h>
  86
  87 #include "coremag.hxx"
  88
  89
  90 static const double pi = 3.14159265358979;
  91 static const double a = 6378.137;       /* semi-major axis (equatorial radius) of WGS84 ellipsoid */
  92 static const double b = 6356.7523142;   /* semi-minor axis referenced to the WGS84 ellipsoid */
  93 static const double r_0 = 6371.2;       /* standard Earth magnetic reference radius  */
  94
  95 static double gnm_wmm2005[13][13] =
  96 {
  97     {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},
  98     {-29556.8, -1671.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  99     {-2340.6, 3046.9, 1657.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 100     {1335.4, -2305.1, 1246.7, 674.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 101     {919.8, 798.1, 211.3, -379.4, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 102     {-227.4, 354.6, 208.7, -136.5, -168.3, -14.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 103     {73.2, 69.7, 76.7, -151.2, -14.9, 14.6, -86.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 104     {80.1, -74.5, -1.4, 38.5, 12.4, 9.5, 5.7, 1.8, 0.0, 0.0, 0.0, 0.0, 0.0},
 105     {24.9, 7.7, -11.6, -6.9, -18.2, 10.0, 9.2, -11.6, -5.2, 0.0, 0.0, 0.0, 0.0},
 106     {5.6, 9.9, 3.5, -7.0, 5.1, -10.8, -1.3, 8.8, -6.7, -9.1, 0.0, 0.0, 0.0},
 107     {-2.3, -6.3, 1.6, -2.6, 0.0, 3.1, 0.4, 2.1, 3.9, -0.1, -2.3, 0.0, 0.0},
 108     {2.8, -1.6, -1.7, 1.7, -0.1, 0.1, -0.7, 0.7, 1.8, 0.0, 1.1, 4.1, 0.0},
 109     {-2.4, -0.4, 0.2, 0.8, -0.3, 1.1, -0.5, 0.4, -0.3, -0.3, -0.1, -0.3, -0.1},
 110 };
 111
 112 static double hnm_wmm2005[13][13]=
 113 {
 114     {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},
 115     {0.0, 5079.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 116     {0.0, -2594.7, -516.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 117     {0.0, -199.9, 269.3, -524.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 118     {0.0, 281.5, -226.0, 145.8, -304.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 119     {0.0, 42.4, 179.8, -123.0, -19.5, 103.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 120     {0.0, -20.3, 54.7, 63.6, -63.4, -0.1, 50.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 121     {0.0, -61.5, -22.4, 7.2, 25.4, 11.0, -26.4, -5.1, 0.0, 0.0, 0.0, 0.0, 0.0},
 122     {0.0, 11.2, -21.0, 9.6, -19.8, 16.1, 7.7, -12.9, -0.2, 0.0, 0.0, 0.0, 0.0},
 123     {0.0, -20.1, 12.9, 12.6, -6.7, -8.1, 8.0, 2.9, -7.9, 6.0, 0.0, 0.0, 0.0},
 124     {0.0, 2.4, 0.2, 4.4, 4.8, -6.5, -1.1, -3.4, -0.8, -2.3, -7.9, 0.0, 0.0},
 125     {0.0, 0.3, 1.2, -0.8, -2.5, 0.9, -0.6, -2.7, -0.9, -1.3, -2.0, -1.2, 0.0},
 126     {0.0, -0.4, 0.3, 2.4, -2.6, 0.6, 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8},
 127 };
 128
 129 static double gtnm_wmm2005[13][13]=
 130 {
 131     {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},
 132     {8.0, 10.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 133     {-15.1, -7.8, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 134     {0.4, -2.6, -1.2, -6.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 135     {-2.5, 2.8, -7.0, 6.2, -3.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 136     {-2.8, 0.7, -3.2, -1.1, 0.1, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 137     {-0.7, 0.4, -0.3, 2.3, -2.1, -0.6, 1.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 138     {0.2, -0.1, -0.3, 1.1, 0.6, 0.5, -0.4, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0},
 139     {0.1, 0.3, -0.4, 0.3, -0.3, 0.2, 0.4, -0.7, 0.4, 0.0, 0.0, 0.0, 0.0},
 140     {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},
 141     {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},
 142     {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},
 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 };
 145
 146 static double htnm_wmm2005[13][13]=
 147 {
 148     {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},
 149     {0.0, -20.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 150     {0.0, -23.2, -14.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 151     {0.0, 5.0, -7.0, -0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 152     {0.0, 2.2, 1.6, 5.8, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 153     {0.0, 0.0, 1.7, 2.1, 4.8, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 154     {0.0, -0.6, -1.9, -0.4, -0.5, -0.3, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 155     {0.0, 0.6, 0.4, 0.2, 0.3, -0.8, -0.2, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
 156     {0.0, -0.2, 0.1, 0.3, 0.4, 0.1, -0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0},
 157     {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},
 158     {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},
 159     {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},
 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 };
 162
 163 static const int nmax = 12;
 164
 165 static double P[13][13];
 166 static double DP[13][13];
 167 static double gnm[13][13];
 168 static double hnm[13][13];
 169 static double sm[13];
 170 static double cm[13];
 171
 172 static double root[13];
 173 static double roots[13][13][2];
 174
 175 /* Convert date to Julian day    1950-2049 */
 176 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
 177 {
 178     unsigned long jd;
 179
 180     yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
 181     jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
 182     jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
 183     jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
 184
 185     /* printf("julian date = %d\n", jd ); */
 186     return jd;
 187 }
 188
 189
 190 /*
 191  * return variation (in radians) given geodetic latitude (radians),
 192  * longitude(radians), height (km) and (Julian) date
 193  * N and E lat and long are positive, S and W negative
 194 */
 195
 196 double calc_magvar( double lat, double lon, double h, long dat, double* field )
 197 {
 198     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 199     int n,m;
 200     /* reference date for current model is 1 januari 2005 */
 201     long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
 202
 203     double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
 204     double sinpsi, cospsi, inv_s;
 205
 206     static int been_here = 0;
 207
 208     double sinlat = sin(lat);
 209     double coslat = cos(lat);
 210
 211     /* convert to geocentric coords: */
 212     // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
 213     sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
 214     /* sr is effective radius */
 215     theta = atan2(coslat * (h*sr + a*a),
 216                   sinlat * (h*sr + b*b));
 217     /* theta is geocentric co-latitude */
 218
 219     r = h*h + 2.0*h * sr +
 220         (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
 221         (a*a - (a*a - b*b) * sinlat*sinlat );
 222
 223     r = sqrt(r);
 224
 225     /* r is geocentric radial distance */
 226     c = cos(theta);
 227     s = sin(theta);
 228     /* protect against zero divide at geographic poles */
 229     inv_s =  1.0 / (s + (s == 0.)*1.0e-8);
 230
 231     /* zero out arrays */
 232     for ( n = 0; n <= nmax; n++ ) {
 233         for ( m = 0; m <= n; m++ ) {
 234             P[n][m] = 0;
 235             DP[n][m] = 0;
 236         }
 237     }
 238
 239     /* diagonal elements */
 240     P[0][0] = 1;
 241     P[1][1] = s;
 242     DP[0][0] = 0;
 243     DP[1][1] = c;
 244     P[1][0] = c ;
 245     DP[1][0] = -s;
 246
 247     // these values will not change for subsequent function calls
 248     if( !been_here ) {
 249         for ( n = 2; n <= nmax; n++ ) {
 250             root[n] = sqrt((2.0*n-1) / (2.0*n));
 251         }
 252
 253         for ( m = 0; m <= nmax; m++ ) {
 254             double mm = m*m;
 255             for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 256                 roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
 257                 roots[m][n][1] = 1.0 / sqrt( n*n - mm);
 258             }
 259         }
 260         been_here = 1;
 261     }
 262
 263     for ( n=2; n <= nmax; n++ ) {
 264         // double root = sqrt((2.0*n-1) / (2.0*n));
 265         P[n][n] = P[n-1][n-1] * s * root[n];
 266         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 267             root[n];
 268     }
 269
 270     /* lower triangle */
 271     for ( m = 0; m <= nmax; m++ ) {
 272         // double mm = m*m;
 273         for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 274             // double root1 = sqrt((n-1)*(n-1) - mm);
 275             // double root2 = 1.0 / sqrt( n*n - mm);
 276             P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
 277                        P[n-2][m] * roots[m][n][0]) *
 278                 roots[m][n][1];
 279
 280             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 281                         (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) *
 282                 roots[m][n][1];
 283         }
 284     }
 285
 286     /* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
 287        achieved by adjusting the coefficients at time t0 for linear secular variation */
 288     /* WMM2005 */
 289     yearfrac = (dat - date0_wmm2005) / 365.25;
 290     for ( n = 1; n <= nmax; n++ ) {
 291         for ( m = 0; m <= nmax; m++ ) {
 292             gnm[n][m] = gnm_wmm2005[n][m] + yearfrac * gtnm_wmm2005[n][m];
 293             hnm[n][m] = hnm_wmm2005[n][m] + yearfrac * htnm_wmm2005[n][m];
 294         }
 295     }
 296
 297     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 298     for ( m = 0; m <= nmax; m++ ) {
 299         sm[m] = sin(m * lon);
 300         cm[m] = cos(m * lon);
 301     }
 302
 303     /* compute B fields */
 304     B_r = 0.0;
 305     B_theta = 0.0;
 306     B_phi = 0.0;
 307     fn_0 = r_0/r;
 308     fn = fn_0 * fn_0;
 309
 310     for ( n = 1; n <= nmax; n++ ) {
 311         double c1_n=0;
 312         double c2_n=0;
 313         double c3_n=0;
 314         for ( m = 0; m <= n; m++ ) {
 315             double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
 316             c1_n=c1_n + tmp * P[n][m];
 317             c2_n=c2_n + tmp * DP[n][m];
 318             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 319         }
 320         // fn=pow(r_0/r,n+2.0);
 321         fn *= fn_0;
 322         B_r = B_r + (n + 1) * c1_n * fn;
 323         B_theta = B_theta - c2_n * fn;
 324         B_phi = B_phi + c3_n * fn * inv_s;
 325     }
 326
 327     /* Find geodetic field components: */
 328     psi = theta - ((pi / 2.0) - lat);
 329     sinpsi = sin(psi);
 330     cospsi = cos(psi);
 331     X = -B_theta * cospsi - B_r * sinpsi;
 332     Y = B_phi;
 333     Z = B_theta * sinpsi - B_r * cospsi;
 334
 335     field[0]=B_r;
 336     field[1]=B_theta;
 337     field[2]=B_phi;
 338     field[3]=X;
 339     field[4]=Y;
 340     field[5]=Z;   /* output fields */
 341
 342     /* find variation in radians */
 343     /* return zero variation at magnetic pole X=Y=0. */
 344     /* E is positive */
 345     return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
 346 }
 347
 348
 349 #ifdef TEST_NHV_HACKS
 350 double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
 351 {
 352     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 353     int n,m;
 354     /* reference dates */
 355     long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
 356
 357     double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
 358
 359     /* convert to geocentric coords: */
 360     sr = sqrt(pow(a*cos(lat),2.0)+pow(b*sin(lat),2.0));
 361     /* sr is effective radius */
 362     theta = atan2(cos(lat) * (h * sr + a * a),
 363                   sin(lat) * (h * sr + b * b));
 364     /* theta is geocentric co-latitude */
 365
 366     r = h * h + 2.0*h * sr +
 367         (pow(a,4.0) - (pow(a,4.0) - pow(b,4.0)) * pow(sin(lat),2.0)) /
 368         (a * a - (a * a - b * b) * pow(sin(lat),2.0));
 369
 370     r = sqrt(r);
 371
 372     /* r is geocentric radial distance */
 373     c = cos(theta);
 374     s = sin(theta);
 375
 376     /* zero out arrays */
 377     for ( n = 0; n <= nmax; n++ ) {
 378         for ( m = 0; m <= n; m++ ) {
 379             P[n][m] = 0;
 380             DP[n][m] = 0;
 381         }
 382     }
 383
 384     /* diagonal elements */
 385     P[0][0] = 1;
 386     P[1][1] = s;
 387     DP[0][0] = 0;
 388     DP[1][1] = c;
 389     P[1][0] = c ;
 390     DP[1][0] = -s;
 391
 392     for ( n = 2; n <= nmax; n++ ) {
 393         P[n][n] = P[n-1][n-1] * s * sqrt((2.0*n-1) / (2.0*n));
 394         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 395             sqrt((2.0*n-1) / (2.0*n));
 396     }
 397
 398     /* lower triangle */
 399     for ( m = 0; m <= nmax; m++ ) {
 400         for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
 401             P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
 402                        sqrt(1.0*(n-1)*(n-1) - m * m)) /
 403                 sqrt(1.0* n * n - m * m);
 404             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 405                         (2.0*n-1) - DP[n-2][m] *
 406                         sqrt(1.0*(n-1) * (n-1) - m * m)) /
 407                 sqrt(1.0* n * n - m * m);
 408         }
 409     }
 410
 411     /* compute gnm, hnm at dat */
 412     /* WMM2005 */
 413     yearfrac = (dat - date0_wmm2005) / 365.25;
 414     for ( n = 1; n <= nmax; n++ ) {
 415         for ( m = 0; m <= nmax; m++ ) {
 416             gnm[n][m] = gnm_wmm2005[n][m] + yearfrac * gtnm_wmm2005[n][m];
 417             hnm[n][m] = hnm_wmm2005[n][m] + yearfrac * htnm_wmm2005[n][m];
 418         }
 419     }
 420
 421     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 422     for ( m = 0; m <= nmax; m++ ) {
 423         sm[m] = sin(m * lon);
 424         cm[m] = cos(m * lon);
 425     }
 426
 427     /* compute B fields */
 428     B_r = 0.0;
 429     B_theta = 0.0;
 430     B_phi = 0.0;
 431
 432     for ( n = 1; n <= nmax; n++ ) {
 433         double c1_n=0;
 434         double c2_n=0;
 435         double c3_n=0;
 436         for ( m = 0; m <= n; m++ ) {
 437             c1_n=c1_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * P[n][m];
 438             c2_n=c2_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * DP[n][m];
 439             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 440         }
 441         fn=pow(r_0/r,n+2.0);
 442         B_r = B_r + (n + 1) * c1_n * fn;
 443         B_theta = B_theta - c2_n * fn;
 444         B_phi = B_phi + c3_n * fn / s;
 445     }
 446
 447     /* Find geodetic field components: */
 448     psi = theta - (pi / 2.0 - lat);
 449     X = -B_theta * cos(psi) - B_r * sin(psi);
 450     Y = B_phi;
 451     Z = B_theta * sin(psi) - B_r * cos(psi);
 452
 453     field[0]=B_r;
 454     field[1]=B_theta;
 455     field[2]=B_phi;
 456     field[3]=X;
 457     field[4]=Y;
 458     field[5]=Z;   /* output fields */
 459
 460     /* find variation, leave in radians! */
 461     return atan2(Y, X);  /* E is positive */
 462 }
 463 #endif // TEST_NHV_HACKS