1 // coremag.cxx -- compute local magnetic variation given position,
4 // This is an implementation of the NIMA (formerly DMA) WMM2000
6 // http://www.nima.mil/GandG/ngdc-wmm2000.html
8 // Copyright (C) 2000 Edward A Williams <Ed_Williams@compuserve.com>
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
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
24 // Put in some bullet-proofing to handle magnetic and geographic poles.
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.
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.
37 // 25/10/2006 Wim Van Hoydonck -- wim.van.hoydonck@gmail.com
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,
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
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.
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.
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.
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.
84 #include <simgear/constants.h>
85 #include <simgear/sg_inlines.h>
87 #include "coremag.hxx"
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 */
95 static double gnm_wmm2005[13][13] =
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},
112 static double hnm_wmm2005[13][13]=
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},
129 static double gtnm_wmm2005[13][13]=
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},
146 static double htnm_wmm2005[13][13]=
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},
163 static const int nmax = 12;
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];
172 static double root[13];
173 static double roots[13][13][2];
175 /* Convert date to Julian day 1950-2049 */
176 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
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;
185 /* printf("julian date = %d\n", jd ); */
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
196 double calc_magvar( double lat, double lon, double h, long dat, double* field )
198 /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
200 /* reference date for current model is 1 januari 2005 */
201 long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
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;
206 static int been_here = 0;
208 double sinlat = sin(lat);
209 double coslat = cos(lat);
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 */
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 );
225 /* r is geocentric radial distance */
228 /* protect against zero divide at geographic poles */
229 inv_s = 1.0 / (s + (s == 0.)*1.0e-8);
231 /* zero out arrays */
232 for ( n = 0; n <= nmax; n++ ) {
233 for ( m = 0; m <= n; m++ ) {
239 /* diagonal elements */
247 // these values will not change for subsequent function calls
249 for ( n = 2; n <= nmax; n++ ) {
250 root[n] = sqrt((2.0*n-1) / (2.0*n));
253 for ( m = 0; m <= nmax; 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);
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) *
271 for ( m = 0; m <= nmax; 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]) *
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]) *
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 */
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];
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);
303 /* compute B fields */
310 for ( n = 1; n <= nmax; n++ ) {
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];
320 // fn=pow(r_0/r,n+2.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;
327 /* Find geodetic field components: */
328 psi = theta - ((pi / 2.0) - lat);
331 X = -B_theta * cospsi - B_r * sinpsi;
333 Z = B_theta * sinpsi - B_r * cospsi;
340 field[5]=Z; /* output fields */
342 /* find variation in radians */
343 /* return zero variation at magnetic pole X=Y=0. */
345 return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
349 #ifdef TEST_NHV_HACKS
350 double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
352 /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
354 /* reference dates */
355 long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
357 double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
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 */
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));
372 /* r is geocentric radial distance */
376 /* zero out arrays */
377 for ( n = 0; n <= nmax; n++ ) {
378 for ( m = 0; m <= n; m++ ) {
384 /* diagonal elements */
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));
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);
411 /* compute gnm, hnm at dat */
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];
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);
427 /* compute B fields */
432 for ( n = 1; n <= nmax; n++ ) {
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];
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;
447 /* Find geodetic field components: */
448 psi = theta - (pi / 2.0 - lat);
449 X = -B_theta * cos(psi) - B_r * sin(psi);
451 Z = B_theta * sin(psi) - B_r * cos(psi);
458 field[5]=Z; /* output fields */
460 /* find variation, leave in radians! */
461 return atan2(Y, X); /* E is positive */
463 #endif // TEST_NHV_HACKS