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 // 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,
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
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.
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.
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.
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.
72 #include <simgear/constants.h>
74 #include "coremag.hxx"
77 #define max(a,b) (((a) > (b)) ? (a) : (b))
79 static const double pi = 3.14159265358979;
80 static const double a = 6378.16; /* major radius (km) IAU66 ellipsoid */
81 static const double f = 1.0 / 298.25; /* inverse flattening IAU66 ellipsoid */
82 static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
83 /* minor radius b=a*(1-f) */
84 static const double r_0 = 6371.2; /* "mean radius" for spherical harmonic expansion */
86 static double gnm_wmm2000[13][13] =
88 {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},
89 {-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},
90 {-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},
91 {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},
92 {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},
93 {-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},
94 {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},
95 {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},
96 {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},
97 {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},
98 {-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},
99 {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},
100 {-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},
103 static double hnm_wmm2000[13][13]=
105 {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},
106 {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},
107 {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},
108 {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},
109 {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},
110 {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},
111 {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},
112 {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},
113 {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},
114 {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},
115 {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},
116 {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},
117 {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},
120 static double gtnm_wmm2000[13][13]=
122 {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},
123 {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},
124 {-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},
125 {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},
126 {-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},
127 {-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},
128 {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},
129 {-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},
130 {-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},
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 {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},
134 {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},
137 static double htnm_wmm2000[13][13]=
139 {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},
140 {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},
141 {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},
142 {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},
143 {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},
144 {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},
145 {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},
146 {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},
147 {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},
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 {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},
151 {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},
154 static const int nmax = 12;
156 static double P[13][13];
157 static double DP[13][13];
158 static double gnm[13][13];
159 static double hnm[13][13];
160 static double sm[13];
161 static double cm[13];
163 static double root[13];
164 static double roots[13][13][2];
166 /* Convert date to Julian day 1950-2049 */
167 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
171 yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
172 jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
173 jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
174 jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
176 /* printf("julian date = %d\n", jd ); */
182 * return variation (in radians) given geodetic latitude (radians),
183 * longitude(radians), height (km) and (Julian) date
184 * N and E lat and long are positive, S and W negative
187 double calc_magvar( double lat, double lon, double h, long dat, double* field )
189 /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
191 /* reference dates */
192 long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
194 double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
195 double sinpsi, cospsi, inv_s;
197 static int been_here = 0;
199 double sinlat = sin(lat);
200 double coslat = cos(lat);
202 /* convert to geocentric coords: */
203 // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
204 sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
205 /* sr is effective radius */
206 theta = atan2(coslat * (h*sr + a*a),
207 sinlat * (h*sr + b*b));
208 /* theta is geocentric co-latitude */
210 r = h*h + 2.0*h * sr +
211 (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
212 (a*a - (a*a - b*b) * sinlat*sinlat );
216 /* r is geocentric radial distance */
219 /* protect against zero divide at geographic poles */
220 inv_s = 1.0 / (s + (s == 0.)*1.0e-8);
222 /* zero out arrays */
223 for ( n = 0; n <= nmax; n++ ) {
224 for ( m = 0; m <= n; m++ ) {
230 /* diagonal elements */
238 // these values will not change for subsequent function calls
240 for ( n = 2; n <= nmax; n++ ) {
241 root[n] = sqrt((2.0*n-1) / (2.0*n));
244 for ( m = 0; m <= nmax; m++ ) {
246 for ( n = max(m + 1, 2); n <= nmax; n++ ) {
247 roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
248 roots[m][n][1] = 1.0 / sqrt( n*n - mm);
254 for ( n=2; n <= nmax; n++ ) {
255 // double root = sqrt((2.0*n-1) / (2.0*n));
256 P[n][n] = P[n-1][n-1] * s * root[n];
257 DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
262 for ( m = 0; m <= nmax; m++ ) {
264 for ( n = max(m + 1, 2); n <= nmax; n++ ) {
265 // double root1 = sqrt((n-1)*(n-1) - mm);
266 // double root2 = 1.0 / sqrt( n*n - mm);
267 P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
268 P[n-2][m] * roots[m][n][0]) *
271 DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
272 (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) *
277 /* compute gnm, hnm at dat */
279 yearfrac = (dat - date0_wmm2000) / 365.25;
280 for ( n = 1; n <= nmax; n++ ) {
281 for ( m = 0; m <= nmax; m++ ) {
282 gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
283 hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
287 /* compute sm (sin(m lon) and cm (cos(m lon)) */
288 for ( m = 0; m <= nmax; m++ ) {
289 sm[m] = sin(m * lon);
290 cm[m] = cos(m * lon);
293 /* compute B fields */
300 for ( n = 1; n <= nmax; n++ ) {
304 for ( m = 0; m <= n; m++ ) {
305 double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
306 c1_n=c1_n + tmp * P[n][m];
307 c2_n=c2_n + tmp * DP[n][m];
308 c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
310 // fn=pow(r_0/r,n+2.0);
312 B_r = B_r + (n + 1) * c1_n * fn;
313 B_theta = B_theta - c2_n * fn;
314 B_phi = B_phi + c3_n * fn * inv_s;
317 /* Find geodetic field components: */
318 psi = theta - ((pi / 2.0) - lat);
321 X = -B_theta * cospsi - B_r * sinpsi;
323 Z = B_theta * sinpsi - B_r * cospsi;
330 field[5]=Z; /* output fields */
332 /* find variation in radians */
333 /* return zero variation at magnetic pole X=Y=0. */
335 return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
339 #ifdef TEST_NHV_HACKS
340 double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
342 /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
344 /* reference dates */
345 long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
347 double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
349 /* convert to geocentric coords: */
350 sr = sqrt(pow(a*cos(lat),2.0)+pow(b*sin(lat),2.0));
351 /* sr is effective radius */
352 theta = atan2(cos(lat) * (h * sr + a * a),
353 sin(lat) * (h * sr + b * b));
354 /* theta is geocentric co-latitude */
356 r = h * h + 2.0*h * sr +
357 (pow(a,4.0) - (pow(a,4.0) - pow(b,4.0)) * pow(sin(lat),2.0)) /
358 (a * a - (a * a - b * b) * pow(sin(lat),2.0));
362 /* r is geocentric radial distance */
366 /* zero out arrays */
367 for ( n = 0; n <= nmax; n++ ) {
368 for ( m = 0; m <= n; m++ ) {
374 /* diagonal elements */
382 for ( n = 2; n <= nmax; n++ ) {
383 P[n][n] = P[n-1][n-1] * s * sqrt((2.0*n-1) / (2.0*n));
384 DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
385 sqrt((2.0*n-1) / (2.0*n));
389 for ( m = 0; m <= nmax; m++ ) {
390 for ( n = max(m + 1, 2); n <= nmax; n++ ) {
391 P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
392 sqrt(1.0*(n-1)*(n-1) - m * m)) /
393 sqrt(1.0* n * n - m * m);
394 DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
395 (2.0*n-1) - DP[n-2][m] *
396 sqrt(1.0*(n-1) * (n-1) - m * m)) /
397 sqrt(1.0* n * n - m * m);
401 /* compute gnm, hnm at dat */
403 yearfrac = (dat - date0_wmm2000) / 365.25;
404 for ( n = 1; n <= nmax; n++ ) {
405 for ( m = 0; m <= nmax; m++ ) {
406 gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
407 hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
411 /* compute sm (sin(m lon) and cm (cos(m lon)) */
412 for ( m = 0; m <= nmax; m++ ) {
413 sm[m] = sin(m * lon);
414 cm[m] = cos(m * lon);
417 /* compute B fields */
422 for ( n = 1; n <= nmax; n++ ) {
426 for ( m = 0; m <= n; m++ ) {
427 c1_n=c1_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * P[n][m];
428 c2_n=c2_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * DP[n][m];
429 c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
432 B_r = B_r + (n + 1) * c1_n * fn;
433 B_theta = B_theta - c2_n * fn;
434 B_phi = B_phi + c3_n * fn / s;
437 /* Find geodetic field components: */
438 psi = theta - (pi / 2.0 - lat);
439 X = -B_theta * cos(psi) - B_r * sin(psi);
441 Z = B_theta * sin(psi) - B_r * cos(psi);
448 field[5]=Z; /* output fields */
450 /* find variation, leave in radians! */
451 return atan2(Y, X); /* E is positive */
453 #endif // TEST_NHV_HACKS