Branch data Line data Source code
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: coremag.cxx,v 1.1 2009/09/03 12:27:07 jberndt Exp $
78 : :
79 : :
80 : : #include <cstdio>
81 : : #include <cstdlib>
82 : : #include <cmath>
83 : :
84 : : #include "coremag.hxx"
85 : :
86 : : #define MAX(X, Y) X>Y?X:Y
87 : :
88 : : static const double a = 6378.137; /* semi-major axis (equatorial radius) of WGS84 ellipsoid */
89 : : static const double b = 6356.7523142; /* semi-minor axis referenced to the WGS84 ellipsoid */
90 : : static const double r_0 = 6371.2; /* standard Earth magnetic reference radius */
91 : :
92 : : static double gnm_wmm2005[13][13] =
93 : : {
94 : : {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},
95 : : {-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},
96 : : {-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},
97 : : {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},
98 : : {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},
99 : : {-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},
100 : : {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},
101 : : {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},
102 : : {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},
103 : : {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},
104 : : {-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},
105 : : {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},
106 : : {-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},
107 : : };
108 : :
109 : : static double hnm_wmm2005[13][13]=
110 : : {
111 : : {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},
112 : : {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},
113 : : {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},
114 : : {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},
115 : : {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},
116 : : {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},
117 : : {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},
118 : : {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},
119 : : {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},
120 : : {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},
121 : : {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},
122 : : {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},
123 : : {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},
124 : : };
125 : :
126 : : static double gtnm_wmm2005[13][13]=
127 : : {
128 : : {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},
129 : : {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},
130 : : {-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},
131 : : {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},
132 : : {-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},
133 : : {-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},
134 : : {-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},
135 : : {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},
136 : : {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},
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, 0.0, 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, 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, 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 : : };
142 : :
143 : : static double htnm_wmm2005[13][13]=
144 : : {
145 : : {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 : : {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},
147 : : {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},
148 : : {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},
149 : : {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},
150 : : {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},
151 : : {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},
152 : : {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},
153 : : {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},
154 : : {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},
155 : : {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},
156 : : {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},
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 : : };
159 : :
160 : : static const int nmax = 12;
161 : :
162 : : static double P[13][13];
163 : : static double DP[13][13];
164 : : static double gnm[13][13];
165 : : static double hnm[13][13];
166 : : static double sm[13];
167 : : static double cm[13];
168 : :
169 : : static double root[13];
170 : : static double roots[13][13][2];
171 : :
172 : : /* Convert date to Julian day 1950-2049 */
173 : 0 : unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
174 : : {
175 : : unsigned long jd;
176 : :
177 [ # # ]: 0 : yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
178 : 0 : jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
179 : 0 : jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
180 : 0 : jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
181 : :
182 : : /* printf("julian date = %d\n", jd ); */
183 : 0 : return jd;
184 : : }
185 : :
186 : :
187 : : /*
188 : : * return variation (in radians) given geodetic latitude (radians),
189 : : * longitude(radians), height (km) and (Julian) date
190 : : * N and E lat and long are positive, S and W negative
191 : : */
192 : :
193 : 0 : double calc_magvar( double lat, double lon, double h, long dat, double* field )
194 : : {
195 : : /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
196 : : int n,m;
197 : : /* reference date for current model is 1 januari 2005 */
198 : 0 : long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
199 : :
200 : : double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
201 : : double sinpsi, cospsi, inv_s;
202 : :
203 : : static int been_here = 0;
204 : :
205 : 0 : double sinlat = sin(lat);
206 : 0 : double coslat = cos(lat);
207 : :
208 : : /* convert to geocentric coords: */
209 : : // sr = sqrt(pow(a*coslat,2.0)+pow(b*sinlat,2.0));
210 : 0 : sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
211 : : /* sr is effective radius */
212 : : theta = atan2(coslat * (h*sr + a*a),
213 : 0 : sinlat * (h*sr + b*b));
214 : : /* theta is geocentric co-latitude */
215 : :
216 : : r = h*h + 2.0*h * sr +
217 : : (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
218 : 0 : (a*a - (a*a - b*b) * sinlat*sinlat );
219 : :
220 : 0 : r = sqrt(r);
221 : :
222 : : /* r is geocentric radial distance */
223 : 0 : c = cos(theta);
224 : 0 : s = sin(theta);
225 : : /* protect against zero divide at geographic poles */
226 [ # # ]: 0 : inv_s = 1.0 / (s + (s == 0.)*1.0e-8);
227 : :
228 : : /* zero out arrays */
229 [ # # ]: 0 : for ( n = 0; n <= nmax; n++ ) {
230 [ # # ]: 0 : for ( m = 0; m <= n; m++ ) {
231 : 0 : P[n][m] = 0;
232 : 0 : DP[n][m] = 0;
233 : : }
234 : : }
235 : :
236 : : /* diagonal elements */
237 : 0 : P[0][0] = 1;
238 : 0 : P[1][1] = s;
239 : 0 : DP[0][0] = 0;
240 : 0 : DP[1][1] = c;
241 : 0 : P[1][0] = c ;
242 : 0 : DP[1][0] = -s;
243 : :
244 : : // these values will not change for subsequent function calls
245 [ # # ]: 0 : if( !been_here ) {
246 [ # # ]: 0 : for ( n = 2; n <= nmax; n++ ) {
247 : 0 : root[n] = sqrt((2.0*n-1) / (2.0*n));
248 : : }
249 : :
250 [ # # ]: 0 : for ( m = 0; m <= nmax; m++ ) {
251 : 0 : double mm = m*m;
252 [ # # ][ # # ]: 0 : for ( n = MAX(m + 1, 2); n <= nmax; n++ ) {
253 : 0 : roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
254 : 0 : roots[m][n][1] = 1.0 / sqrt( n*n - mm);
255 : : }
256 : : }
257 : 0 : been_here = 1;
258 : : }
259 : :
260 [ # # ]: 0 : for ( n=2; n <= nmax; n++ ) {
261 : : // double root = sqrt((2.0*n-1) / (2.0*n));
262 : 0 : P[n][n] = P[n-1][n-1] * s * root[n];
263 : : DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
264 : 0 : root[n];
265 : : }
266 : :
267 : : /* lower triangle */
268 [ # # ]: 0 : for ( m = 0; m <= nmax; m++ ) {
269 : : // double mm = m*m;
270 [ # # ][ # # ]: 0 : for ( n = MAX(m + 1, 2); n <= nmax; n++ ) {
271 : : // double root1 = sqrt((n-1)*(n-1) - mm);
272 : : // double root2 = 1.0 / sqrt( n*n - mm);
273 : : P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
274 : : P[n-2][m] * roots[m][n][0]) *
275 : 0 : roots[m][n][1];
276 : :
277 : : DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
278 : : (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) *
279 : 0 : roots[m][n][1];
280 : : }
281 : : }
282 : :
283 : : /* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
284 : : achieved by adjusting the coefficients at time t0 for linear secular variation */
285 : : /* WMM2005 */
286 : 0 : yearfrac = (dat - date0_wmm2005) / 365.25;
287 [ # # ]: 0 : for ( n = 1; n <= nmax; n++ ) {
288 [ # # ]: 0 : for ( m = 0; m <= nmax; m++ ) {
289 : 0 : gnm[n][m] = gnm_wmm2005[n][m] + yearfrac * gtnm_wmm2005[n][m];
290 : 0 : hnm[n][m] = hnm_wmm2005[n][m] + yearfrac * htnm_wmm2005[n][m];
291 : : }
292 : : }
293 : :
294 : : /* compute sm (sin(m lon) and cm (cos(m lon)) */
295 [ # # ]: 0 : for ( m = 0; m <= nmax; m++ ) {
296 : 0 : sm[m] = sin(m * lon);
297 : 0 : cm[m] = cos(m * lon);
298 : : }
299 : :
300 : : /* compute B fields */
301 : 0 : B_r = 0.0;
302 : 0 : B_theta = 0.0;
303 : 0 : B_phi = 0.0;
304 : 0 : fn_0 = r_0/r;
305 : 0 : fn = fn_0 * fn_0;
306 : :
307 [ # # ]: 0 : for ( n = 1; n <= nmax; n++ ) {
308 : 0 : double c1_n=0;
309 : 0 : double c2_n=0;
310 : 0 : double c3_n=0;
311 [ # # ]: 0 : for ( m = 0; m <= n; m++ ) {
312 : 0 : double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
313 : 0 : c1_n=c1_n + tmp * P[n][m];
314 : 0 : c2_n=c2_n + tmp * DP[n][m];
315 : 0 : c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
316 : : }
317 : : // fn=pow(r_0/r,n+2.0);
318 : 0 : fn *= fn_0;
319 : 0 : B_r = B_r + (n + 1) * c1_n * fn;
320 : 0 : B_theta = B_theta - c2_n * fn;
321 : 0 : B_phi = B_phi + c3_n * fn * inv_s;
322 : : }
323 : :
324 : : /* Find geodetic field components: */
325 : 0 : psi = theta - ((M_PI / 2.0) - lat);
326 : 0 : sinpsi = sin(psi);
327 : 0 : cospsi = cos(psi);
328 : 0 : X = -B_theta * cospsi - B_r * sinpsi;
329 : 0 : Y = B_phi;
330 : 0 : Z = B_theta * sinpsi - B_r * cospsi;
331 : :
332 : 0 : field[0]=B_r;
333 : 0 : field[1]=B_theta;
334 : 0 : field[2]=B_phi;
335 : 0 : field[3]=X;
336 : 0 : field[4]=Y;
337 : 0 : field[5]=Z; /* output fields */
338 : :
339 : : /* find variation in radians */
340 : : /* return zero variation at magnetic pole X=Y=0. */
341 : : /* E is positive */
342 [ # # ][ # # ]: 0 : return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
343 : : }
|
