simgear/magvar/magvar.cxx

   1 // magvar.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 //  The routine uses a spherical harmonic expansion of the magnetic
  28 // potential up to twelfth order, together with its time variation, as
  29 // described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
  30 // Academic Press (London 1987)". The program first converts geodetic
  31 // coordinates (lat/long on elliptic earth and altitude) to spherical
  32 // geocentric (spherical lat/long and radius) coordinates. Using this,
  33 // the spherical (B_r, B_theta, B_phi) magnetic field components are
  34 // computed from the model. These are finally referred to surface (X, Y,
  35 // Z) coordinates.
  36 //
  37 //   Fields are accurate to better than 200nT, variation and dip to
  38 // better than 0.5 degrees, with the exception of the declination near
  39 // the magnetic poles (where it is ill-defined) where the error may reach
  40 // 4 degrees or more.
  41 //
  42 //   Variation is undefined at both the geographic and
  43 // magnetic poles, even though the field itself is well-behaved. To
  44 // avoid the routine blowing up, latitude entries corresponding to
  45 // the geographic poles are slightly offset. At the magnetic poles,
  46 // the routine returns zero variation.
  47
  48
  49 //
  50 // This library is free software; you can redistribute it and/or
  51 // modify it under the terms of the GNU Library General Public
  52 // License as published by the Free Software Foundation; either
  53 // version 2 of the License, or (at your option) any later version.
  54 //
  55 // This library is distributed in the hope that it will be useful,
  56 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  57 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  58 // Library General Public License for more details.
  59 //
  60 // You should have received a copy of the GNU Library General Public
  61 // License along with this library; if not, write to the
  62 // Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  63 // Boston, MA  02111-1307, USA.
  64 //
  65 // $Id$
  66
  67
  68 #include <stdio.h>
  69 #include <stdlib.h>
  70 #include <math.h>
  71
  72 #include "magvar.hxx"
  73
  74
  75 #define max(a,b)        (((a) > (b)) ? (a) : (b))
  76
  77 static const double pi = 3.14159265358979;
  78 static const double a = 6378.16;        /* major radius (km) IAU66 ellipsoid */
  79 static const double f = 1.0 / 298.25;   /* inverse flattening IAU66 ellipsoid */
  80 static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
  81 /* minor radius b=a*(1-f) */
  82 static const double r_0 = 6371.2;       /* "mean radius" for spherical harmonic expansion */
  83
  84 static double gnm_wmm2000[13][13] =
  85 {
  86     {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},
  87     {-29616.0, -1722.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  88     {-2266.7, 3070.2, 1677.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  89     {1322.4, -2291.5, 1255.9, 724.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  90     {932.1, 786.3, 250.6, -401.5, 106.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  91     {-211.9, 351.6, 220.8, -134.5, -168.8, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  92     {73.8, 68.2, 74.1, -163.5, -3.8, 17.1, -85.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  93     {77.4, -73.9, 2.2, 35.7, 7.3, 5.2, 8.4, -1.5, 0.0, 0.0, 0.0, 0.0, 0.0},
  94     {23.3, 7.3, -8.5, -6.6, -16.9, 8.6, 4.9, -7.8, -7.6, 0.0, 0.0, 0.0, 0.0},
  95     {5.7, 8.5, 2.0, -9.8, 7.6, -7.0, -2.0, 9.2, -2.2, -6.6, 0.0, 0.0, 0.0},
  96     {-2.2, -5.7, 1.6, -3.7, -0.6, 4.1, 2.2, 2.2, 4.6, 2.3, 0.1, 0.0, 0.0},
  97     {3.3, -1.1, -2.4, 2.6, -1.3, -1.7, -0.6, 0.4, 0.7, -0.3, 2.3, 4.2, 0.0},
  98     {-1.5, -0.2, -0.3, 0.5, 0.2, 0.9, -1.4, 0.6, -0.6, -1.0, -0.3, 0.3, 0.4},
  99 };
 100
 101 static double hnm_wmm2000[13][13]=
 102 {
 103     {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},
 104     {0.0, 5194.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 105     {0.0, -2484.8, -467.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 106     {0.0, -224.7, 293.0, -486.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 107     {0.0, 273.3, -227.9, 120.9, -302.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 108     {0.0, 42.0, 173.8, -135.0, -38.6, 105.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 109     {0.0, -17.4, 61.2, 63.2, -62.9, 0.2, 43.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 110     {0.0, -62.3, -24.5, 8.9, 23.4, 15.0, -27.6, -7.8, 0.0, 0.0, 0.0, 0.0, 0.0},
 111     {0.0, 12.4, -20.8, 8.4, -21.2, 15.5, 9.1, -15.5, -5.4, 0.0, 0.0, 0.0, 0.0},
 112     {0.0, -20.4, 13.9, 12.0, -6.2, -8.6, 9.4, 5.0, -8.4, 3.2, 0.0, 0.0, 0.0},
 113     {0.0, 0.9, -0.7, 3.9, 4.8, -5.3, -1.0, -2.4, 1.3, -2.3, -6.4, 0.0, 0.0},
 114     {0.0, -1.5, 0.7, -1.1, -2.3, 1.3, -0.6, -2.8, -1.6, -0.1, -1.9, 1.4, 0.0},
 115     {0.0, -1.0, 0.7, 2.2, -2.5, -0.2, 0.0, -0.2, 0.0, 0.2, -0.9, -0.2, 1.0},
 116 };
 117
 118 static double gtnm_wmm2000[13][13]=
 119 {
 120     {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},
 121     {14.7, 11.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 122     {-13.6, -0.7, -1.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 123     {0.3, -4.3, 0.9, -8.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 124     {-1.6, 0.9, -7.6, 2.2, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 125     {-0.9, -0.2, -2.5, -2.7, -0.9, 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 126     {1.2, 0.2, 1.7, 1.6, -0.1, -0.3, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 127     {-0.4, -0.8, -0.2, 1.1, 0.4, 0.0, -0.2, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
 128     {-0.3, 0.6, -0.8, 0.3, -0.2, 0.5, 0.0, -0.6, 0.1, 0.0, 0.0, 0.0, 0.0},
 129     {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},
 130     {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},
 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     {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},
 133 };
 134
 135 static double htnm_wmm2000[13][13]=
 136 {
 137     {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},
 138     {0.0, -20.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 139     {0.0, -21.5, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 140     {0.0, 6.4, -1.3, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 141     {0.0, 2.3, 0.7, 3.7, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 142     {0.0, 0.0, 2.1, 2.3, 3.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 143     {0.0, -0.3, -1.7, -0.9, -1.0, -0.1, 1.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 144     {0.0, 1.4, 0.2, 0.7, 0.4, -0.3, -0.8, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
 145     {0.0, -0.5, 0.1, -0.2, 0.0, 0.1, -0.1, 0.3, 0.2, 0.0, 0.0, 0.0, 0.0},
 146     {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},
 147     {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},
 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 150 };
 151
 152 static const int nmax = 12;
 153
 154 static double P[13][13];
 155 static double DP[13][13];
 156 static double gnm[13][13];
 157 static double hnm[13][13];
 158 static double sm[13];
 159 static double cm[13];
 160
 161 static double root[13];
 162 static double roots[13][13][2];
 163
 164 /* Convert date to Julian day    1950-2049 */
 165 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
 166 {
 167     unsigned long jd;
 168
 169     yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
 170     jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
 171     jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
 172     jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
 173
 174     /* printf("julian date = %d\n", jd ); */
 175     return jd;
 176 }
 177
 178
 179 /* Convert degrees to radians */
 180 double deg_to_rad( double deg )
 181 {
 182     return deg*pi/180.;
 183 }
 184
 185
 186 /* Convert radians to degrees */
 187 double rad_to_deg( double rad )
 188 {
 189     return rad*180./pi;
 190 }
 191
 192
 193 /*
 194  * return variation (in radians) given geodetic latitude (radians),
 195  * longitude(radians), height (km) and (Julian) date
 196  * N and E lat and long are positive, S and W negative
 197 */
 198
 199 double SGMagVar( double lat, double lon, double h, long dat, double* field )
 200 {
 201     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 202     int n,m;
 203     /* reference dates */
 204     long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
 205
 206     double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
 207     double sinpsi, cospsi, inv_s;
 208
 209     static int been_here = 0;
 210
 211     double sinlat = sin(lat);
 212     double coslat = cos(lat);
 213
 214     /* convert to geocentric coords: */
 215     // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
 216     sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
 217     /* sr is effective radius */
 218     theta = atan2(coslat * (h*sr + a*a),
 219                   sinlat * (h*sr + b*b));
 220     /* theta is geocentric co-latitude */
 221
 222     r = h*h + 2.0*h * sr +
 223         (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
 224         (a*a - (a*a - b*b) * sinlat*sinlat );
 225
 226     r = sqrt(r);
 227
 228     /* r is geocentric radial distance */
 229     c = cos(theta);
 230     s = sin(theta);
 231     /* protect against zero divide at geographic poles */
 232     inv_s =  1.0 / (s + (s == 0.)*1.0e-8);
 233
 234     /* zero out arrays */
 235     for ( n = 0; n <= nmax; n++ ) {
 236         for ( m = 0; m <= n; m++ ) {
 237             P[n][m] = 0;
 238             DP[n][m] = 0;
 239         }
 240     }
 241
 242     /* diagonal elements */
 243     P[0][0] = 1;
 244     P[1][1] = s;
 245     DP[0][0] = 0;
 246     DP[1][1] = c;
 247     P[1][0] = c ;
 248     DP[1][0] = -s;
 249
 250     // these values will not change for subsequent function calls
 251     if( !been_here ) {
 252         for ( n = 2; n <= nmax; n++ ) {
 253             root[n] = sqrt((2.0*n-1) / (2.0*n));
 254         }
 255
 256         for ( m = 0; m <= nmax; m++ ) {
 257             double mm = m*m;
 258             for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 259                 roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
 260                 roots[m][n][1] = 1.0 / sqrt( n*n - mm);
 261             }
 262         }
 263         been_here = 1;
 264     }
 265
 266     for ( n=2; n <= nmax; n++ ) {
 267         // double root = sqrt((2.0*n-1) / (2.0*n));
 268         P[n][n] = P[n-1][n-1] * s * root[n];
 269         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 270             root[n];
 271     }
 272
 273     /* lower triangle */
 274     for ( m = 0; m <= nmax; m++ ) {
 275         // double mm = m*m;
 276         for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 277             // double root1 = sqrt((n-1)*(n-1) - mm);
 278             // double root2 = 1.0 / sqrt( n*n - mm);
 279             P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
 280                        P[n-2][m] * roots[m][n][0]) *
 281                 roots[m][n][1];
 282
 283             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 284                         (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) *
 285                 roots[m][n][1];
 286         }
 287     }
 288
 289     /* compute gnm, hnm at dat */
 290     /* WMM2000 */
 291     yearfrac = (dat - date0_wmm2000) / 365.25;
 292     for ( n = 1; n <= nmax; n++ ) {
 293         for ( m = 0; m <= nmax; m++ ) {
 294             gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
 295             hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
 296         }
 297     }
 298
 299     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 300     for ( m = 0; m <= nmax; m++ ) {
 301         sm[m] = sin(m * lon);
 302         cm[m] = cos(m * lon);
 303     }
 304
 305     /* compute B fields */
 306     B_r = 0.0;
 307     B_theta = 0.0;
 308     B_phi = 0.0;
 309     fn_0 = r_0/r;
 310     fn = fn_0 * fn_0;
 311
 312     for ( n = 1; n <= nmax; n++ ) {
 313         double c1_n=0;
 314         double c2_n=0;
 315         double c3_n=0;
 316         for ( m = 0; m <= n; m++ ) {
 317             double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
 318             c1_n=c1_n + tmp * P[n][m];
 319             c2_n=c2_n + tmp * DP[n][m];
 320             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 321         }
 322         // fn=pow(r_0/r,n+2.0);
 323         fn *= fn_0;
 324         B_r = B_r + (n + 1) * c1_n * fn;
 325         B_theta = B_theta - c2_n * fn;
 326         B_phi = B_phi + c3_n * fn * inv_s;
 327     }
 328
 329     /* Find geodetic field components: */
 330     psi = theta - ((pi / 2.0) - lat);
 331     sinpsi = sin(psi);
 332     cospsi = cos(psi);
 333     X = -B_theta * cospsi - B_r * sinpsi;
 334     Y = B_phi;
 335     Z = B_theta * sinpsi - B_r * cospsi;
 336
 337     field[0]=B_r;
 338     field[1]=B_theta;
 339     field[2]=B_phi;
 340     field[3]=X;
 341     field[4]=Y;
 342     field[5]=Z;   /* output fields */
 343
 344     /* find variation in radians */
 345     /* return zero variation at magnetic pole X=Y=0. */
 346     /* E is positive */
 347     return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
 348 }
 349
 350
 351 #ifdef TEST_NHV_HACKS
 352 double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
 353 {
 354     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 355     int n,m;
 356     /* reference dates */
 357     long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
 358
 359     double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
 360
 361     /* convert to geocentric coords: */
 362     sr = sqrt(pow(a*cos(lat),2.0)+pow(b*sin(lat),2.0));
 363     /* sr is effective radius */
 364     theta = atan2(cos(lat) * (h * sr + a * a),
 365                   sin(lat) * (h * sr + b * b));
 366     /* theta is geocentric co-latitude */
 367
 368     r = h * h + 2.0*h * sr +
 369         (pow(a,4.0) - (pow(a,4.0) - pow(b,4.0)) * pow(sin(lat),2.0)) /
 370         (a * a - (a * a - b * b) * pow(sin(lat),2.0));
 371
 372     r = sqrt(r);
 373
 374     /* r is geocentric radial distance */
 375     c = cos(theta);
 376     s = sin(theta);
 377
 378     /* zero out arrays */
 379     for ( n = 0; n <= nmax; n++ ) {
 380         for ( m = 0; m <= n; m++ ) {
 381             P[n][m] = 0;
 382             DP[n][m] = 0;
 383         }
 384     }
 385
 386     /* diagonal elements */
 387     P[0][0] = 1;
 388     P[1][1] = s;
 389     DP[0][0] = 0;
 390     DP[1][1] = c;
 391     P[1][0] = c ;
 392     DP[1][0] = -s;
 393
 394     for ( n = 2; n <= nmax; n++ ) {
 395         P[n][n] = P[n-1][n-1] * s * sqrt((2.0*n-1) / (2.0*n));
 396         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 397             sqrt((2.0*n-1) / (2.0*n));
 398     }
 399
 400     /* lower triangle */
 401     for ( m = 0; m <= nmax; m++ ) {
 402         for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 403             P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
 404                        sqrt(1.0*(n-1)*(n-1) - m * m)) /
 405                 sqrt(1.0* n * n - m * m);
 406             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 407                         (2.0*n-1) - DP[n-2][m] *
 408                         sqrt(1.0*(n-1) * (n-1) - m * m)) /
 409                 sqrt(1.0* n * n - m * m);
 410         }
 411     }
 412
 413     /* compute gnm, hnm at dat */
 414     /* WMM2000 */
 415     yearfrac = (dat - date0_wmm2000) / 365.25;
 416     for ( n = 1; n <= nmax; n++ ) {
 417         for ( m = 0; m <= nmax; m++ ) {
 418             gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
 419             hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
 420         }
 421     }
 422
 423     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 424     for ( m = 0; m <= nmax; m++ ) {
 425         sm[m] = sin(m * lon);
 426         cm[m] = cos(m * lon);
 427     }
 428
 429     /* compute B fields */
 430     B_r = 0.0;
 431     B_theta = 0.0;
 432     B_phi = 0.0;
 433
 434     for ( n = 1; n <= nmax; n++ ) {
 435         double c1_n=0;
 436         double c2_n=0;
 437         double c3_n=0;
 438         for ( m = 0; m <= n; m++ ) {
 439             c1_n=c1_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * P[n][m];
 440             c2_n=c2_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * DP[n][m];
 441             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 442         }
 443         fn=pow(r_0/r,n+2.0);
 444         B_r = B_r + (n + 1) * c1_n * fn;
 445         B_theta = B_theta - c2_n * fn;
 446         B_phi = B_phi + c3_n * fn / s;
 447     }
 448
 449     /* Find geodetic field components: */
 450     psi = theta - (pi / 2.0 - lat);
 451     X = -B_theta * cos(psi) - B_r * sin(psi);
 452     Y = B_phi;
 453     Z = B_theta * sin(psi) - B_r * cos(psi);
 454
 455     field[0]=B_r;
 456     field[1]=B_theta;
 457     field[2]=B_phi;
 458     field[3]=X;
 459     field[4]=Y;
 460     field[5]=Z;   /* output fields */
 461
 462     /* find variation, leave in radians! */
 463     return atan2(Y, X);  /* E is positive */
 464 }
 465 #endif // TEST_NHV_HACKS