source: Sophya/trunk/SophyaExt/XAstroPack/magfield.c@ 1679

Last change on this file since 1679 was 1672, checked in by cmv, 24 years ago

modeles de champ magnetiques terrestres cmv 06/10/01
File size: 33.5 KB
RevLine 
[1672]1/* module magfield.c */
2
3/* Module to calculate magnetic variation and field given position,
4** altitude, and date
5** Implements the NIMA (formerly DMA) WMM and IGRF models
6**
7** http://www.nima.mil/GandG/ngdc-wmm2000.html
8** For WMM2000 coefficients:
9** ftp://ftp.ngdc.noaa.gov/Solid_Earth/Mainfld_Mag/DoD_Model/wmm.cof
10** For IGRF/DGRF coefficients:
11** http://swdcdb.kugi.kyoto-u.ac.jp/igrf/coef/igrfall.d
12**
13** Copyright (C) 2000 Edward A Williams <Ed_Williams@compuserve.com>
14**
15** The routine uses a spherical harmonic expansion of the magnetic
16** potential up to twelfth order, together with its time variation, as
17** described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
18** Academic Press (London 1987)". The program first converts geodetic
19** coordinates (lat/long on elliptic earth and altitude) to spherical
20** geocentric (spherical lat/long and radius) coordinates. Using this,
21** the spherical (B_r, B_theta, B_phi) magnetic field components are
22** computed from the model. These are finally referred to surface (X, Y,
23** Z) coordinates.
24**
25** Fields are accurate to better than 200nT, variation and dip to
26** better than 0.5 degrees, with the exception of the declination near
27** the magnetic poles (where it is ill-defined) where the error may reach
28** 4 degrees or more.
29**
30** Variation is undefined at both the geographic and
31** magnetic poles, even though the field itself is well-behaved. To
32** avoid the routine blowing up, latitude entries corresponding to
33** the geographic poles are slightly offset. At the magnetic poles,
34** the routine returns zero variation.
35**
36** HISTORY
37** Adapted from EAW Excel 3.0 version 3/27/94 EAW
38** Recoded in C++ by Starry Chan
39** WMM95 added and rearranged in ANSI-C EAW 7/9/95
40** Put shell around program and made Borland & GCC compatible EAW 11/22/95
41** IGRF95 added 2/96 EAW
42** WMM2000 IGR2000 added 2/00 EAW
43** Released under GPL 3/26/00 EAW
44** Adaptions and modifications for the SimGear project 3/27/2000 CLO
45** Removed all pow() calls and made static roots[][] arrays to
46** save many sqrt() calls on subsequent invocations
47** 3/28/2000 Norman Vine -- nhv@yahoo.com
48** Put in some bullet-proofing to handle magnetic and geographic poles.
49** 3/28/2000 EAW
50** Added missing comment close, the lack of which caused the altitude
51** correction to be omitted.
52** 01/31/01 Jim Seymour (jseymour@LinxNet.com)
53**
54** This program is free software; you can redistribute it and/or
55** modify it under the terms of the GNU General Public License as
56** published by the Free Software Foundation; either version 2 of the
57** License, or (at your option) any later version.
58**
59** This program is distributed in the hope that it will be useful, but
60** WITHOUT ANY WARRANTY; without even the implied warranty of
61** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
62** General Public License for more details.
63**
64** You should have received a copy of the GNU General Public License
65** along with this program; if not, write to the Free Software
66** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
67**
68*/
69
70#include <stdio.h>
71#include <stdlib.h>
72#include <math.h>
73#include "magfield.h"
74
75#define max(a,b) (((a) > (b)) ? (a) : (b))
76
77static const double pi = 3.14159265358979;
78static const double a = 6378.16; /* major radius (km) IAU66 ellipsoid */
79static const double f = 1.0 / 298.25; /* inverse flattening IAU66 ellipsoid */
80static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
81 /* minor radius b=a*(1-f) */
82static const double r_0 = 6371.2; /* "mean radius" for spherical harmonic expansion */
83
84/*IGRF90 constants */
85static const double gnm_igrf90[13][13] =
86 {
87 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
88 {-29775.4, -1851, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
89 {-2135.8, 3058.2, 1693.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
90 {1314.6, -2240.2, 1245.6, 806.5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
91 {938.9, 782.3, 323.9, -422.7, 141.7, 0, 0, 0, 0, 0, 0, 0, 0},
92 {-211, 352.5, 243.8, -110.8, -165.6, -37, 0, 0, 0, 0, 0, 0, 0},
93 {60.7, 63.9, 60.4, -177.5, 2, 16.7, -96.3, 0, 0, 0, 0, 0, 0},
94 {76.6, -64.2, 3.7, 27.5, 0.9, 5.7, 9.8, -0.5, 0, 0, 0, 0, 0},
95 {22.4, 5.1, -0.9, -10.8, -12.4, 3.8, 3.8, 2.6, -6, 0, 0, 0, 0},
96 {4.4, 9.9, 0.8, -12, 9.3, -3.9, -1.4, 7.3, 1.5, -5.5, 0, 0, 0},
97 {-3.6, -3.9, 2.4, -5.3, -2.4, 4.4, 3, 1.2, 2.2, 2.9, 0, 0, 0},
98 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
99 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
100static const double hnm_igrf90[13][13] =
101 {
102 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
103 {0, 5410.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
104 {0, -2277.7, -380, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
105 {0, -286.5, 293.3, -348.5, 0, 0, 0, 0, 0, 0, 0, 0, 0},
106 {0, 248.1, -239.5, 87, -299.4, 0, 0, 0, 0, 0, 0, 0, 0},
107 {0, 47.2, 153.5, -154.4, -69.2, 97.7, 0, 0, 0, 0, 0, 0, 0},
108 {0, -15.8, 82.7, 68.3, -52.5, 1.8, 26.9, 0, 0, 0, 0, 0, 0},
109 {0, -81.1, -27.3, 0.6, 20.4, 16.4, -22.6, -5, 0, 0, 0, 0, 0},
110 {0, 9.7, -19.9, 7.1, -22.1, 11.9, 11, -16, -10.7, 0, 0, 0, 0},
111 {0, -20.8, 15.4, 9.5, -5.7, -6.4, 8.6, 9.1, -6.6, 1.9, 0, 0, 0},
112 {0, 1.3, 0.4, 3.1, 5.6, -4.2, -0.5, -1.5, 3.8, -0.5, -6.2, 0, 0},
113 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
114 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
115static const double gtnm_igrf90[13][13] =
116 {
117 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
118 {18, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
119 {-12.9, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
120 {3.3, -6.7, 0, -5.9, 0, 0, 0, 0, 0, 0, 0, 0, 0},
121 {0.5, 0.6, -7, 0.5, -5.5, 0, 0, 0, 0, 0, 0, 0, 0},
122 {0.6, -0.1, -1.6, -3.1, 0, 2.3, 0, 0, 0, 0, 0, 0, 0},
123 {1.3, -0.2, 1.8, 1.3, -0.2, 0.1, 1.2, 0, 0, 0, 0, 0, 0},
124 {0.6, -0.5, -0.3, 0.6, 1.6, 0.2, 0.2, 0.3, 0, 0, 0, 0, 0},
125 {0.2, -0.7, -0.2, 0.1, -1.1, 0, 0, -0.5, -0.6, 0, 0, 0, 0},
126 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
127 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
128 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
129 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
130static const double htnm_igrf90[13][13] =
131 {
132 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
133 {0, -16.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
134 {0, -15.8, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
135 {0, 4.4, 1.6, -10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0},
136 {0, 2.6, 1.8, 3.1, -1.4, 0, 0, 0, 0, 0, 0, 0, 0},
137 {0, -0.1, 0.5, 0.4, 1.7, 0.4, 0, 0, 0, 0, 0, 0, 0},
138 {0, 0.2, -1.3, 0, -0.9, 0.5, 1.2, 0, 0, 0, 0, 0, 0},
139 {0, 0.6, 0.2, 0.8, -0.5, -0.2, 0, 0, 0, 0, 0, 0, 0},
140 {0, 0.5, -0.2, 0.3, 0.3, 0.4, -0.5, -0.3, 0.6, 0, 0, 0, 0},
141 {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},
143 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
144 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
145/* IGRF95 */
146static const double gnm_igrf95[13][13] =
147{
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 {-29682.0, -1789.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
150 {-2197.0, 3074.0, 1685.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
151 {1329.0, -2268.0, 1249.0, 769.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
152 {941.0, 782.0, 291.0, -421.0, 116.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
153 {-210.0, 352.0, 237.0, -122.0, -167.0, -26.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0},
154 {66.0, 64.0, 65.0, -172.0, 2.0, 17.0, -94.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
155 {78.0, -67.0, 1.0, 29.0, 4.0, 8.0, 10.0, -2.0, 0.0,0.0, 0.0, 0.0, 0.0},
156 {24.0, 4.0, -1.0, -9.0, -14.0, 4.0, 5.0, 0.0, -7.0, 0.0,0.0, 0.0, 0.0},
157 {4.0, 9.0, 1.0, -12.0, 9.0, -4.0, -2.0, 7.0, 0.0, -6.0, 0.0, 0.0, 0.0},
158 {-3.0, -4.0, 2.0, -5.0, -2.0, 4.0, 3.0, 1.0, 3.0, 3.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}
161 };
162
163static const double hnm_igrf95[13][13]=
164{
165 {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},
166 {0.0, 5318.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
167 {0.0, -2356.0, -425.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
168 {0.0, -263.0, 302.0, -406.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
169 {0.0, 262.0, -232.0, 98.0, -301.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
170 {0.0, 44.0, 157.0, -152.0, -64.0, 99.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
171 {0.0, -16.0, 77.0, 67.0, -57.0, 4.0, 28.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
172 {0.0, -77.0, -25.0, 3.0, 22.0, 16.0, -23.0, -3.0, 0.0,0.0, 0.0, 0.0, 0.0},
173 {0.0, 12.0, -20.0, 7.0, -21.0, 12.0, 10.0, -17.0, -10.0,0.0, 0.0, 0.0, 0.0},
174 {0.0, -19.0, 15.0, 11.0, -7.0, -7.0, 9.0, 7.0, -8.0, 1.0, 0.0, 0.0, 0.0},
175 {0.0, 2.0, 1.0, 3.0, 6.0, -4.0, 0.0, -2.0, 3.0, -1.0, -6.0, 0.0, 0.0},
176 {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},
177 {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}
178};
179
180static const double gtnm_igrf95[13][13]=
181{
182 {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},
183 {17.6, 13.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
184 {-13.2, 3.7, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
185 {1.5, -6.4, -0.2, -8.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
186 {0.8, 0.9, -6.9, 0.5, -4.6, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
187 {0.8, 0.1, -1.5, -2.0, -0.1, 2.3, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
188 {0.5, -0.4, 0.6, 1.9, -0.2, -0.2, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
189 {-0.2, -0.8, -0.6, 0.6, 1.2, 0.1, 0.2, -0.6, 0.0,0.0, 0.0, 0.0, 0.0},
190 {0.3, -0.2, 0.1, 0.4, -1.1, 0.3, 0.2, -0.9, -0.3, 0.0,0.0, 0.0, 0.0},
191 {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},
192 {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},
193 {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},
194 {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}
195 };
196
197static const double htnm_igrf95[13][13]=
198{
199 {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},
200 {0.0, -18.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
201 {0.0, -15.0, -8.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
202 {0.0, 4.1, 2.2, -12.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
203 {0.0, 1.8, 1.2, 2.7, -1.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
204 {0.0, 0.2, 1.2, 0.3, 1.8, 0.9, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
205 {0.0, 0.3, -1.6, -0.2, -0.9, 1.0, 2.2, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
206 {0.0, 0.8, 0.2, 0.6, -0.4, 0.0, -0.3, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
207 {0.0, 0.4, -0.2, 0.2, 0.7, 0.0, -1.2, -0.7, -0.6, 0.0,0.0, 0.0, 0.0},
208 {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},
209 {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},
210 {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},
211 {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}
212};
213
214/*WMM85 constants */
215static const double gnm_wmm85[13][13] =
216 {
217 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
218 {-29879.8, -1903.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
219 {-2070.6, 3040.8, 1696.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
220 {1303.9, -2203, 1241.7, 839.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
221 {933.8, 781.8, 359, -424.5, 164.5, 0, 0, 0, 0, 0, 0, 0, 0},
222 {-216.4, 353, 254.3, -93.7, -157.5, -45.2, 0, 0, 0, 0, 0, 0, 0},
223 {53.2, 63.8, 51.3, -188.4, 3.3, 20.3, -101.7, 0, 0, 0, 0, 0, 0},
224 {76.9, -60.7, 0.7, 25.4, -8.1, 6.9, 7, -4.4, 0, 0, 0, 0, 0},
225 {18.4, 5.1, 1.2, -12, -9.1, 0.1, 4.7, 6.5, -9.5, 0, 0, 0, 0},
226 {5.7, 10.9, 0.9, -12.2, 9.5, -3.3, -1, 6.5, 1.5, -4.8, 0, 0, 0},
227 {-3.4, -4.7, 2.5, -5.5, -2.1, 4.6, 3.2, 0.6, 1.9, 2.8, -0.2, 0, 0},
228 {2.3, -0.8, -2, 2.1, 0.2, -0.4, -0.4, 1.6, 1.5, -0.7, 2.3, 3.5, 0},
229 {-1.8, 0, 0.1, -0.3, 0.5, 0.5, -0.6, -0.4, 0, -0.5, 0, 0.7, -0.2}};
230static const double hnm_wmm85[13][13] =
231 {
232 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
233 {0, 5490.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
234 {0, -2189.1, -312, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
235 {0, -310.3, 282.6, -299.2, 0, 0, 0, 0, 0, 0, 0, 0, 0},
236 {0, 227.2, -246.7, 72.5, -299.1, 0, 0, 0, 0, 0, 0, 0, 0},
237 {0, 43.4, 148.2, -154.8, -71.8, 91.5, 0, 0, 0, 0, 0, 0, 0},
238 {0, -12.3, 87.9, 67.8, -51.1, -4, 20.8, 0, 0, 0, 0, 0, 0},
239 {0, -80.1, -25.9, -0.9, 21.6, 18.5, -20, -7.7, 0, 0, 0, 0, 0},
240 {0, 3.8, -20.2, 5, -24.2, 12.2, 7.6, -16.3, -10.9, 0, 0, 0, 0},
241 {0, -20.8, 15.8, 9, -5, -6.4, 9.1, 9.9, -5.8, 2.3, 0, 0, 0},
242 {0, 1.2, 0.4, 2.5, 5.6, -4.4, -0.5, -1.6, 3.7, -0.5, -6.1, 0, 0},
243 {0, 1.3, 2, -1.1, -2.8, 0.7, -0.1, -2.4, -0.4, -1.5, -1.5, 0.7, 0},
244 {0, 0.3, 0.6, 2.5, -1.7, 0.3, 0.2, -0.1, 0.1, 0.1, -1.4, 0.4, 0.7}};
245static const double gtnm_wmm85[13][13] =
246 {
247 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
248 {21.9, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
249 {-11.2, 1.8, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
250 {8.3, -2, -0.6, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
251 {-1.2, 0.1, -9.7, -1.7, -9.3, 0, 0, 0, 0, 0, 0, 0, 0},
252 {1.4, -0.5, -1.2, -2.2, 0.9, 0, 0, 0, 0, 0, 0, 0, 0},
253 {3.1, 0, 1.8, -0.2, -0.4, 2.4, 1.8, 0, 0, 0, 0, 0, 0},
254 {-0.1, -0.8, -1.2, 1.1, 0, 0.6, -1.8, -1.2, 0, 0, 0, 0, 0},
255 {0.2, 0, 0.7, 0.1, 0.2, -0.3, -0.1, 0.2, -2.2, 0, 0, 0, 0},
256 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
257 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
258 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
259 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
260static const double htnm_wmm85[13][13] =
261 {
262 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
263 {0, -31.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
264 {0, -9.7, -19.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
265 {0, 6.1, 1.3, -13, 0, 0, 0, 0, 0, 0, 0, 0, 0},
266 {0, 1.3, 3.6, 2.5, 0.6, 0, 0, 0, 0, 0, 0, 0, 0},
267 {0, -0.9, 0.6, 0.3, 2.4, -1.4, 0, 0, 0, 0, 0, 0, 0},
268 {0, 0.7, -2.1, -1.4, -4.3, -0.7, 0, 0, 0, 0, 0, 0, 0},
269 {0, 0, 1.2, 2, 2.6, 0.9, 0.8, 0.4, 0, 0, 0, 0, 0},
270 {0, -0.6, -1.5, 0.1, -1.1, 0.4, -2, 0.9, 1.5, 0, 0, 0, 0},
271 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
272 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
273 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
274 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
275
276/* wmm90 constants */
277static const double gnm_wmm90[13][13] =
278 {
279 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
280 {-29780.5, -1851.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
281 {-2134.3, 3062.2, 1691.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
282 {1312.9, -2244.7, 1246.8, 808.6, 0, 0, 0, 0, 0, 0, 0, 0, 0},
283 {933.5, 784.9, 323.5, -421.7, 139.2, 0, 0, 0, 0, 0, 0, 0, 0},
284 {-208.3, 352.2, 246.5, -110.8, -162.3, -37.2, 0, 0, 0, 0, 0, 0, 0},
285 {59, 63.7, 60, -181.3, 0.4, 15.4, -96, 0, 0, 0, 0, 0, 0},
286 {76.1, -62.1, 1.3, 30.2, 4.7, 7.9, 10.1, 1.9, 0, 0, 0, 0, 0},
287 {22.9, 2.3, -1.2, -11.7, -17.5, 2.2, 5.7, 3, -7, 0, 0, 0, 0},
288 {3.6, 9.5, -0.9, -10.7, 10.7, -3.2, -1.4, 6.3, 0.8, -5.5, 0, 0, 0},
289 {-3.3, -2.6, 4.5, -5.6, -3.6, 3.9, 3.2, 1.7, 3, 3.7, 0.7, 0, 0},
290 {1.3, -1.4, -2.5, 3.2, 0.2, -1.1, 0.3, -0.3, 0.9, -1.1, 2.4, 3, 0},
291 {-1.3, 0.1, 0.5, 0.7, 0.4, -0.2, -1.1, 0.9, -0.6, 0.8, 0.2, 0.4, 0.2}};
292static const double hnm_wmm90[13][13] =
293 {
294 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
295 {0, 5407.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
296 {0, -2278.3, -384.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
297 {0, -284.9, 291.7, -352.4, 0, 0, 0, 0, 0, 0, 0, 0, 0},
298 {0, 249.4, -232.7, 91.3, -296.5, 0, 0, 0, 0, 0, 0, 0, 0},
299 {0, 40.8, 148.7, -154.6, -67.6, 97.4, 0, 0, 0, 0, 0, 0, 0},
300 {0, -14.7, 82.2, 70, -56.2, -1.4, 24.6, 0, 0, 0, 0, 0, 0},
301 {0, -78.6, -26.7, 0.1, 19.9, 17.9, -21.5, -6.8, 0, 0, 0, 0, 0},
302 {0, 9.7, -19.3, 6.6, -20.1, 13.4, 9.8, -19, -9.1, 0, 0, 0, 0},
303 {0, -21.9, 14.3, 9.5, -6.7, -6.4, 9.1, 8.9, -8, 2.1, 0, 0, 0},
304 {0, 2.6, 1.2, 2.6, 5.7, -4, -0.4, -1.7, 3.8, -0.8, -6.5, 0, 0},
305 {0, 0, 1, -1.6, -2.2, 1.1, -0.7, -1.7, -1.5, -1.3, -1.1, 0.6, 0},
306 {0, 0.7, 0.7, 1.3, -1.5, 0.3, 0.2, -1.1, 1.2, -0.2, -1.3, 0.6, 0.6}};
307static const double gtnm_wmm90[13][13] =
308 {
309 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
310 {16, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
311 {-11.7, 3.7, 1.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
312 {2.1, -7.6, 0, -5.8, 0, 0, 0, 0, 0, 0, 0, 0, 0},
313 {-0.8, 1, -7.4, 0.8, -6.4, 0, 0, 0, 0, 0, 0, 0, 0},
314 {1.7, 0, 0, -2.7, 0, 3, 0, 0, 0, 0, 0, 0, 0},
315 {0.8, 0, 1.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
316 {0.5, 0, -0.9, 1.5, 2.7, -1, 0, 0, 0, 0, 0, 0, 0},
317 {0, -1.1, 0, 0, -2.1, 0, 1, 0, 0, 0, 0, 0, 0},
318 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
319 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
320 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
321 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
322static const double htnm_wmm90[13][13] =
323 {
324 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
325 {0, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
326 {0, -12.8, -14.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
327 {0, 3.1, 0.8, -11.3, 0, 0, 0, 0, 0, 0, 0, 0, 0},
328 {0, 3.3, 3.7, 2.8, 0, 0, 0, 0, 0, 0, 0, 0, 0},
329 {0, 0, -2.1, 1.2, 1.2, 0.6, 0, 0, 0, 0, 0, 0, 0},
330 {0, -0.6, -0.6, 0, -2.3, 0, 0, 0, 0, 0, 0, 0, 0},
331 {0, 0.6, 0.8, 0, 0, 0, 0.4, 0, 0, 0, 0, 0, 0},
332 {0, 0.4, -0.8, 0.5, 0.3, 0.5, 0, -0.7, 0, 0, 0, 0, 0},
333 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
334 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
335 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
336 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
337/* wmm95 constants */
338static const double gnm_wmm95[13][13] =
339 {
340 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
341 {-29682.1,-1782.2,0,0,0,0,0,0,0,0,0,0,0},
342 {-2194.7,3078.6,1685.7,0,0,0,0,0,0,0,0,0,0},
343 {1318.8,-2273.6,1246.9,766.3,0,0,0,0,0,0,0,0,0},
344 {940,782.9,290.9,-418.9,113.8,0,0,0,0,0,0,0,0},
345 {-209.5,354,238.2,-122.1,-162.8,-23.3,0,0,0,0,0,0,0},
346 {68.5,65.6,64.1,-169.1,-0.5,16.5,-91,0,0,0,0,0,0},
347 {78,-68.1,0.1,29.6,6,8.7,9.2,-2.4,0,0,0,0,0},
348 {24.7,3.4,-1.5,-9.6,-16.5,2.6,3.6,-4.9,-8.5,0,0,0,0},
349 {2.9,7.5,0.4,-10.3,9.7,-2.3,-2.4,6.8,-0.5,-6.5,0,0,0},
350 {-2.9,-3.3,2.8,-4.3,-3.1,2.4,2.8,0.7,4.1,3.6,0.6,0,0},
351 {1.7,-1.6,-3.6,1.2,-0.6,0.1,-0.7,-0.8,1.3,-0.3,2.2,4.2,0},
352 {-1.8,0.9,-0.1,-0.5,0.8,0.2,0.5,0.4,-0.4,0.3,0.2,0.4,0.6}};
353static const double hnm_wmm95[13][13] =
354 {
355 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
356 {0,5315.6,0,0,0,0,0,0,0,0,0,0,0},
357 {0,-2359.1,-418.6,0,0,0,0,0,0,0,0,0,0},
358 {0,-261.1,301,-416.5,0,0,0,0,0,0,0,0,0},
359 {0,259.4,-230.9,99.8,-306.1,0,0,0,0,0,0,0,0},
360 {0,43.7,157.6,-150.1,-59.2,104.4,0,0,0,0,0,0,0},
361 {0,-15.2,74.3,69.4,-55.3,3,33.3,0,0,0,0,0,0},
362 {0,-76.1,-24.5,1.6,20,16.5,-23.6,-6.8,0,0,0,0,0},
363 {0,14.9,-19.5,6.3,-20.4,12.2,7,-19,-8.8,0,0,0,0},
364 {0,-19.8,14.6,10.9,-7.5,-6.8,9.3,7.7,-8.1,2.6,0,0,0},
365 {0,3.2,1.7,2.9,5.6,-3.4,-0.7,-2.9,2.3,-1.6,-6.6,0,0},
366 {0,0.3,1,-3.6,-1.4,1.9,0.2,-1.3,-2.4,-0.6,-2.2,1.3,0},
367 {0,0.3,1.4,0.8,-3,0.7,0.5,-0.8,0.6,0.1,-1.3,-0.4,0.9}};
368static const double gtnm_wmm95[13][13] =
369 {
370 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
371 {17.6,13.2,0,0,0,0,0,0,0,0,0,0,0},
372 {-13.7,4,-0.3,0,0,0,0,0,0,0,0,0,0},
373 {0.8,-6.6,-0.5,-8.5,0,0,0,0,0,0,0,0,0},
374 {1.2,1.1,-6.8,0.3,-4.5,0,0,0,0,0,0,0,0},
375 {0.9,0.5,-1.4,-1.7,0,2.1,0,0,0,0,0,0,0},
376 {0.4,-0.3,0.3,2.1,0,-0.4,-0.4,0,0,0,0,0,0},
377 {-0.3,-1.1,-0.5,0.5,1.3,0.1,0,-0.9,0,0,0,0,0},
378 {0.1,0,0.4,0.3,-1.3,0.5,0.4,-0.9,0.1,0,0,0,0},
379 {0,0,0,0,0,0,0,0,0,0,0,0,0},
380 {0,0,0,0,0,0,0,0,0,0,0,0,0},
381 {0,0,0,0,0,0,0,0,0,0,0,0,0},
382 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
383static const double htnm_wmm95[13][13] =
384 {
385 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
386 {0,-18,0,0,0,0,0,0,0,0,0,0,0},
387 {0,-14.6,-7.2,0,0,0,0,0,0,0,0,0,0},
388 {0,4,2.2,-12.6,0,0,0,0,0,0,0,0,0},
389 {0,1.3,1,2.5,-1.2,0,0,0,0,0,0,0,0},
390 {0,0.5,1.5,0.6,1.7,0.6,0,0,0,0,0,0,0},
391 {0,0.7,-1.5,-0.5,-0.7,1.1,2.6,0,0,0,0,0,0},
392 {0,0.3,0,0.7,-0.6,0.1,-0.6,-0.4,0,0,0,0,0},
393 {0,0.4,-0.3,0.1,0.8,-0.1,-1.3,-0.9,-1.1,0,0,0,0},
394 {0,0,0,0,0,0,0,0,0,0,0,0,0},
395 {0,0,0,0,0,0,0,0,0,0,0,0,0},
396 {0,0,0,0,0,0,0,0,0,0,0,0,0},
397 {0,0,0,0,0,0,0,0,0,0,0,0,0}};
398
399static const double gnm_wmm2000[13][13] =
400 {
401 {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},
402 {-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},
403 {-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},
404 {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},
405 {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},
406 {-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},
407 {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},
408 {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},
409 {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},
410 {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},
411 {-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},
412 {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},
413 {-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},
414 };
415
416static const double hnm_wmm2000[13][13]=
417 {
418 {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},
419 {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},
420 {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},
421 {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},
422 {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},
423 {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},
424 {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},
425 {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},
426 {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},
427 {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},
428 {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},
429 {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},
430 {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},
431 };
432
433static const double gtnm_wmm2000[13][13]=
434 {
435 {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},
436 {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},
437 {-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},
438 {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},
439 {-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},
440 {-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},
441 {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},
442 {-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},
443 {-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},
444 {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},
445 {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},
446 {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},
447 {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},
448 };
449
450static const double htnm_wmm2000[13][13]=
451 {
452 {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},
453 {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},
454 {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},
455 {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},
456 {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},
457 {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},
458 {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},
459 {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},
460 {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},
461 {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},
462 {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},
463 {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},
464 {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},
465 };
466
467static const double gnm_igrf2000[13][13] =
468 {
469 {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},
470 {-29615.0, -1728.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
471 {-2267.0, 3072.0, 1672.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
472 {1341.0, -2290.0, 1253.0, 715.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
473 {935.0, 787.0, 251.0, -405.0, 110.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
474 {-217.0, 351.0, 222.0, -131.0, -169.0, -12.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
475 {72.0, 68.0, 74.0, -161.0, -5.0, 17.0, -91.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
476 {79.0, -74.0, 0.0, 33.0, 9.0, 7.0, 8.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
477 {25.0, 6.0, -9.0, -8.0, -17.0, 9.0, 7.0, -8.0, -7.0, 0.0, 0.0, 0.0, 0.0},
478 {5.0, 9.0, 3.0, -8.0, 6.0, -9.0, -2.0, 9.0, -4.0, -8.0, 0.0, 0.0, 0.0},
479 {-2.0, -6.0, 2.0, -3.0, 0.0, 4.0, 1.0, 2.0, 4.0, 0.0, -1.0, 0.0, 0.0},
480 {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},
481 {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},
482 };
483
484static const double hnm_igrf2000[13][13]=
485 {
486 {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},
487 {0.0, 5186.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
488 {0.0, -2478.0, -458.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
489 {0.0, -227.0, 296.0, -492.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
490 {0.0, 272.0, -232.0, 119.0, -304.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
491 {0.0, 44.0, 172.0, -134.0, -40.0, 107.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
492 {0.0, -17.0, 64.0, 65.0, -61.0, 1.0, 44.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
493 {0.0, -65.0, -24.0, 6.0, 24.0, 15.0, -25.0, -6.0, 0.0, 0.0, 0.0, 0.0, 0.0},
494 {0.0, 12.0, -22.0, 8.0, -21.0, 15.0, 9.0, -16.0, -3.0, 0.0, 0.0, 0.0, 0.0},
495 {0.0, -20.0, 13.0, 12.0, -6.0, -8.0, 9.0, 4.0, -8.0, 5.0, 0.0, 0.0, 0.0},
496 {0.0, 1.0, 0.0, 4.0, 5.0, -6.0, -1.0, -3.0, 0.0, -2.0, -8.0, 0.0, 0.0},
497 {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},
498 {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},
499 };
500
501static const double gtnm_igrf2000[13][13]=
502 {
503 {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},
504 {14.6, 10.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
505 {-12.4, 1.1, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
506 {0.7, -5.4, 0.9, -7.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
507 {-1.3, 1.6, -7.3, 2.9, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
508 {0.0, -0.7, -2.1, -2.8, -0.8, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
509 {1.0, -0.4, 0.9, 2.0, -0.6, -0.3, 1.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
510 {-0.4, -0.4, -0.3, 1.1, 1.1, -0.2, 0.6, -0.9, 0.0, 0.0, 0.0, 0.0, 0.0},
511 {-0.3, 0.2, -0.3, 0.4, -1.0, 0.3, -0.5, -0.7, -0.4, 0.0, 0.0, 0.0, 0.0},
512 {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},
513 {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},
514 {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},
515 {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},
516 };
517
518static const double htnm_igrf2000[13][13]=
519 {
520 {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},
521 {0.0, -22.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
522 {0.0, -20.6, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
523 {0.0, 6.0, -0.1, -14.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
524 {0.0, 2.1, 1.3, 5.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
525 {0.0, -0.1, 0.6, 1.7, 1.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
526 {0.0, -0.2, -1.4, 0.0, -0.8, 0.0, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
527 {0.0, 1.1, 0.0, 0.3, -0.1, -0.6, -0.7, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
528 {0.0, 0.1, 0.0, 0.0, 0.3, 0.6, -0.4, 0.3, 0.7, 0.0, 0.0, 0.0, 0.0},
529 {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},
530 {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},
531 {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},
532 {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},
533 };
534
535
536
537static const int nmax = 12;
538
539double P[13][13];
540double DP[13][13];
541double gnm[13][13];
542double hnm[13][13];
543double sm[13];
544double cm[13];
545
546static double root[13];
547static double roots[13][13][2];
548
549/* Convert date to Julian day 1950-2049 */
550unsigned long int yymmdd_to_julian_days(int yy,int mm,int dd)
551{
552 unsigned long jd;
553
554 yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
555 jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
556 jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
557 jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
558
559 return(jd);
560}
561
562/* Convert degrees to radians */
563double deg_to_rad(double deg)
564{
565return deg*pi/180.;
566}
567
568/* Convert radians to degrees */
569double rad_to_deg(double rad)
570{
571return rad*180./pi;
572}
573
574
575
576
577/*
578 * return variation (in radians) given geodetic latitude (radians), longitude
579 * (radians) ,height (km) and (Julian) date
580 * model=1 is IGRF90, 2 is WMM85, 3 is WMM90, 4 is WMM95,
581 * 5 is IGRF95, 6 is WMM2000, 7 is IGRF2000
582 * N and E lat and long are positive, S and W negative
583*/
584
585double SGMagVar( double lat, double lon, double h, long dat, int model, double* field )
586{
587 /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
588 int n,m;
589 double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
590 double sinpsi, cospsi, inv_s;
591
592 static int been_here = 0;
593
594 double sinlat = sin(lat);
595 double coslat = cos(lat);
596
597 /* convert to geocentric */
598 sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
599 /* sr is effective radius */
600 theta = atan2(coslat * (h*sr + a*a), sinlat * (h*sr + b*b));
601
602 /* theta is geocentric co-latitude */
603
604 r = h*h + 2.0*h * sr +
605 (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
606 (a*a - (a*a - b*b) * sinlat*sinlat );
607
608 r = sqrt(r);
609
610 /* r is geocentric radial distance */
611 c = cos(theta);
612 s = sin(theta);
613 /* protect against zero divide at geographic poles */
614 inv_s = 1.0 / (s + (s == 0.)*1.0e-8);
615
616 /*zero out arrays */
617 for ( n = 0; n <= nmax; n++ ) {
618 for ( m = 0; m <= n; m++ ) {
619 P[n][m] = 0;
620 DP[n][m] = 0;
621 }
622 }
623
624 /* diagonal elements */
625 P[0][0] = 1;
626 P[1][1] = s;
627 DP[0][0] = 0;
628 DP[1][1] = c;
629 P[1][0] = c ;
630 DP[1][0] = -s;
631
632 /* these values will not change for subsequent function calls */
633 if( !been_here ) {
634 for ( n = 2; n <= nmax; n++ ) {
635 root[n] = sqrt((2.0*n-1) / (2.0*n));
636 }
637
638 for ( m = 0; m <= nmax; m++ ) {
639 double mm = m*m;
640 for ( n = max(m + 1, 2); n <= nmax; n++ ) {
641 roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
642 roots[m][n][1] = 1.0 / sqrt( n*n - mm);
643 }
644 }
645 been_here = 1;
646 }
647
648 for ( n=2; n <= nmax; n++ ) {
649 /* double root = sqrt((2.0*n-1) / (2.0*n)); */
650 P[n][n] = P[n-1][n-1] * s * root[n];
651 DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) * root[n];
652 }
653
654 /* lower triangle */
655 for ( m = 0; m <= nmax; m++ ) {
656 /* double mm = m*m; */
657 for ( n = max(m + 1, 2); n <= nmax; n++ ) {
658 /* double root1 = sqrt((n-1)*(n-1) - mm); */
659 /* double root2 = 1.0 / sqrt( n*n - mm); */
660 P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
661 P[n-2][m] * roots[m][n][0]) * roots[m][n][1];
662 DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
663 (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) * roots[m][n][1];
664 }
665 }
666
667 /* compute gnm, hnm at dat */
668 switch(model) {
669 case 1: /* IGRF90 */
670 yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
671 for (n=1;n<=nmax;n++)
672 for (m = 0;m<=nmax;m++) {
673 gnm[n][m] = gnm_igrf90[n][m] + yearfrac * gtnm_igrf90[n][m];
674 hnm[n][m] = hnm_igrf90[n][m] + yearfrac * htnm_igrf90[n][m];
675 }
676 break;
677
678 case 2: /* WMM85 */
679 yearfrac = (dat - yymmdd_to_julian_days(85,1,1)) / 365.25;
680 for (n=1;n<=nmax;n++)
681 for (m = 0;m<=nmax;m++) {
682 gnm[n][m] = gnm_wmm85[n][m] + yearfrac * gtnm_wmm85[n][m];
683 hnm[n][m] = hnm_wmm85[n][m] + yearfrac * htnm_wmm85[n][m];
684 }
685 break;
686
687 case 3: /* WMM90 */
688 yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
689 for (n=1;n<=nmax;n++)
690 for (m = 0;m<=nmax;m++) {
691 gnm[n][m] = gnm_wmm90[n][m] + yearfrac * gtnm_wmm90[n][m];
692 hnm[n][m] = hnm_wmm90[n][m] + yearfrac * htnm_wmm90[n][m];
693 }
694 break;
695
696 case 4: /* WMM95 */
697 yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
698 for (n=1;n<=nmax;n++)
699 for (m = 0;m<=nmax;m++) {
700 gnm[n][m] = gnm_wmm95[n][m] + yearfrac * gtnm_wmm95[n][m];
701 hnm[n][m] = hnm_wmm95[n][m] + yearfrac * htnm_wmm95[n][m];
702 }
703 break;
704 case 5: /* IGRF95 */
705 yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
706 for (n=1;n<=nmax;n++)
707 for (m = 0;m<=nmax;m++) {
708 gnm[n][m] = gnm_igrf95[n][m] + yearfrac * gtnm_igrf95[n][m];
709 hnm[n][m] = hnm_igrf95[n][m] + yearfrac * htnm_igrf95[n][m];
710 }
711 break;
712 case 6: /* WMM2000 */
713 yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
714 for (n=1;n<=nmax;n++)
715 for (m = 0;m<=nmax;m++) {
716 gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
717 hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
718 }
719 break;
720 case 7: /* IGRF2000 */
721 yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
722 for (n=1;n<=nmax;n++)
723 for (m = 0;m<=nmax;m++) {
724 gnm[n][m] = gnm_igrf2000[n][m] + yearfrac * gtnm_igrf2000[n][m];
725 hnm[n][m] = hnm_igrf2000[n][m] + yearfrac * htnm_igrf2000[n][m];
726 }
727 break;
728 }
729
730 /* compute sm (sin(m lon) and cm (cos(m lon)) */
731 for (m = 0;m<=nmax;m++) {
732 sm[m] = sin(m * lon);
733 cm[m] = cos(m * lon);
734 }
735
736 /* compute B fields */
737 B_r = 0.0;
738 B_theta = 0.0;
739 B_phi = 0.0;
740 fn_0 = r_0/r;
741 fn = fn_0 * fn_0;
742
743 for ( n = 1; n <= nmax; n++ ) {
744 double c1_n=0;
745 double c2_n=0;
746 double c3_n=0;
747 for ( m = 0; m <= n; m++ ) {
748 double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
749 c1_n += tmp * P[n][m];
750 c2_n += tmp * DP[n][m];
751 c3_n += m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
752 }
753 /* fn=pow(r_0/r,n+2.0); */
754 fn *= fn_0;
755 B_r += (n + 1) * c1_n * fn;
756 B_theta -= c2_n * fn;
757 B_phi += c3_n * fn * inv_s;
758 }
759
760
761
762 /* Find geodetic field components: */
763 psi = theta - (pi / 2.0 - lat);
764 sinpsi = sin(psi);
765 cospsi = cos(psi);
766 X = -B_theta * cospsi - B_r * sinpsi;
767 Y = B_phi;
768 Z = B_theta * sinpsi - B_r * cospsi;
769
770 field[0]=B_r;
771 field[1]=B_theta;
772 field[2]=B_phi;
773 field[3]=X;
774 field[4]=Y;
775 field[5]=Z; /* output fields */
776 /* find variation in radians */
777 /* return zero variation at magnetic pole X=Y=0. */
778 /* E is positive */
779 return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
780}
781
Note: See TracBrowser for help on using the repository browser.