Geo Stars Library  0.9.3
Geodetic and Astrometry library
geoMag.c
Go to the documentation of this file.
1 /*!
2 \file geoMag.c
3 \brief This file contains the source for deriving geomagnetic variables
4  from geodetic coordinates
5 */
6 #include <math.h>
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <time.h>
10 
11 #include "geoStars.h"
12 
13 
14 #define wmmData wmmDta2010
15 #define EPOCH 2010.0
16 
17 /*!
18 \brief This structure array contains the all of the coefficients for WMM-2000 used in this library
19 
20 
21 This is documentation that came with the WMM. It is included here for completeness.
22 \verbatim
23 C MODEL: THE WMM SERIES GEOMAGNETIC MODELS ARE COMPOSED
24 C OF TWO PARTS: THE MAIN FIELD MODEL, WHICH IS
25 C VALID AT THE BASE EPOCH OF THE CURRENT MODEL AND
26 C A SECULAR VARIATION MODEL, WHICH ACCOUNTS FOR SLOW
27 C TEMPORAL VARIATIONS IN THE MAIN GEOMAGNETIC FIELD
28 C FROM THE BASE EPOCH TO A MAXIMUM OF 5 YEARS BEYOND
29 C THE BASE EPOCH. FOR EXAMPLE, THE BASE EPOCH OF
30 C THE WMM-2005 MODEL IS 2005.0. THIS MODEL IS THEREFORE
31 C CONSIDERED VALID BETWEEN 2005.0 AND 2010.0. THE
32 C COMPUTED MAGNETIC PARAMETERS ARE REFERENCED TO THE
33 C WGS-84 ELLIPSOID.
34 C
35 C***********************************************************************
36 C
37 C ACCURACY: IN OCEAN AREAS AT THE EARTH'S SURFACE OVER THE
38 C ENTIRE 5 YEAR LIFE OF THE DEGREE AND ORDER 12
39 C SPHERICAL HARMONIC MODEL WMM-2005, THE ESTIMATED
40 C MAXIMUM RMS ERRORS FOR THE VARIOUS MAGNETIC COMPONENTS
41 C ARE:
42 C
43 C DEC - 0.5 Degrees
44 C DIP - 0.5 Degrees
45 C TI - 280.0 nanoTeslas (nT)
46 C GV - 0.5 Degrees
47 C
48 C OTHER MAGNETIC COMPONENTS THAT CAN BE DERIVED FROM
49 C THESE FOUR BY SIMPLE TRIGONOMETRIC RELATIONS WILL
50 C HAVE THE FOLLOWING APPROXIMATE ERRORS OVER OCEAN AREAS:
51 C
52 C X - 140 nT (North)
53 C Y - 140 nT (East)
54 C Z - 200 nT (Vertical) Positive is down
55 C H - 200 nT (Horizontal)
56 C
57 C OVER LAND THE MAXIMUM RMS ERRORS ARE EXPECTED TO BE
58 C HIGHER, ALTHOUGH THE RMS ERRORS FOR DEC, DIP, AND GV
59 C ARE STILL ESTIMATED TO BE LESS THAN 1.0 DEGREE, FOR
60 C THE ENTIRE 5-YEAR LIFE OF THE MODEL AT THE EARTH's
61 C SURFACE. THE OTHER COMPONENT ERRORS OVER LAND ARE
62 C MORE DIFFICULT TO ESTIMATE AND SO ARE NOT GIVEN.
63 C
64 C THE ACCURACY AT ANY GIVEN TIME FOR ALL OF THESE
65 C GEOMAGNETIC PARAMETERS DEPENDS ON THE GEOMAGNETIC
66 C LATITUDE. THE ERRORS ARE LEAST FROM THE EQUATOR TO
67 C MID-LATITUDES AND GREATEST NEAR THE MAGNETIC POLES.
68 C
69 C IT IS VERY IMPORTANT TO NOTE THAT A DEGREE AND
70 C ORDER 12 MODEL, SUCH AS WMM-2005, DESCRIBES ONLY
71 C THE LONG WAVELENGTH SPATIAL MAGNETIC FLUCTUATIONS
72 C DUE TO EARTH'S CORE. NOT INCLUDED IN THE WMM SERIES
73 C MODELS ARE INTERMEDIATE AND SHORT WAVELENGTH
74 C SPATIAL FLUCTUATIONS OF THE GEOMAGNETIC FIELD
75 C WHICH ORIGINATE IN THE EARTH'S MANTLE AND CRUST.
76 C CONSEQUENTLY, ISOLATED ANGULAR ERRORS AT VARIOUS
77 C POSITIONS ON THE SURFACE (PRIMARILY OVER LAND, IN
78 C CONTINENTAL MARGINS AND OVER OCEANIC SEAMOUNTS,
79 C RIDGES AND TRENCHES) OF SEVERAL DEGREES MAY BE
80 C EXPECTED. ALSO NOT INCLUDED IN THE MODEL ARE
81 C NONSECULAR TEMPORAL FLUCTUATIONS OF THE GEOMAGNETIC
82 C FIELD OF MAGNETOSPHERIC AND IONOSPHERIC ORIGIN.
83 C DURING MAGNETIC STORMS, TEMPORAL FLUCTUATIONS CAN
84 C CAUSE SUBSTANTIAL DEVIATIONS OF THE GEOMAGNETIC
85 C FIELD FROM MODEL VALUES. IN ARCTIC AND ANTARCTIC
86 C REGIONS, AS WELL AS IN EQUATORIAL REGIONS, DEVIATIONS
87 C FROM MODEL VALUES ARE BOTH FREQUENT AND PERSISTENT.
88 C
89 C IF THE REQUIRED DECLINATION ACCURACY IS MORE
90 C STRINGENT THAN THE WMM SERIES OF MODELS PROVIDE, THEN
91 C THE USER IS ADVISED TO REQUEST SPECIAL (REGIONAL OR
92 C LOCAL) SURVEYS BE PERFORMED AND MODELS PREPARED.
93 C REQUESTS OF THIS NATURE SHOULD BE MADE TO NIMA
94 C AT THE ADDRESS ABOVE.
95 \endverbatim
96 */
97 
98 static WMM_DATA wmmDta2000[] =
99 {
100 /*-----n---m------gnm--------hnm---------dgnm---------dhnm---*/
101 { 1, 0, -29616.0, 0.0, 14.7, 0.0 },
102 { 1, 1, -1722.7, 5194.5, 11.1, -20.4 },
103 { 2, 0, -2266.7, 0.0, -13.6, 0.0 },
104 { 2, 1, 3070.2, -2484.8, -0.7, -21.5 },
105 { 2, 2, 1677.6, -467.9, -1.8, -9.6 },
106 { 3, 0, 1322.4, 0.0, 0.3, 0.0 },
107 { 3, 1, -2291.5, -224.7, -4.3, 6.4 },
108 { 3, 2, 1255.9, 293.0, 0.9, -1.3 },
109 { 3, 3, 724.8, -486.5, -8.4, -13.3 },
110 { 4, 0, 932.1, 0.0, -1.6, 0.0 },
111 { 4, 1, 786.3, 273.3, 0.9, 2.3 },
112 { 4, 2, 250.6, -227.9, -7.6, 0.7 },
113 { 4, 3, -401.5, 120.9, 2.2, 3.7 },
114 { 4, 4, 106.2, -302.7, -3.2, -0.5 },
115 { 5, 0, -211.9, 0.0, -0.9, 0.0 },
116 { 5, 1, 351.6, 42.0, -0.2, 0.0 },
117 { 5, 2, 220.8, 173.8, -2.5, 2.1 },
118 { 5, 3, -134.5, -135.0, -2.7, 2.3 },
119 { 5, 4, -168.8, -38.6, -0.9, 3.1 },
120 { 5, 5, -13.3, 105.2, 1.7, 0.0 },
121 { 6, 0, 73.8, 0.0, 1.2, 0.0 },
122 { 6, 1, 68.2, -17.4, 0.2, -0.3 },
123 { 6, 2, 74.1, 61.2, 1.7, -1.7 },
124 { 6, 3, -163.5, 63.2, 1.6, -0.9 },
125 { 6, 4, -3.8, -62.9, -0.1, -1.0 },
126 { 6, 5, 17.1, 0.2, -0.3, -0.1 },
127 { 6, 6, -85.1, 43.0, 0.8, 1.9 },
128 { 7, 0, 77.4, 0.0, -0.4, 0.0 },
129 { 7, 1, -73.9, -62.3, -0.8, 1.4 },
130 { 7, 2, 2.2, -24.5, -0.2, 0.2 },
131 { 7, 3, 35.7, 8.9, 1.1, 0.7 },
132 { 7, 4, 7.3, 23.4, 0.4, 0.4 },
133 { 7, 5, 5.2, 15.0, 0.0, -0.3 },
134 { 7, 6, 8.4, -27.6, -0.2, -0.8 },
135 { 7, 7, -1.5, -7.8, -0.2, -0.1 },
136 { 8, 0, 23.3, 0.0, -0.3, 0.0 },
137 { 8, 1, 7.3, 12.4, 0.6, -0.5 },
138 { 8, 2, -8.5, -20.8, -0.8, 0.1 },
139 { 8, 3, -6.6, 8.4, 0.3, -0.2 },
140 { 8, 4, -16.9, -21.2, -0.2, 0.0 },
141 { 8, 5, 8.6, 15.5, 0.5, 0.1 },
142 { 8, 6, 4.9, 9.1, 0.0, -0.1 },
143 { 8, 7, -7.8, -15.5, -0.6, 0.3 },
144 { 8, 8, -7.6, -5.4, 0.1, 0.2 },
145 { 9, 0, 5.7, 0.0, 0.0, 0.0 },
146 { 9, 1, 8.5, -20.4, 0.0, 0.0 },
147 { 9, 2, 2.0, 13.9, 0.0, 0.0 },
148 { 9, 3, -9.8, 12.0, 0.0, 0.0 },
149 { 9, 4, 7.6, -6.2, 0.0, 0.0 },
150 { 9, 5, -7.0, -8.6, 0.0, 0.0 },
151 { 9, 6, -2.0, 9.4, 0.0, 0.0 },
152 { 9, 7, 9.2, 5.0, 0.0, 0.0 },
153 { 9, 8, -2.2, -8.4, 0.0, 0.0 },
154 { 9, 9, -6.6, 3.2, 0.0, 0.0 },
155 { 10, 0, -2.2, 0.0, 0.0, 0.0 },
156 { 10, 1, -5.7, 0.9, 0.0, 0.0 },
157 { 10, 2, 1.6, -0.7, 0.0, 0.0 },
158 { 10, 3, -3.7, 3.9, 0.0, 0.0 },
159 { 10, 4, -0.6, 4.8, 0.0, 0.0 },
160 { 10, 5, 4.1, -5.3, 0.0, 0.0 },
161 { 10, 6, 2.2, -1.0, 0.0, 0.0 },
162 { 10, 7, 2.2, -2.4, 0.0, 0.0 },
163 { 10, 8, 4.6, 1.3, 0.0, 0.0 },
164 { 10, 9, 2.3, -2.3, 0.0, 0.0 },
165 { 10, 10, 0.1, -6.4, 0.0, 0.0 },
166 { 11, 0, 3.3, 0.0, 0.0, 0.0 },
167 { 11, 1, -1.1, -1.5, 0.0, 0.0 },
168 { 11, 2, -2.4, 0.7, 0.0, 0.0 },
169 { 11, 3, 2.6, -1.1, 0.0, 0.0 },
170 { 11, 4, -1.3, -2.3, 0.0, 0.0 },
171 { 11, 5, -1.7, 1.3, 0.0, 0.0 },
172 { 11, 6, -0.6, -0.6, 0.0, 0.0 },
173 { 11, 7, 0.4, -2.8, 0.0, 0.0 },
174 { 11, 8, 0.7, -1.6, 0.0, 0.0 },
175 { 11, 9, -0.3, -0.1, 0.0, 0.0 },
176 { 11, 10, 2.3, -1.9, 0.0, 0.0 },
177 { 11, 11, 4.2, 1.4, 0.0, 0.0 },
178 { 12, 0, -1.5, 0.0, 0.0, 0.0 },
179 { 12, 1, -0.2, -1.0, 0.0, 0.0 },
180 { 12, 2, -0.3, 0.7, 0.0, 0.0 },
181 { 12, 3, 0.5, 2.2, 0.0, 0.0 },
182 { 12, 4, 0.2, -2.5, 0.0, 0.0 },
183 { 12, 5, 0.9, -0.2, 0.0, 0.0 },
184 { 12, 6, -1.4, 0.0, 0.0, 0.0 },
185 { 12, 7, 0.6, -0.2, 0.0, 0.0 },
186 { 12, 8, -0.6, 0.0, 0.0, 0.0 },
187 { 12, 9, -1.0, 0.2, 0.0, 0.0 },
188 { 12, 10, -0.3, -0.9, 0.0, 0.0 },
189 { 12, 11, 0.3, -0.2, 0.0, 0.0 },
190 { 12, 12, 0.4, 1.0, 0.0, 0.0 }
191 };
192 
193 
194 /*!
195 \brief This structure array contains the all of the coefficients for WMM-2005 used in this library
196 
197 */
198 
199 static WMM_DATA wmmDta2005[] =
200 {
201 /*-----n---m------gnm--------hnm---------dgnm---------dhnm---*/
202 { 1, 0, -29556.8, 0.0, 8.0, 0.0},
203 { 1, 1, -1671.7, 5079.8, 10.6, -20.9},
204 { 2, 0, -2340.6, 0.0, -15.1, 0.0},
205 { 2, 1, 3046.9, -2594.7, -7.8, -23.2},
206 { 2, 2, 1657.0, -516.7, -0.8, -14.6},
207 { 3, 0, 1335.4, 0.0, 0.4, 0.0},
208 { 3, 1, -2305.1, -199.9, -2.6, 5.0},
209 { 3, 2, 1246.7, 269.3, -1.2, -7.0},
210 { 3, 3, 674.0, -524.2, -6.5, -0.6},
211 { 4, 0, 919.8, 0.0, -2.5, 0.0},
212 { 4, 1, 798.1, 281.5, 2.8, 2.2},
213 { 4, 2, 211.3, -226.0, -7.0, 1.6},
214 { 4, 3, -379.4, 145.8, 6.2, 5.8},
215 { 4, 4, 100.0, -304.7, -3.8, 0.1},
216 { 5, 0, -227.4, 0.0, -2.8, 0.0},
217 { 5, 1, 354.6, 42.4, 0.7, 0.0},
218 { 5, 2, 208.7, 179.8, -3.2, 1.7},
219 { 5, 3, -136.5, -123.0, -1.1, 2.1},
220 { 5, 4, -168.3, -19.5, 0.1, 4.8},
221 { 5, 5, -14.1, 103.6, -0.8, -1.1},
222 { 6, 0, 73.2, 0.0, -0.7, 0.0},
223 { 6, 1, 69.7, -20.3, 0.4, -0.6},
224 { 6, 2, 76.7, 54.7, -0.3, -1.9},
225 { 6, 3, -151.2, 63.6, 2.3, -0.4},
226 { 6, 4, -14.9, -63.4, -2.1, -0.5},
227 { 6, 5, 14.6, -0.1, -0.6, -0.3},
228 { 6, 6, -86.3, 50.4, 1.4, 0.7},
229 { 7, 0, 80.1, 0.0, 0.2, 0.0},
230 { 7, 1, -74.5, -61.5, -0.1, 0.6},
231 { 7, 2, -1.4, -22.4, -0.3, 0.4},
232 { 7, 3, 38.5, 7.2, 1.1, 0.2},
233 { 7, 4, 12.4, 25.4, 0.6, 0.3},
234 { 7, 5, 9.5, 11.0, 0.5, -0.8},
235 { 7, 6, 5.7, -26.4, -0.4, -0.2},
236 { 7, 7, 1.8, -5.1, 0.6, 0.1},
237 { 8, 0, 24.9, 0.0, 0.1, 0.0},
238 { 8, 1, 7.7, 11.2, 0.3, -0.2},
239 { 8, 2, -11.6, -21.0, -0.4, 0.1},
240 { 8, 3, -6.9, 9.6, 0.3, 0.3},
241 { 8, 4, -18.2, -19.8, -0.3, 0.4},
242 { 8, 5, 10.0, 16.1, 0.2, 0.1},
243 { 8, 6, 9.2, 7.7, 0.4, -0.2},
244 { 8, 7, -11.6, -12.9, -0.7, 0.4},
245 { 8, 8, -5.2, -0.2, 0.4, 0.4},
246 { 9, 0, 5.6, 0.0, 0.0, 0.0},
247 { 9, 1, 9.9, -20.1, 0.0, 0.0},
248 { 9, 2, 3.5, 12.9, 0.0, 0.0},
249 { 9, 3, -7.0, 12.6, 0.0, 0.0},
250 { 9, 4, 5.1, -6.7, 0.0, 0.0},
251 { 9, 5, -10.8, -8.1, 0.0, 0.0},
252 { 9, 6, -1.3, 8.0, 0.0, 0.0},
253 { 9, 7, 8.8, 2.9, 0.0, 0.0},
254 { 9, 8, -6.7, -7.9, 0.0, 0.0},
255 { 9, 9, -9.1, 6.0, 0.0, 0.0},
256 { 10, 0, -2.3, 0.0, 0.0, 0.0},
257 { 10, 1, -6.3, 2.4, 0.0, 0.0},
258 { 10, 2, 1.6, 0.2, 0.0, 0.0},
259 { 10, 3, -2.6, 4.4, 0.0, 0.0},
260 { 10, 4, 0.0, 4.8, 0.0, 0.0},
261 { 10, 5, 3.1, -6.5, 0.0, 0.0},
262 { 10, 6, 0.4, -1.1, 0.0, 0.0},
263 { 10, 7, 2.1, -3.4, 0.0, 0.0},
264 { 10, 8, 3.9, -0.8, 0.0, 0.0},
265 { 10, 9, -0.1, -2.3, 0.0, 0.0},
266 { 10, 10, -2.3, -7.9, 0.0, 0.0},
267 { 11, 0, 2.8, 0.0, 0.0, 0.0},
268 { 11, 1, -1.6, 0.3, 0.0, 0.0},
269 { 11, 2, -1.7, 1.2, 0.0, 0.0},
270 { 11, 3, 1.7, -0.8, 0.0, 0.0},
271 { 11, 4, -0.1, -2.5, 0.0, 0.0},
272 { 11, 5, 0.1, 0.9, 0.0, 0.0},
273 { 11, 6, -0.7, -0.6, 0.0, 0.0},
274 { 11, 7, 0.7, -2.7, 0.0, 0.0},
275 { 11, 8, 1.8, -0.9, 0.0, 0.0},
276 { 11, 9, 0.0, -1.3, 0.0, 0.0},
277 { 11, 10, 1.1, -2.0, 0.0, 0.0},
278 { 11, 11, 4.1, -1.2, 0.0, 0.0},
279 { 12, 0, -2.4, 0.0, 0.0, 0.0},
280 { 12, 1, -0.4, -0.4, 0.0, 0.0},
281 { 12, 2, 0.2, 0.3, 0.0, 0.0},
282 { 12, 3, 0.8, 2.4, 0.0, 0.0},
283 { 12, 4, -0.3, -2.6, 0.0, 0.0},
284 { 12, 5, 1.1, 0.6, 0.0, 0.0},
285 { 12, 6, -0.5, 0.3, 0.0, 0.0},
286 { 12, 7, 0.4, 0.0, 0.0, 0.0},
287 { 12, 8, -0.3, 0.0, 0.0, 0.0},
288 { 12, 9, -0.3, 0.3, 0.0, 0.0},
289 { 12, 10, -0.1, -0.9, 0.0, 0.0},
290 { 12, 11, -0.3, -0.4, 0.0, 0.0},
291 { 12, 12, -0.1, 0.8, 0.0, 0.0}
292 };
293 
294 
295 
296 
297 static WMM_DATA wmmDta2010[] =
298 {
299 /*-----n---m------gnm--------hnm---------dgnm---------dhnm---*/
300 { 1, 0, -29496.6, 0.0, 11.6, 0.0},
301 { 1, 1, -1586.3, 4944.4, 16.5, -25.9},
302 { 2, 0, -2396.6, 0.0, -12.1, 0.0},
303 { 2, 1, 3026.1, -2707.7, -4.4, -22.5},
304 { 2, 2, 1668.6, -576.1, 1.9, -11.8},
305 { 3, 0, 1340.1, 0.0, 0.4, 0.0},
306 { 3, 1, -2326.2, -160.2, -4.1, 7.3},
307 { 3, 2, 1231.9, 251.9, -2.9, -3.9},
308 { 3, 3, 634.0, -536.6, -7.7, -2.6},
309 { 4, 0, 912.6, 0.0, -1.8, 0.0},
310 { 4, 1, 808.9, 286.4, 2.3, 1.1},
311 { 4, 2, 166.7, -211.2, -8.7, 2.7},
312 { 4, 3, -357.1, 164.3, 4.6, 3.9},
313 { 4, 4, 89.4, -309.1, -2.1, -0.8},
314 { 5, 0, -230.9, 0.0, -1.0, 0.0},
315 { 5, 1, 357.2, 44.6, 0.6, 0.4},
316 { 5, 2, 200.3, 188.9, -1.8, 1.8},
317 { 5, 3, -141.1, -118.2, -1.0, 1.2},
318 { 5, 4, -163.0, 0.0, 0.9, 4.0},
319 { 5, 5, -7.8, 100.9, 1.0, -0.6},
320 { 6, 0, 72.8, 0.0, -0.2, 0.0},
321 { 6, 1, 68.6, -20.8, -0.2, -0.2},
322 { 6, 2, 76.0, 44.1, -0.1, -2.1},
323 { 6, 3, -141.4, 61.5, 2.0, -0.4},
324 { 6, 4, -22.8, -66.3, -1.7, -0.6},
325 { 6, 5, 13.2, 3.1, -0.3, 0.5},
326 { 6, 6, -77.9, 55.0, 1.7, 0.9},
327 { 7, 0, 80.5, 0.0, 0.1, 0.0},
328 { 7, 1, -75.1, -57.9, -0.1, 0.7},
329 { 7, 2, -4.7, -21.1, -0.6, 0.3},
330 { 7, 3, 45.3, 6.5, 1.3, -0.1},
331 { 7, 4, 13.9, 24.9, 0.4, -0.1},
332 { 7, 5, 10.4, 7.0, 0.3, -0.8},
333 { 7, 6, 1.7, -27.7, -0.7, -0.3},
334 { 7, 7, 4.9, -3.3, 0.6, 0.3},
335 { 8, 0, 24.4, 0.0, -0.1, 0.0},
336 { 8, 1, 8.1, 11.0, 0.1, -0.1},
337 { 8, 2, -14.5, -20.0, -0.6, 0.2},
338 { 8, 3, -5.6, 11.9, 0.2, 0.4},
339 { 8, 4, -19.3, -17.4, -0.2, 0.4},
340 { 8, 5, 11.5, 16.7, 0.3, 0.1},
341 { 8, 6, 10.9, 7.0, 0.3, -0.1},
342 { 8, 7, -14.1, -10.8, -0.6, 0.4},
343 { 8, 8, -3.7, 1.7, 0.2, 0.3},
344 { 9, 0, 5.4, 0.0, 0.0, 0.0},
345 { 9, 1, 9.4, -20.5, -0.1, 0.0},
346 { 9, 2, 3.4, 11.5, 0.0, -0.2},
347 { 9, 3, -5.2, 12.8, 0.3, 0.0},
348 { 9, 4, 3.1, -7.2, -0.4, -0.1},
349 { 9, 5, -12.4, -7.4, -0.3, 0.1},
350 { 9, 6, -0.7, 8.0, 0.1, 0.0},
351 { 9, 7, 8.4, 2.1, -0.1, -0.2},
352 { 9, 8, -8.5, -6.1, -0.4, 0.3},
353 { 9, 9, -10.1, 7.0, -0.2, 0.2},
354 { 10, 0, -2.0, 0.0, 0.0, 0.0},
355 { 10, 1, -6.3, 2.8, 0.0, 0.1},
356 { 10, 2, 0.9, -0.1, -0.1, -0.1},
357 { 10, 3, -1.1, 4.7, 0.2, 0.0},
358 { 10, 4, -0.2, 4.4, 0.0, -0.1},
359 { 10, 5, 2.5, -7.2, -0.1, -0.1},
360 { 10, 6, -0.3, -1.0, -0.2, 0.0},
361 { 10, 7, 2.2, -3.9, 0.0, -0.1},
362 { 10, 8, 3.1, -2.0, -0.1, -0.2},
363 { 10, 9, -1.0, -2.0, -0.2, 0.0},
364 { 10, 10, -2.8, -8.3, -0.2, -0.1},
365 { 11, 0, 3.0, 0.0, 0.0, 0.0},
366 { 11, 1, -1.5, 0.2, 0.0, 0.0},
367 { 11, 2, -2.1, 1.7, 0.0, 0.1},
368 { 11, 3, 1.7, -0.6, 0.1, 0.0},
369 { 11, 4, -0.5, -1.8, 0.0, 0.1},
370 { 11, 5, 0.5, 0.9, 0.0, 0.0},
371 { 11, 6, -0.8, -0.4, 0.0, 0.1},
372 { 11, 7, 0.4, -2.5, 0.0, 0.0},
373 { 11, 8, 1.8, -1.3, 0.0, -0.1},
374 { 11, 9, 0.1, -2.1, 0.0, -0.1},
375 { 11, 10, 0.7, -1.9, -0.1, 0.0},
376 { 11, 11, 3.8, -1.8, 0.0, -0.1},
377 { 12, 0, -2.2, 0.0, 0.0, 0.0},
378 { 12, 1, -0.2, -0.9, 0.0, 0.0},
379 { 12, 2, 0.3, 0.3, 0.1, 0.0},
380 { 12, 3, 1.0, 2.1, 0.1, 0.0},
381 { 12, 4, -0.6, -2.5, -0.1, 0.0},
382 { 12, 5, 0.9, 0.5, 0.0, 0.0},
383 { 12, 6, -0.1, 0.6, 0.0, 0.1},
384 { 12, 7, 0.5, 0.0, 0.0, 0.0},
385 { 12, 8, -0.4, 0.1, 0.0, 0.0},
386 { 12, 9, -0.4, 0.3, 0.0, 0.0},
387 { 12, 10, 0.2, -0.9, 0.0, 0.0},
388 { 12, 11, -0.8, -0.2, -0.1, 0.0},
389 { 12, 12, 0.0, 0.9, 0.1, 0.0}
390 };
391 
392 static int started = 0; //! \internal Flag for limiting multiple initializations
393 
394 #define snorm p //! \internal redefinition to fit original fortran code
395 #define MAXDEG 12 //! \internal Maximum degrees
396 
397 
398 // Static Prototypes
399 static void geoMagInit(void);
400 
401 
402 
403 static double
404 c[13][13], //!< C - GAUSS COEFFICIENTS OF MAIN GEOMAGNETIC MODEL (NT)
405 cd[13][13], //!< CD - GAUSS COEFFICIENTS OF SECULAR GEOMAGNETIC MODEL (NT/YR)
406 tc[13][13], //!< TC - TIME ADJUSTED GEOMAGNETIC GAUSS COEFFICIENTS (NT)
407 dp[13][13], //!< DP(N,M)- THETA DERIVATIVE OF P(N,M) (UNNORMALIZED)
408 snorm[13][13], //!< SNORM - SCHMIDT NORMALIZATION FACTORS
409 sp[13], //!< SP(M) - SINE OF (M*SPHERICAL COORD. LONGITUDE)
410 cp[13], //!< CP(M) - COSINE OF (M*SPHERICAL COORD. LONGITUDE)
411 fn[13], fm[13],
412 pp[13], //!< PP(N) - ASSOCIATED LEGENDRE POLYNOMIALS FOR M=1 (UNNORMALIZED)
413 k[13][13],
414 epoch;
415 // otime, //!< OTIME - TIME ON PREVIOUS CALL TO GEOMAG (YRS)
416 // oalt, //!< OALT - GEODETIC ALTITUDE ON PREVIOUS CALL TO GEOMAG (YRS)
417 // olat, //!< OLAT - GEODETIC LATITUDE ON PREVIOUS CALL TO GEOMAG (DEG.)
418 // olon; //!< OLON - GEODETIC LONGITUDE ON PREVIOUS CALL TO GEOMAG (DEG.)
419 
420 static int maxord; //!< MAXORD - MAXIMUM ORDER OF SPHERICAL HARMONIC MODEL
421 
422 /*************************************************************************/
423 
424 static void geoMagInit(void)
425 {
426  int n,m,i,j;
427  double flnmj,n2m1,nm12,n2m1stuff;
428 
429  /* Only do this once */
430  if(started) return;
431  started = 1;
432 
433  // INITIALIZE CONSTANTS
434  maxord = MAXDEG;
435  sp[0] = 0.0;
436  cp[0] = p[0][0] = pp[0] = 1.0;
437  dp[0][0] = 0.0;
438 
439  // READ WORLD MAGNETIC MODEL SPHERICAL HARMONIC COEFFICIENTS
440  c[0][0] = 0.0;
441  cd[0][0] = 0.0;
442 
443  // Load the WMM Coefficients
444  epoch=EPOCH;
445  for( i=0;i<90;i++)
446  {
447  n=wmmData[i].n;
448  m=wmmData[i].m;
449  if (m <= n)
450  {
451  c[m][n] = wmmData[i].gnm;
452  cd[m][n] = wmmData[i].dgnm;
453  if (m != 0)
454  {
455  c[n][m-1] = wmmData[i].hnm;
456  cd[n][m-1] = wmmData[i].dhnm;
457  }
458  }
459  } //end for
460 
461 
462  /* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */
463 
464  snorm[0][0] = 1.0;
465  for (n=1; n<=maxord; n++)
466  {
467  nm12 = (n-1)*(n-1);
468  n2m1 = (double)2*n-1;
469  n2m1stuff = (double)(n2m1*(2*n-3));
470  snorm[0][n] = snorm[0][n-1]*n2m1/(double)n;
471  j = 2;
472  for (m=0;m<=n;m++)
473  {
474  k[m][n] = (double)(nm12-(m*m))/n2m1stuff;
475  if (m > 0)
476  {
477  flnmj = (double)((n-m+1)*j)/(double)(n+m);
478  snorm[m][n] = snorm[m-1][n]*sqrt(flnmj);
479  j = 1;
480  c[n][m-1] = snorm[m][n]*c[n][m-1];
481  cd[n][m-1] = snorm[m][n]*cd[n][m-1];
482  }
483  c[m][n] = snorm[m][n]*c[m][n];
484  cd[m][n] = snorm[m][n]*cd[m][n];
485  } //end for m
486  fn[n] = (double)(n+1);
487  fm[n] = (double)n;
488  } //end for n
489  k[1][1] = 0.0;
490 
491 
492 }
493 
494 
495 
496 // Main GeoMag routine
497 /*! \brief This routine computes all of the relevant geomagnetic data
498 
499 \param double alt : in: altitude in meters
500 \param double glat : in: latitude in decimal degrees
501 \param double glon : in: longitude in decimal degrees
502 \param double time : in: time in decimal years
503 \param double *dec :out: declination in degrees
504 \param double *dip :out: dip in degrees
505 \param double *ti :out: total intensity in nT
506 \param double *gv :out: geomagnetic grid variation
507 \param double *adec :out: Annual declination in degrees
508 \param double *adip :out: Annual dip in degrees
509 \param double *ati :out: Annual total intensity in nT
510 \param double *x :out: X Component of the magnetic field
511 \param double *y :out: Y Component of the magnetic field
512 \param double *z :out: Z Component of the magnetic field
513 \param double *h :out: h Component of the magnetic field
514 \param double *ax :out: Annual X Component of the magnetic field
515 \param double *ay :out: Annual Y Component of the magnetic field
516 \param double *az :out: Annual Z Component of the magnetic field
517 \param double *ah :out: Annual H Component of the magnetic field
518 \retval GEO_OK on success
519 \retval GEO_ERROR on error
520 
521 {
522  */
523 int DLL_CALLCONV geoMag( double alt, double glat, double glon, double time,
524 double *dec, double *dip, double *ti, double *gv,
525 double *adec, double *adip, double *ati,
526 double *x, double *y, double *z, double *h,
527 double *ax, double *ay, double *az, double *ah)
528 {
529  double time1;
530  double dec2, dip2, ti2, rdec, rdip, rdec2, rdip2;
531  double x2,y2,z2,h2;
532  double c_rdip,c_rdip2;
533  // double rTd=0.017453292;
534 
535  geoMagInit();
536 
537  geomg1(alt,glat,glon,time,dec,dip,ti,gv);
538  time1 = time + 1.0; // add a year to time
539  geomg1(alt,glat,glon,time1,&dec2,&dip2,&ti2,gv);
540 
541 
542 
543  /*COMPUTE X, Y, Z, AND H COMPONENTS OF THE MAGNETIC FIELD*/
544  rdec = *dec * DEG_TO_RAD;
545  rdip = *dip * DEG_TO_RAD;
546  rdec2 = dec2 * DEG_TO_RAD;
547  rdip2 = dip2 * DEG_TO_RAD;
548  c_rdip = cos(rdip);
549  c_rdip2 = cos(rdip2);
550  /*
551  *x= *ti * (cos(rdec) * cos(rdip));
552  *y= *ti * (cos(rdip) * sin(rdec));
553  *z= *ti * (sin(rdip));
554  *h= *ti * (cos(rdip));
555 
556  x2=ti2*(cos(rdec2) * cos(rdip2));
557  y2=ti2*(cos(rdip2) * sin(rdec2));
558  z2=ti2*(sin(rdip2));
559  h2=ti2*(cos(rdip2));
560 */
561  *x= *ti * (cos(rdec) * c_rdip);
562  *y= *ti * (c_rdip * sin(rdec));
563  *z= *ti * (sin(rdip));
564  *h= *ti * (c_rdip);
565 
566  x2=ti2*(cos(rdec2) * c_rdip2);
567  y2=ti2*(c_rdip2 * sin(rdec2));
568  z2=ti2*(sin(rdip2));
569  h2=ti2*(c_rdip2);
570 
571 
572  /* COMPUTE ANNUAL CHANGE FOR TOTAL INTENSITY */
573  *ati = ti2 - *ti;
574 
575  /* COMPUTE ANNUAL CHANGE FOR DIP & DEC (in minutes) */
576  *adip = (dip2 - *dip) * DEG_TO_MIN;
577  *adec = (dec2 - *dec) * DEG_TO_MIN;
578 
579 
580  /* COMPUTE ANNUAL CHANGE FOR X, Y, Z, AND H */
581  *ax = x2 - *x;
582  *ay = y2 - *y;
583  *az = z2 - *z;
584  *ah = h2 - *h;
585 
586  return(GEO_OK);
587 }
588 
589 
590 
591 
592 
593 /*************************************************************************/
594 
595 int DLL_CALLCONV geomg1( double alt, double glat, double glon, double time,
596 double *dec, double *dip, double *ti, double *gv)
597 {
598 
599  /*!
600 This is documentation that came with the WMM. It is included here for completeness.
601 \verbatim
602 // MAXDEG - MAXIMUM DEGREE OF SPHERICAL HARMONIC MODEL (INPUT)
603 // TIME - COMPUTATION TIME (YRS) (INPUT)
604 // (EG. 1 JULY 1995 = 1995.500)
605 // ALT - GEODETIC ALTITUDE (M) (INPUT)
606 // GLAT - GEODETIC LATITUDE (DEG.) (INPUT)
607 // GLON - GEODETIC LONGITUDE (DEG.) (INPUT)
608 // EPOCH - BASE TIME OF GEOMAGNETIC MODEL (YRS)
609 // P(N,M) - ASSOCIATED LEGENDRE POLYNOMIALS (UNNORMALIZED)
610 // DEC - GEOMAGNETIC DECLINATION (DEG.) (OUTPUT)
611 // EAST=POSITIVE ANGLES
612 // WEST=NEGATIVE ANGLES
613 // DIP - GEOMAGNETIC INCLINATION (DEG.) (OUTPUT)
614 // DOWN=POSITIVE ANGLES
615 // UP=NEGATIVE ANGLES
616 // TI - GEOMAGNETIC TOTAL INTENSITY (NT) (OUTPUT)
617 // GV - GEOMAGNETIC GRID VARIATION (DEG.) (OUTPUT)
618 // REFERENCED TO GRID NORTH
619 // GRID NORTH REFERENCED TO 0 MERIDIAN
620 // OF A POLAR STEREOGRAPHIC PROJECTION
621 // (ARCTIC/ANTARCTIC ONLY)
622 // A - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM)
623 // B - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM)
624 // RE - MEAN RADIUS OF IAU-66 ELLIPSOID (KM)
625 // a2, // a * a
626 // b2, // b * b
627 // c2, // a2 - b2
628 // a4, // a2 * a2
629 // b4, // b2 * b2
630 // c4, // a4 - b4
631 // ct, //C CT - COSINE OF (SPHERICAL COORD. LATITUDE)
632 // st, //C ST - SINE OF (SPHERICAL COORD. LATITUDE)
633 // r2,
634 // r, //C R - SPHERICAL COORDINATE RADIAL POSITION (KM)
635 // ca, //C CA - COSINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE
636 // sa, //C SA - SINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE
637 // br, //C BR - RADIAL COMPONENT OF GEOMAGNETIC FIELD (NT)
638 // bt, //C BT - THETA COMPONENT OF GEOMAGNETIC FIELD (NT)
639 // bp, //C BP - PHI COMPONENT OF GEOMAGNETIC FIELD (NT)
640 // bx, //C BX - NORTH GEOMAGNETIC COMPONENT (NT)
641 // by, //C BY - EAST GEOMAGNETIC COMPONENT (NT)
642 // bz, //C BZ - VERTICALLY DOWN GEOMAGNETIC COMPONENT (NT)
643 // bh; //C BH - HORIZONTAL GEOMAGNETIC COMPONENT (NT)
644 \endverbatim
645 */
646 
647  int m,n;
648  double temp1, temp2,d,ar,aor,dt,bpp,
649  par, parp,r,r2,sa,ca,st,ct,br,bt,bp,bx,by,bz,bh;
650 
651  double rlon,rlat,
652  srlon, srlat, crlon, crlat,
653  srlat2, crlat2,
654  a, // A - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM)
655  b, // B - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM)
656  re, // RE - MEAN RADIUS OF IAU-66 ELLIPSOID (KM)
657  a2, // a * a
658  b2, // b * b
659  c2, // a2 - b2
660  a4, // a2 * a2
661  b4, // b2 * b2
662  c4, // a4 - b4
663 
664  q, q1, q2;
665 
666  a = ellips[GEO_DATUM_WE].a; //6378137.0;
667  b = GEO_B(ellips[GEO_DATUM_WE].a, (ellips[GEO_DATUM_WE].f1)); //6356752.3142; //wgs-84
668  re = 6371200.0;
669  a2 = a*a;
670  b2 = b*b;
671  c2 = a2-b2;
672  a4 = a2*a2;
673  b4 = b2*b2;
674  c4 = a4 - b4;
675 
676 
677  dt = time - epoch;
678  // if (dt < 0.0 || dt > 5.0)
679  // printf("\n WARNING geoMag - TIME EXTENDS BEYOND MODEL 5-YEAR LIFE SPAN (dt=%f, time=%f)\n",dt,time);
680 
681  rlon = glon*DEG_TO_RAD;
682  rlat = glat*DEG_TO_RAD;
683  srlon = sin(rlon);
684  srlat = sin(rlat);
685  crlon = cos(rlon);
686  crlat = cos(rlat);
687  srlat2 = srlat*srlat;
688  crlat2 = crlat*crlat;
689  sp[1] = srlon;
690  cp[1] = crlon;
691 
692  /* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */
693 
694  q = sqrt(a2-c2*srlat2);
695  q1 = alt*q;
696  q2 = ((q1+a2)/(q1+b2))*((q1+a2)/(q1+b2));
697  ct = srlat/sqrt(q2*crlat2+srlat2);
698  st = sqrt(1.0-(ct*ct));
699  r2 = (alt*alt)+2.0*q1+(a4-c4*srlat2)/(q*q);
700  r = sqrt(r2);
701  d = sqrt(a2*crlat2+b2*srlat2);
702  ca = (alt+d)/r;
703  sa = c2*crlat*srlat/(r*d);
704 
705  for (m=2; m<=maxord; m++)
706  {
707  sp[m] = sp[1]*cp[m-1]+cp[1]*sp[m-1];
708  cp[m] = cp[1]*cp[m-1]-sp[1]*sp[m-1];
709  }
710 
711  aor = re/r;
712  ar = aor*aor;
713  br = bt = bp = bpp = 0.0;
714  for (n=1; n<=maxord; n++)
715  {
716  ar = ar*aor;
717  for (m=0;m<=n;m++)
718  {
719  /*
720  COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
721  AND DERIVATIVES VIA RECURSION RELATIONS
722  */
723  if (n == m)
724  {
725  p[m][n] = st * p[m-1][n-1];
726  dp[m][n] = st * dp[m-1][n-1] + ct * p[m-1][n-1];
727  goto S50;
728  }
729  if (n == 1 && m == 0)
730  {
731  p[m][n] = ct * p[m][n-1];
732  dp[m][n] = ct * dp[m][n-1] - st * p[m][n-1];
733  goto S50;
734  }
735  if (n > 1 && n != m)
736  {
737  if (m > n-2) p[m][n-2] = 0.0;
738  if (m > n-2) dp[m][n-2] = 0.0;
739  p[m][n] = ct * p[m][n-1] - k[m][n] * p[m][n-2];
740  dp[m][n] = ct * dp[m][n-1] - st * p[m][n-1] - k[m][n] * dp[m][n-2];
741  }
742 S50:
743  // TIME ADJUST THE GAUSS COEFFICIENTS
744  tc[m][n] = c[m][n]+dt*cd[m][n];
745  if (m != 0) tc[n][m-1] = c[n][m-1]+dt*cd[n][m-1];
746 
747  // ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
748  par = ar * p[m][n];
749  if (m == 0)
750  {
751  temp1 = tc[m][n]*cp[m];
752  temp2 = tc[m][n]*sp[m];
753  }
754  else
755  {
756  temp1 = tc[m][n]*cp[m]+tc[n][m-1]*sp[m];
757  temp2 = tc[m][n]*sp[m]-tc[n][m-1]*cp[m];
758  }
759  bt = bt-ar*temp1*dp[m][n];
760  bp += (fm[m]*temp2*par);
761  br += (fn[n]*temp1*par);
762 
763  // SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES
764  if (st == 0.0 && m == 1)
765  {
766  if (n == 1) pp[n] = pp[n-1];
767  else pp[n] = ct*pp[n-1]-k[m][n]*pp[n-2];
768  parp = ar*pp[n];
769  bpp += (fm[m]*temp2*parp);
770  }
771  }
772  }
773  if (st == 0.0) bp = bpp;
774  else bp /= st;
775  /*
776  ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO
777  GEODETIC COORDINATES
778 */
779  bx = -bt*ca-br*sa;
780  by = bp;
781  bz = bt*sa-br*ca;
782  /*
783 ** COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND
784 ** TOTAL INTENSITY (TI)
785 */
786  bh = sqrt((bx*bx)+(by*by));
787  *ti = sqrt((bh*bh)+(bz*bz));
788  *dec = atan2(by,bx)*RAD_TO_DEG;
789  *dip = atan2(bz,bh)*RAD_TO_DEG;
790 
791  /*
792 ** COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT
793 ** GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC
794 ** (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES)
795 **
796 ** OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0
797 */
798  *gv = -999.0;
799  if (fabs(glat) >= 55.)
800  {
801  if (glat > 0.0 && glon >= 0.0) *gv = *dec-glon;
802  if (glat > 0.0 && glon < 0.0) *gv = *dec+fabs(glon);
803  if (glat < 0.0 && glon >= 0.0) *gv = *dec+glon;
804  if (glat < 0.0 && glon < 0.0) *gv = *dec-fabs(glon);
805  if (*gv > +HALF_CIRCLE) *gv -= CIRCLE;
806  if (*gv < -HALF_CIRCLE) *gv += CIRCLE;
807  }
808 
809  return (GEO_OK);
810 
811 }
812 
813 /*************************************************************************/
814 // Short GeoMag routine
815 /*! \brief This routine computes magnetic declination
816 
817 \param double lat : in: latitude in decimal degrees
818 \param double lon : in: longitude in decimal degrees
819 \param double hgt : in: altitude in meters
820 \param int month : in: Month of the year
821 \param int day : in: Day of the month
822 \param int year : in: Year
823 \param double dec :out: declination in degrees
824 */
825 int DLL_CALLCONV geoMagGetDec(double lat, double lon, double hgt,
826 int month, int day, int year, double *dec)
827 {
828  /* Days in each month 1-12 */
829  int jday[] = {0,0,31,59,90,120,151,180,211,242,272,303,333 };
830  double time,dip,ti,gv;
831 
832  /* Determin decimal years from date */
833  time = ((jday[month]+day) / 365.0) + year;
834 
835  /* Determin declination */
836  geoMagInit();
837  geomg1(hgt,lat,lon,time,dec,&dip,&ti,&gv);
838 
839  return (GEO_OK);
840 
841 }
842 
843 /*************************************************************************/
844 // Short GeoMag routine
845 /*! \brief This routine computes and returns magnetic declination for a specific date
846 
847 \param double lat : in: latitude in decimal degrees
848 \param double lon : in: longitude in decimal degrees
849 \param double hgt : in: altitude in meters
850 \param int month : in: Month of the year
851 \param int day : in: Day of the month
852 \param int year : in: Year
853 \retval double dec :out: declination in degrees
854 */
855 double DLL_CALLCONV geoMagGetDecRet(double lat, double lon, double hgt,
856 int month, int day, int year)
857 {
858  /* Days in each month 1-12 */
859  double dec;
860  int jday[] = {0,0,31,59,90,120,151,180,211,242,272,303,333 };
861  double time,dip,ti,gv;
862 
863  /* Determin decimal years from date */
864  time = ((jday[month]+day) / 365.0) + year;
865 
866  /* Determin declination */
867  geoMagInit();
868  geomg1(hgt,lat,lon,time,&dec,&dip,&ti,&gv);
869 
870  return (dec);
871 
872 }
873 /*************************************************************************/
874 // Short GeoMag routine
875 /*! \brief This routine computes magnetic declination for now time
876 
877 \param double lat : in: latitude in decimal degrees
878 \param double lon : in: longitude in decimal degrees
879 \param double hgt : in: altitude in meters
880 \retval double dec :out: declination in degrees
881 */
882 double DLL_CALLCONV geoMagGetDecNow(double lat, double lon, double hgt)
883 {
884  struct tm *newtime;
885  time_t aclock;
886  double dec;
887 
888  /* Compute the Magnetic Declination */
889  time( &aclock ); /* Get time in seconds */
890  newtime = localtime( &aclock ); /* Convert time to struct tm form */
891 
892  if(newtime->tm_year < 200 ) newtime->tm_year += 1900;
893 
894  // printf("%d / %d / %d \n",newtime->tm_mon+1, newtime->tm_mday, newtime->tm_year);
895  geoMagGetDec(lat,lon,hgt,
896  newtime->tm_mon+1, newtime->tm_mday, newtime->tm_year, &dec);
897 
898  return (dec);
899 
900 }
901 
902 /* Fill the geo struct with the magnetic declination using the computer's time/date */
903 int DLL_CALLCONV geoMagFillDec(GEO_LOCATION *l, double *dec)
904 {
905  struct tm *newtime;
906  time_t aclock;
907 
908  /* Compute the Magnetic Declination */
909  time( &aclock ); /* Get time in seconds */
910  newtime = localtime( &aclock ); /* Convert time to struct tm form */
911 
912  if(newtime->tm_year < 200 ) newtime->tm_year += 1900;
913 
914  // printf("%d / %d / %d \n",newtime->tm_mon+1, newtime->tm_mday, newtime->tm_year);
915  geoMagGetDec(l->lat,l->lon,l->hgt,
916  newtime->tm_mon+1, newtime->tm_mday, newtime->tm_year, dec);
917  l->Declination = *dec;
918 
919  return (GEO_OK);
920 }
921 
922 // Main test routine
923 #ifdef GEOMAG_MAIN
924 void main(void)
925 {
926  int i;
927  double dlat, dlon;
928  double ati, adec, adip;
929  double alt, gv;
930  double time, dec, dip, ti;
931  double x,y,z,h;
932  double ax,ay,az,ah;
933 
934 
935  dlat = 40.0;
936  dlon = -105.0;
937  alt = 1600.0;
938  time = EPOCH;
939 
940  // use the geoMag routine to get all of the specific values
941  geoMag(alt,dlat,dlon,time,&dec,&dip,&ti,&gv,&adec,&adip,&ati,
942  &x, &y, &z, &h, &ax, &ay, &az, &ah);
943 
944 
945  printf("\nLATITUDE: %7.2f Deg",dlat);
946  printf( " LONGITUDE: %7.2f Deg",dlon);
947  printf( " ALTITUDE: %.2f M",alt);
948  printf("\n DATE : %5.1f\n",time);
949  printf("\n\t\t\t OUTPUT\n\t\t\t ------");
950  printf("\n\nMAIN FIELD COMPONENTS\t\t\t ANNUAL CHANGE");
951  printf("\n---------------------\t\t\t -------------\n");
952  printf("\n TI = %-7.0f nT\t\t TI = %-6.0f nT/yr",ti,ati);
953  printf("\n HI = %-7.0f nT\t\t HI = %-6.0f nT/yr",h,ah);
954  printf("\n X = %-7.0f nT\t\t X = %-6.0f nT/yr ",x,ax);
955  printf("\n Y = %-7.0f nT\t\t Y = %-6.0f nT/yr ",y,ay);
956  printf("\n Z = %-7.0f nT\t\t Z = %-6.0f nT/yr ",z,az);
957  printf("\n DEC = %-7.2f DEG\t\t DEC = %-6.2f MIN/yr [%-6.2f deg/yr] ",dec,adec,adec/60.0);
958  printf("\n DIP = %-7.2f DEG\t\t DIP = %-6.2f MIN/yr\n",dip,adip);
959 
960 
961  // Print the monthly declination values
962  printf("----------------------------------------------------------\n");
963  printf(" Yearly Declination Values \n");
964  printf("----------------------------------------------------------\n");
965  for(i=0;i<=4;i++)
966  {
967  geoMagGetDec(dlat,dlon,alt,1,1,2000+i, &dec);
968  printf("%04d : %f deg\n",i+2000,dec);
969  }
970 
971  // Print the simple declination values for a location and time
972  printf("\n\n----------------------------------------------------------\n");
973  printf(" Lat, Lon, Hgt, and time (MM/DD/YY) declination example \n");
974  printf("----------------------------------------------------------\n");
975 
976  geoMagGetDec(dlat, dlon, alt, 6, 1, 2000, &dec);
977  printf(" Declination for %7.1f latitude, %5.1f longitude, %5.1f altitude \n is %4.2f degrees on 6/1/2000\n",
978  dlat, dlon, alt, dec);
979 
980 }
981 
982 #endif
983 /*************************************************************************/
  • Generated by doxygen 1.8.4