java/com.sap.sailing.declination/src/com/sap/sailing/declination/impl/Geomagnetism.java
... ...
@@ -0,0 +1,448 @@
1
+package com.sap.sailing.declination.impl;
2
+
3
+/* License Statement from the NOAA
4
+* The WMM source code is in the public domain and not licensed or
5
+* under copyright. The information and software may be used freely
6
+* by the public. As required by 17 U.S.C. 403, third parties producing
7
+* copyrighted works consisting predominantly of the material produced
8
+* by U.S. government agencies must provide notice with such work(s)
9
+* identifying the U.S. Government material incorporated and stating
10
+* that such material is not subject to copyright protection.*/
11
+
12
+import java.util.GregorianCalendar;
13
+
14
+/**
15
+ * <p>
16
+ * Class to calculate magnetic declination, magnetic field strength, inclination etc. for any point on the earth.
17
+ * </p>
18
+ * <p>
19
+ * Adapted from the geomagc software and World Magnetic Model of the NOAA Satellite and Information Service, National
20
+ * Geophysical Data Center
21
+ * </p>
22
+ * http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
23
+ * <p>
24
+ * © Deep Pradhan, 2017
25
+ * </p>
26
+ */
27
+class Geomagnetism {
28
+
29
+ /** Initialise the instance without calculations */
30
+ Geomagnetism() {
31
+ // Initialize constants
32
+ maxord = MAX_DEG;
33
+ sp[0] = 0;
34
+ cp[0] = snorm[0] = pp[0] = 1;
35
+ dp[0][0] = 0;
36
+
37
+ c[0][0] = 0;
38
+ cd[0][0] = 0;
39
+
40
+ epoch = Double.parseDouble(WMM_COF[0].trim().split("\\s+")[0]);
41
+
42
+ String[] tokens;
43
+
44
+ double gnm, hnm, dgnm, dhnm;
45
+ for (int i = 1, m, n; i < WMM_COF.length; i++) {
46
+ tokens = WMM_COF[i].trim().split("\\s+");
47
+ n = Integer.parseInt(tokens[0]);
48
+ m = Integer.parseInt(tokens[1]);
49
+ gnm = Double.parseDouble(tokens[2]);
50
+ hnm = Double.parseDouble(tokens[3]);
51
+ dgnm = Double.parseDouble(tokens[4]);
52
+ dhnm = Double.parseDouble(tokens[5]);
53
+ if (m <= n) {
54
+ c[m][n] = gnm;
55
+ cd[m][n] = dgnm;
56
+ if (m != 0) {
57
+ c[n][m - 1] = hnm;
58
+ cd[n][m - 1] = dhnm;
59
+ }
60
+ }
61
+ }
62
+ // Convert schmidt normalized gauss coefficients to unnormalized
63
+ snorm[0] = 1;
64
+ double flnmj;
65
+ for (int j, n = 1; n <= maxord; n++) {
66
+ snorm[n] = snorm[n - 1] * (2 * n - 1) / n;
67
+ j = 2;
68
+ for (int m = 0, d1 = 1, d2 = (n - m + d1) / d1; d2 > 0; d2--, m += d1) {
69
+ k[m][n] = (double) (((n - 1) * (n - 1)) - (m * m)) / (double) ((2 * n - 1) * (2 * n - 3));
70
+ if (m > 0) {
71
+ flnmj = ((n - m + 1) * j) / (double) (n + m);
72
+ snorm[n + m * 13] = snorm[n + (m - 1) * 13] * Math.sqrt(flnmj);
73
+ j = 1;
74
+ c[n][m - 1] = snorm[n + m * 13] * c[n][m - 1];
75
+ cd[n][m - 1] = snorm[n + m * 13] * cd[n][m - 1];
76
+ }
77
+ c[m][n] = snorm[n + m * 13] * c[m][n];
78
+ cd[m][n] = snorm[n + m * 13] * cd[m][n];
79
+ }
80
+ fn[n] = (n + 1);
81
+ fm[n] = n;
82
+ }
83
+ k[1][1] = 0;
84
+ fm[0] = 0;
85
+ otime = oalt = olat = olon = -1000;
86
+ }
87
+
88
+ /**
89
+ * Initialise the instance and calculate for given location, altitude and date
90
+ *
91
+ * @param longitude
92
+ * Longitude in decimal degrees
93
+ * @param latitude
94
+ * Latitude in decimal degrees
95
+ * @param altitude
96
+ * Altitude in metres (with respect to WGS-1984 ellipsoid)
97
+ * @param calendar
98
+ * Calendar for date of calculation
99
+ */
100
+ Geomagnetism(double longitude, double latitude, double altitude, GregorianCalendar calendar) {
101
+ this();
102
+ calculate(longitude, latitude, altitude, calendar);
103
+ }
104
+
105
+ /**
106
+ * Initialise the instance and calculate for given location, altitude and current date
107
+ *
108
+ * @param longitude
109
+ * Longitude in decimal degrees
110
+ * @param latitude
111
+ * Latitude in decimal degrees
112
+ * @param altitude
113
+ * Altitude in metres (with respect to WGS-1984 ellipsoid)
114
+ */
115
+ Geomagnetism(double longitude, double latitude, double altitude) {
116
+ this();
117
+ calculate(longitude, latitude, altitude);
118
+ }
119
+
120
+ /**
121
+ * Initialise the instance and calculate for given location, zero altitude and current date
122
+ *
123
+ * @param longitude
124
+ * Longitude in decimal degrees
125
+ * @param latitude
126
+ * Latitude in decimal degrees
127
+ */
128
+ Geomagnetism(double longitude, double latitude) {
129
+ this();
130
+ calculate(longitude, latitude);
131
+ }
132
+
133
+ /**
134
+ * Calculate for given location, altitude and date
135
+ *
136
+ * @param longitude
137
+ * Longitude in decimal degrees
138
+ * @param latitude
139
+ * Latitude in decimal degrees
140
+ * @param altitude
141
+ * Altitude in metres (with respect to WGS-1984 ellipsoid)
142
+ * @param calendar
143
+ * Calendar for date of calculation
144
+ */
145
+ void calculate(double longitude, double latitude, double altitude, GregorianCalendar calendar) {
146
+ double rlon = Math.toRadians(longitude), rlat = Math.toRadians(latitude),
147
+ altitudeKm = Double.isNaN(altitude) ? 0 : altitude / 1000,
148
+ yearFraction = calendar.get(GregorianCalendar.YEAR)
149
+ + (double) calendar.get(GregorianCalendar.DAY_OF_YEAR)
150
+ / calendar.getActualMaximum(GregorianCalendar.DAY_OF_YEAR),
151
+ dt = yearFraction - epoch, srlon = Math.sin(rlon), srlat = Math.sin(rlat), crlon = Math.cos(rlon),
152
+ crlat = Math.cos(rlat), srlat2 = srlat * srlat, crlat2 = crlat * crlat, a2 = WGS84_A * WGS84_A,
153
+ b2 = WGS84_B * WGS84_B, c2 = a2 - b2, a4 = a2 * a2, b4 = b2 * b2, c4 = a4 - b4;
154
+
155
+ sp[1] = srlon;
156
+ cp[1] = crlon;
157
+
158
+ // Convert from geodetic coords. to spherical coords.
159
+ if (altitudeKm != oalt || latitude != olat) {
160
+ double q = Math.sqrt(a2 - c2 * srlat2), q1 = altitudeKm * q,
161
+ q2 = ((q1 + a2) / (q1 + b2)) * ((q1 + a2) / (q1 + b2)),
162
+ r2 = ((altitudeKm * altitudeKm) + 2 * q1 + (a4 - c4 * srlat2) / (q * q));
163
+ ct = srlat / Math.sqrt(q2 * crlat2 + srlat2);
164
+ st = Math.sqrt(1 - (ct * ct));
165
+ r = Math.sqrt(r2);
166
+ d = Math.sqrt(a2 * crlat2 + b2 * srlat2);
167
+ ca = (altitudeKm + d) / r;
168
+ sa = c2 * crlat * srlat / (r * d);
169
+ }
170
+ if (longitude != olon) {
171
+ for (int m = 2; m <= maxord; m++) {
172
+ sp[m] = sp[1] * cp[m - 1] + cp[1] * sp[m - 1];
173
+ cp[m] = cp[1] * cp[m - 1] - sp[1] * sp[m - 1];
174
+ }
175
+ }
176
+ double aor = IAU66_RADIUS / r, ar = aor * aor, br = 0, bt = 0, bp = 0, bpp = 0, par, parp, temp1, temp2;
177
+
178
+ for (int n = 1; n <= maxord; n++) {
179
+ ar = ar * aor;
180
+ for (int m = 0, d3 = 1, d4 = (n + m + d3) / d3; d4 > 0; d4--, m += d3) {
181
+
182
+ // Compute unnormalized associated legendre polynomials and derivatives via recursion relations
183
+ if (altitudeKm != oalt || latitude != olat) {
184
+ if (n == m) {
185
+ snorm[n + m * 13] = st * snorm[n - 1 + (m - 1) * 13];
186
+ dp[m][n] = st * dp[m - 1][n - 1] + ct * snorm[n - 1 + (m - 1) * 13];
187
+ }
188
+ if (n == 1 && m == 0) {
189
+ snorm[n + m * 13] = ct * snorm[n - 1 + m * 13];
190
+ dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13];
191
+ }
192
+ if (n > 1 && n != m) {
193
+ if (m > n - 2)
194
+ snorm[n - 2 + m * 13] = 0;
195
+ if (m > n - 2)
196
+ dp[m][n - 2] = 0;
197
+ snorm[n + m * 13] = ct * snorm[n - 1 + m * 13] - k[m][n] * snorm[n - 2 + m * 13];
198
+ dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13] - k[m][n] * dp[m][n - 2];
199
+ }
200
+ }
201
+
202
+ // Time adjust the gauss coefficients
203
+ if (yearFraction != otime) {
204
+ tc[m][n] = c[m][n] + dt * cd[m][n];
205
+
206
+ if (m != 0)
207
+ tc[n][m - 1] = c[n][m - 1] + dt * cd[n][m - 1];
208
+ }
209
+
210
+ // Accumulate terms of the spherical harmonic expansions
211
+ par = ar * snorm[n + m * 13];
212
+ if (m == 0) {
213
+ temp1 = tc[m][n] * cp[m];
214
+ temp2 = tc[m][n] * sp[m];
215
+ } else {
216
+ temp1 = tc[m][n] * cp[m] + tc[n][m - 1] * sp[m];
217
+ temp2 = tc[m][n] * sp[m] - tc[n][m - 1] * cp[m];
218
+ }
219
+
220
+ bt = bt - ar * temp1 * dp[m][n];
221
+ bp += (fm[m] * temp2 * par);
222
+ br += (fn[n] * temp1 * par);
223
+
224
+ // Special case: north/south geographic poles
225
+ if (st == 0 && m == 1) {
226
+ if (n == 1)
227
+ pp[n] = pp[n - 1];
228
+ else
229
+ pp[n] = ct * pp[n - 1] - k[m][n] * pp[n - 2];
230
+ parp = ar * pp[n];
231
+ bpp += (fm[m] * temp2 * parp);
232
+ }
233
+ }
234
+ }
235
+
236
+ if (st == 0)
237
+ bp = bpp;
238
+ else
239
+ bp /= st;
240
+
241
+ // Rotate magnetic vector components from spherical to geodetic coordinates
242
+ // bx must be the east-west field component
243
+ // by must be the north-south field component
244
+ // bz must be the vertical field component.
245
+ bx = -bt * ca - br * sa;
246
+ by = bp;
247
+ bz = bt * sa - br * ca;
248
+
249
+ // Compute declination (dec), inclination (dip) and total intensity (ti)
250
+ bh = Math.sqrt((bx * bx) + (by * by));
251
+ intensity = Math.sqrt((bh * bh) + (bz * bz));
252
+ // Calculate the declination.
253
+ declination = Math.toDegrees(Math.atan2(by, bx));
254
+ inclination = Math.toDegrees(Math.atan2(bz, bh));
255
+
256
+ otime = yearFraction;
257
+ oalt = altitudeKm;
258
+ olat = latitude;
259
+ olon = longitude;
260
+ }
261
+
262
+ /**
263
+ * Calculate for given location, altitude and current date
264
+ *
265
+ * @param longitude
266
+ * Longitude in decimal degrees
267
+ * @param latitude
268
+ * Latitude in decimal degrees
269
+ * @param altitude
270
+ * Altitude in metres (with respect to WGS-1984 ellipsoid)
271
+ */
272
+ void calculate(double longitude, double latitude, double altitude) {
273
+ calculate(longitude, latitude, altitude, new GregorianCalendar());
274
+ }
275
+
276
+ /**
277
+ * Calculate for given location, zero altitude and current date
278
+ *
279
+ * @param longitude
280
+ * Longitude in decimal degrees
281
+ * @param latitude
282
+ * Latitude in decimal degrees
283
+ */
284
+ void calculate(double longitude, double latitude) {
285
+ calculate(longitude, latitude, 0);
286
+ }
287
+
288
+ /** @return Geomagnetic declination (degrees) [opposite of variation, positive Eastward/negative Westward] */
289
+ double getDeclination() {
290
+ return declination;
291
+ }
292
+
293
+ /** @return Geomagnetic inclination/dip angle (degrees) [positive downward] */
294
+ double getInclination() {
295
+ return inclination;
296
+ }
297
+
298
+ /** @return Geomagnetic field intensity/strength (nano Teslas) */
299
+ double getIntensity() {
300
+ return intensity;
301
+ }
302
+
303
+ /** @return Geomagnetic horizontal field intensity/strength (nano Teslas) */
304
+ double getHorizontalIntensity() {
305
+ return bh;
306
+ }
307
+
308
+ /** @return Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
309
+ double getVerticalIntensity() {
310
+ return bz;
311
+ }
312
+
313
+ /** @return Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
314
+ double getNorthIntensity() {
315
+ return bx;
316
+ }
317
+
318
+ /** @return Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
319
+ double getEastIntensity() {
320
+ return by;
321
+ }
322
+
323
+ /**
324
+ * The input string array which contains each line of input for the wmm.cof input file. The columns in this file are
325
+ * as follows: n, m, gnm, hnm, dgnm, dhnm
326
+ */
327
+ private final static String[] WMM_COF = { " 2020.0 WMM-2020 12/10/2019",
328
+ " 1 0 -29404.5 0.0 6.7 0.0", " 1 1 -1450.7 4652.9 7.7 -25.1",
329
+ " 2 0 -2500.0 0.0 -11.5 0.0", " 2 1 2982.0 -2991.6 -7.1 -30.2",
330
+ " 2 2 1676.8 -734.8 -2.2 -23.9", " 3 0 1363.9 0.0 2.8 0.0",
331
+ " 3 1 -2381.0 -82.2 -6.2 5.7", " 3 2 1236.2 241.8 3.4 -1.0",
332
+ " 3 3 525.7 -542.9 -12.2 1.1", " 4 0 903.1 0.0 -1.1 0.0",
333
+ " 4 1 809.4 282.0 -1.6 0.2", " 4 2 86.2 -158.4 -6.0 6.9",
334
+ " 4 3 -309.4 199.8 5.4 3.7", " 4 4 47.9 -350.1 -5.5 -5.6",
335
+ " 5 0 -234.4 0.0 -0.3 0.0", " 5 1 363.1 47.7 0.6 0.1",
336
+ " 5 2 187.8 208.4 -0.7 2.5", " 5 3 -140.7 -121.3 0.1 -0.9",
337
+ " 5 4 -151.2 32.2 1.2 3.0", " 5 5 13.7 99.1 1.0 0.5",
338
+ " 6 0 65.9 0.0 -0.6 0.0", " 6 1 65.6 -19.1 -0.4 0.1",
339
+ " 6 2 73.0 25.0 0.5 -1.8", " 6 3 -121.5 52.7 1.4 -1.4",
340
+ " 6 4 -36.2 -64.4 -1.4 0.9", " 6 5 13.5 9.0 -0.0 0.1",
341
+ " 6 6 -64.7 68.1 0.8 1.0", " 7 0 80.6 0.0 -0.1 0.0",
342
+ " 7 1 -76.8 -51.4 -0.3 0.5", " 7 2 -8.3 -16.8 -0.1 0.6",
343
+ " 7 3 56.5 2.3 0.7 -0.7", " 7 4 15.8 23.5 0.2 -0.2",
344
+ " 7 5 6.4 -2.2 -0.5 -1.2", " 7 6 -7.2 -27.2 -0.8 0.2",
345
+ " 7 7 9.8 -1.9 1.0 0.3", " 8 0 23.6 0.0 -0.1 0.0",
346
+ " 8 1 9.8 8.4 0.1 -0.3", " 8 2 -17.5 -15.3 -0.1 0.7",
347
+ " 8 3 -0.4 12.8 0.5 -0.2", " 8 4 -21.1 -11.8 -0.1 0.5",
348
+ " 8 5 15.3 14.9 0.4 -0.3", " 8 6 13.7 3.6 0.5 -0.5",
349
+ " 8 7 -16.5 -6.9 0.0 0.4", " 8 8 -0.3 2.8 0.4 0.1",
350
+ " 9 0 5.0 0.0 -0.1 0.0", " 9 1 8.2 -23.3 -0.2 -0.3",
351
+ " 9 2 2.9 11.1 -0.0 0.2", " 9 3 -1.4 9.8 0.4 -0.4",
352
+ " 9 4 -1.1 -5.1 -0.3 0.4", " 9 5 -13.3 -6.2 -0.0 0.1",
353
+ " 9 6 1.1 7.8 0.3 -0.0", " 9 7 8.9 0.4 -0.0 -0.2",
354
+ " 9 8 -9.3 -1.5 -0.0 0.5", " 9 9 -11.9 9.7 -0.4 0.2",
355
+ " 10 0 -1.9 0.0 0.0 0.0", " 10 1 -6.2 3.4 -0.0 -0.0",
356
+ " 10 2 -0.1 -0.2 -0.0 0.1", " 10 3 1.7 3.5 0.2 -0.3",
357
+ " 10 4 -0.9 4.8 -0.1 0.1", " 10 5 0.6 -8.6 -0.2 -0.2",
358
+ " 10 6 -0.9 -0.1 -0.0 0.1", " 10 7 1.9 -4.2 -0.1 -0.0",
359
+ " 10 8 1.4 -3.4 -0.2 -0.1", " 10 9 -2.4 -0.1 -0.1 0.2",
360
+ " 10 10 -3.9 -8.8 -0.0 -0.0", " 11 0 3.0 0.0 -0.0 0.0",
361
+ " 11 1 -1.4 -0.0 -0.1 -0.0", " 11 2 -2.5 2.6 -0.0 0.1",
362
+ " 11 3 2.4 -0.5 0.0 0.0", " 11 4 -0.9 -0.4 -0.0 0.2",
363
+ " 11 5 0.3 0.6 -0.1 -0.0", " 11 6 -0.7 -0.2 0.0 0.0",
364
+ " 11 7 -0.1 -1.7 -0.0 0.1", " 11 8 1.4 -1.6 -0.1 -0.0",
365
+ " 11 9 -0.6 -3.0 -0.1 -0.1", " 11 10 0.2 -2.0 -0.1 0.0",
366
+ " 11 11 3.1 -2.6 -0.1 -0.0", " 12 0 -2.0 0.0 0.0 0.0",
367
+ " 12 1 -0.1 -1.2 -0.0 -0.0", " 12 2 0.5 0.5 -0.0 0.0",
368
+ " 12 3 1.3 1.3 0.0 -0.1", " 12 4 -1.2 -1.8 -0.0 0.1",
369
+ " 12 5 0.7 0.1 -0.0 -0.0", " 12 6 0.3 0.7 0.0 0.0",
370
+ " 12 7 0.5 -0.1 -0.0 -0.0", " 12 8 -0.2 0.6 0.0 0.1",
371
+ " 12 9 -0.5 0.2 -0.0 -0.0", " 12 10 0.1 -0.9 -0.0 -0.0",
372
+ " 12 11 -1.1 -0.0 -0.0 0.0", " 12 12 -0.3 0.5 -0.1 -0.1" };
373
+
374
+ /** Mean radius of IAU-66 ellipsoid, in km */
375
+ private static final double IAU66_RADIUS = 6371.2;
376
+
377
+ /** Semi-major axis of WGS-1984 ellipsoid, in km */
378
+ private static final double WGS84_A = 6378.137;
379
+
380
+ /** Semi-minor axis of WGS-1984 ellipsoid, in km */
381
+ private static final double WGS84_B = 6356.7523142;
382
+
383
+ /** The maximum number of degrees of the spherical harmonic model */
384
+ private static final int MAX_DEG = 12;
385
+
386
+ /** Geomagnetic declination (decimal degrees) [opposite of variation, positive Eastward/negative Westward] */
387
+ private double declination = 0;
388
+
389
+ /** Geomagnetic inclination/dip angle (degrees) [positive downward] */
390
+ private double inclination = 0;
391
+
392
+ /** Geomagnetic field intensity/strength (nano Teslas) */
393
+ private double intensity = 0;
394
+
395
+ /** Geomagnetic horizontal field intensity/strength (nano Teslas) */
396
+ private double bh;
397
+
398
+ /** Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
399
+ private double bz;
400
+
401
+ /** Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
402
+ private double bx;
403
+
404
+ /** Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
405
+ private double by;
406
+
407
+ /** The maximum order of spherical harmonic model */
408
+ private int maxord;
409
+
410
+ /** The Gauss coefficients of main geomagnetic model (nt) */
411
+ private double c[][] = new double[13][13];
412
+
413
+ /** The Gauss coefficients of secular geomagnetic model (nt/yr) */
414
+ private double cd[][] = new double[13][13];
415
+
416
+ /** The time adjusted geomagnetic gauss coefficients (nt) */
417
+ private double tc[][] = new double[13][13];
418
+
419
+ /** The theta derivative of p(n,m) (unnormalized) */
420
+ private double dp[][] = new double[13][13];
421
+
422
+ /** The Schmidt normalization factors */
423
+ private double snorm[] = new double[169];
424
+
425
+ /** The sine of (m*spherical coordinate longitude) */
426
+ private double sp[] = new double[13];
427
+
428
+ /** The cosine of (m*spherical coordinate longitude) */
429
+ private double cp[] = new double[13];
430
+ private double fn[] = new double[13];
431
+ private double fm[] = new double[13];
432
+
433
+ /** The associated Legendre polynomials for m = 1 (unnormalized) */
434
+ private double pp[] = new double[13];
435
+
436
+ private double k[][] = new double[13][13];
437
+
438
+ /**
439
+ * The variables otime (old time), oalt (old altitude), olat (old latitude), olon (old longitude), are used to store
440
+ * the values used from the previous calculation to save on calculation time if some inputs don't change
441
+ */
442
+ private double otime, oalt, olat, olon;
443
+
444
+ /** The date in years, for the start of the valid time of the fit coefficients */
445
+ private double epoch;
446
+
447
+ private double r, d, ca, sa, ct, st;
448
+}
... ...
\ No newline at end of file