bug6178: added WMM coefficients files for 2010-2025

This commit is contained in:
Axel Uhl
2025-11-04 15:53:25 +01:00
parent c165041109
commit 80a48393b7
15 changed files with 574 additions and 161 deletions
@@ -9,6 +9,7 @@ import java.io.InputStreamReader;
import java.text.ParseException;
import java.util.GregorianCalendar;
import org.junit.jupiter.api.BeforeAll;
import org.junit.jupiter.api.Test;
import com.sap.sailing.declination.Declination;
@@ -17,22 +18,55 @@ import com.sap.sailing.domain.common.impl.DegreePosition;
import com.sap.sse.common.TimePoint;
public class WMMFileTest {
private static Geomagnetism g;
@BeforeAll
public static void setUpForClass() throws IOException, ParseException {
g = new Geomagnetism(new BufferedReader(new InputStreamReader(WMMFileTest.class.getClassLoader().getResourceAsStream("WMM2025.COF"))));
}
@Test
public void testWMM2025Reading() throws IOException {
final Geomagnetism g = new Geomagnetism(new BufferedReader(new InputStreamReader(getClass().getClassLoader().getResourceAsStream("WMM2025.COF"))));
assertNotNull(g);
}
@Test
public void testWMM2025Content() throws IOException, ParseException {
final Geomagnetism g = new Geomagnetism(new BufferedReader(new InputStreamReader(getClass().getClassLoader().getResourceAsStream("WMMHR2025.COF"))));
public void testWMM2025Content1() throws IOException, ParseException {
final double lat = 44.123;
final double lng = 8.234;
final GregorianCalendar cal = new GregorianCalendar(2025, 10, 4, 0, 46, 9);
g.calculate(lng, lat, /* altitude */ 0, cal);
final double declinationWMM2025 = g.getDeclination();
assertWMMEqualToOneTenthOfADegreeToDeclinationService(g, lat, lng, cal);
}
@Test
public void testWMM2025Content2() throws IOException, ParseException {
final double lat = 49.234;
final double lng = 9.777;
final GregorianCalendar cal = new GregorianCalendar(2025, 2, 7, 10, 22, 17);
assertWMMEqualToOneTenthOfADegreeToDeclinationService(g, lat, lng, cal);
}
@Test
public void testWMM2025ContentSouthernHemisphere() throws IOException, ParseException {
final double lat = -44.123;
final double lng = -8.234;
final GregorianCalendar cal = new GregorianCalendar(2026, 8, 1, 3, 26, 7);
assertWMMEqualToOneTenthOfADegreeToDeclinationService(g, lat, lng, cal);
}
@Test
public void testWMM2025ContentFarEastEquator() throws IOException, ParseException {
final double lat = 0.000;
final double lng = 170.444;
final GregorianCalendar cal = new GregorianCalendar(2027, 1, 24, 18, 0, 0);
assertWMMEqualToOneTenthOfADegreeToDeclinationService(g, lat, lng, cal);
}
private void assertWMMEqualToOneTenthOfADegreeToDeclinationService(final Geomagnetism g, final double lat,
final double lng, final GregorianCalendar cal) throws IOException, ParseException {
final double declinationWMM2025 = g.calculate(lng, lat, /* altitude */ 0, cal).getDeclination();
final DeclinationService s = DeclinationService.INSTANCE;
final Declination declinationFromService = s.getDeclination(TimePoint.of(cal.getTimeInMillis()), new DegreePosition(lat, lng), /* timeout */ 10000);
assertEquals(declinationWMM2025, declinationFromService.getBearing().getDegrees());
assertEquals(declinationFromService.getBearing().getDegrees(), declinationWMM2025, /* delta in degrees */ 0.15, "Deviation of more than 0.15 degrees; that's too much!");
}
}
@@ -13,7 +13,7 @@ import org.junit.jupiter.api.Test;
import com.sap.sailing.declination.Declination;
import com.sap.sailing.declination.DeclinationService;
import com.sap.sailing.declination.impl.DeclinationImporter;
import com.sap.sailing.declination.impl.DeclinationServiceImpl;
import com.sap.sailing.declination.impl.DeclinationServiceImplWithStore;
import com.sap.sailing.domain.common.impl.CentralAngleDistance;
import com.sap.sailing.domain.common.impl.DegreePosition;
import com.sap.sse.common.impl.MillisecondsTimePoint;
@@ -23,7 +23,7 @@ public abstract class DeclinationServiceTest<I extends DeclinationImporter> exte
@BeforeEach
public void setUp() {
service = new DeclinationServiceImpl(new CentralAngleDistance(1./180.*Math.PI), importer);
service = new DeclinationServiceImplWithStore(new CentralAngleDistance(1./180.*Math.PI), importer);
}
@Test
@@ -0,0 +1,41 @@
package com.sap.sailing.declination.test;
import static org.junit.jupiter.api.Assertions.assertEquals;
import static org.junit.jupiter.api.Assertions.assertNotNull;
import java.io.IOException;
import java.text.ParseException;
import java.text.SimpleDateFormat;
import org.junit.jupiter.api.BeforeAll;
import org.junit.jupiter.api.Test;
import com.sap.sailing.declination.Declination;
import com.sap.sailing.declination.impl.WMMCalculatorDeclinationService;
import com.sap.sailing.domain.common.impl.DegreePosition;
import com.sap.sse.common.impl.MillisecondsTimePoint;
public class WMMDeclinationServiceTest {
private static WMMCalculatorDeclinationService service;
private final SimpleDateFormat simpleDateFormat = new SimpleDateFormat("yyyy-MM-dd");
@BeforeAll
public static void setUp() {
service = new WMMCalculatorDeclinationService();
}
@Test
public void testSimpleDeclinationQueryMatchedInStore() throws IOException, ClassNotFoundException, ParseException {
Declination result = service.getDeclination(new MillisecondsTimePoint(simpleDateFormat.parse("2011-02-03").getTime()),
new DegreePosition(51, -5), /* timeoutForOnlineFetchInMilliseconds */ 3000);
assertEquals(-3.-14./60., result.getBearing().getDegrees(), 0.25);
}
@Test
public void testDeclinationQueryNotMatchedInStore() throws IOException, ClassNotFoundException, ParseException {
Declination result = service.getDeclination(new MillisecondsTimePoint(simpleDateFormat.parse("2020-02-03").getTime()),
new DegreePosition(51, -5), /* timeoutForOnlineFetchInMilliseconds */ 5000);
assertNotNull(result);
assertEquals(-1.531, result.getBearing().getDegrees(), 0.1);
}
}
@@ -0,0 +1,93 @@
2010.0 WMM-2010 12/10/2009
1 0 -29496.6 0.0 11.6 0.0
1 1 -1586.3 4944.4 16.5 -25.9
2 0 -2396.6 0.0 -12.1 0.0
2 1 3026.1 -2707.7 -4.4 -27.2
2 2 1668.6 -576.1 -10.0 -2.2
3 0 1351.1 0.0 -1.4 0.0
3 1 -2352.3 -214.9 -9.4 6.7
3 2 1225.6 245.0 7.2 -3.1
3 3 581.9 -538.3 -2.6 -1.1
4 0 907.3 0.0 0.8 0.0
4 1 813.8 283.9 -0.4 0.6
4 2 120.3 -188.1 -3.1 0.1
4 3 -335.0 180.9 0.8 1.7
4 4 70.3 -329.5 -1.8 2.2
5 0 -232.6 0.0 -0.3 0.0
5 1 360.1 47.0 0.1 0.4
5 2 192.4 196.9 -0.5 -0.6
5 3 -141.1 -119.0 0.1 0.5
5 4 -157.4 16.1 0.6 0.6
5 5 13.7 93.8 0.0 0.0
6 0 67.4 0.0 0.1 0.0
6 1 65.9 -20.6 -0.3 0.1
6 2 73.3 23.6 -0.1 -0.2
6 3 -121.5 54.3 0.0 0.1
6 4 -36.2 -64.4 -0.1 0.0
6 5 13.5 9.0 0.0 0.0
6 6 -64.7 68.0 0.0 0.0
7 0 80.6 0.0 0.0 0.0
7 1 -76.8 -51.4 0.0 0.1
7 2 -8.3 -16.8 0.0 0.0
7 3 56.5 2.3 0.0 0.0
7 4 15.8 23.5 0.0 0.0
7 5 6.4 -2.2 0.0 0.0
7 6 -7.2 8.6 0.0 0.0
7 7 9.8 -3.9 0.0 0.0
8 0 23.6 0.0 0.0 0.0
8 1 9.8 8.4 0.0 0.0
8 2 -17.5 -15.3 0.0 0.0
8 3 -0.4 12.8 0.0 0.0
8 4 -21.1 -7.4 0.0 0.0
8 5 15.3 9.3 0.0 0.0
8 6 13.7 6.7 0.0 0.0
8 7 -16.5 -0.3 0.0 0.0
8 8 -0.3 3.0 0.0 0.0
9 0 5.0 0.0 0.0 0.0
9 1 8.2 -23.3 0.0 0.0
9 2 2.9 11.1 0.0 0.0
9 3 -1.4 9.8 0.0 0.0
9 4 -1.1 -6.4 0.0 0.0
9 5 -13.3 1.5 0.0 0.0
9 6 1.1 4.2 0.0 0.0
9 7 8.9 0.1 0.0 0.0
9 8 -9.3 -1.4 0.0 0.0
9 9 -11.9 9.9 0.0 0.0
10 0 -1.9 0.0 0.0 0.0
10 1 -6.2 3.4 0.0 0.0
10 2 0.3 -0.4 0.0 0.0
10 3 0.6 -0.8 0.0 0.0
10 4 1.7 -0.7 0.0 0.0
10 5 -0.5 0.3 0.0 0.0
10 6 0.2 -1.7 0.0 0.0
10 7 1.7 1.3 0.0 0.0
10 8 -0.2 0.2 0.0 0.0
10 9 0.3 0.4 0.0 0.0
10 10 -1.2 0.0 0.0 0.0
11 0 3.1 0.0 0.0 0.0
11 1 -1.5 -0.9 0.0 0.0
11 2 -2.3 -2.2 0.0 0.0
11 3 2.1 -0.3 0.0 0.0
11 4 -0.9 0.3 0.0 0.0
11 5 0.6 0.9 0.0 0.0
11 6 0.5 0.2 0.0 0.0
11 7 -0.4 0.1 0.0 0.0
11 8 0.1 0.5 0.0 0.0
11 9 -0.1 0.4 0.0 0.0
11 10 -0.3 -0.2 0.0 0.0
11 11 0.2 -0.2 0.0 0.0
12 0 -2.0 0.0 0.0 0.0
12 1 -0.1 0.6 0.0 0.0
12 2 0.5 -0.6 0.0 0.0
12 3 1.3 0.1 0.0 0.0
12 4 -0.5 0.3 0.0 0.0
12 5 0.4 -0.1 0.0 0.0
12 6 0.1 0.1 0.0 0.0
12 7 0.0 -0.1 0.0 0.0
12 8 0.2 0.2 0.0 0.0
12 9 -0.3 0.0 0.0 0.0
12 10 -0.1 0.0 0.0 0.0
12 11 0.0 0.0 0.0 0.0
12 12 0.1 0.0 0.0 0.0
999999999999999999999999999999999999999999999999
999999999999999999999999999999999999999999999999
@@ -0,0 +1,93 @@
2015.0 WMM-2015 12/10/2014
1 0 -29404.8 0.0 6.7 0.0
1 1 -1450.9 4652.5 7.7 -25.9
2 0 -2500.0 0.0 -11.5 0.0
2 1 2982.0 -2991.6 -7.1 -30.2
2 2 1677.0 -734.8 -2.2 -23.9
3 0 1363.9 0.0 2.8 0.0
3 1 -2381.0 -82.2 -6.2 5.2
3 2 1236.2 241.8 3.4 -1.0
3 3 525.7 -542.9 -12.2 1.1
4 0 903.1 0.0 -1.1 0.0
4 1 809.4 282.0 -1.6 0.2
4 2 86.2 -158.4 -6.0 6.9
4 3 -309.4 199.8 5.4 3.7
4 4 47.9 -350.1 -5.5 8.6
5 0 -234.4 0.0 -0.5 0.0
5 1 363.1 47.7 0.5 1.0
5 2 187.8 208.4 -0.7 -0.8
5 3 -140.7 -121.3 0.8 1.0
5 4 -151.2 32.3 1.2 3.0
5 5 13.7 99.1 1.0 0.5
6 0 65.9 0.0 -0.6 0.0
6 1 65.6 -19.1 -0.4 0.1
6 2 73.0 25.0 -0.2 -0.4
6 3 -121.5 52.7 0.4 1.2
6 4 -36.2 -64.4 -0.7 0.6
6 5 13.5 9.0 0.1 0.1
6 6 -64.7 68.0 0.8 -0.1
7 0 80.6 0.0 -0.1 0.0
7 1 -76.8 -51.4 -0.3 0.5
7 2 -8.3 -16.8 0.1 0.1
7 3 56.5 2.3 0.4 -0.7
7 4 15.8 23.5 0.0 -0.5
7 5 6.4 -2.2 -0.3 0.3
7 6 -7.2 8.6 0.0 0.0
7 7 9.8 -3.9 0.3 0.2
8 0 23.6 0.0 -0.3 0.0
8 1 9.8 8.4 -0.1 0.1
8 2 -17.5 -15.3 0.3 0.1
8 3 -0.4 12.8 0.1 -0.3
8 4 -21.1 -7.4 0.0 0.0
8 5 15.3 9.3 -0.1 0.2
8 6 13.7 6.7 0.1 0.0
8 7 -16.5 -0.3 0.2 0.0
8 8 -0.3 3.0 0.0 0.0
9 0 5.0 0.0 0.0 0.0
9 1 8.2 -23.3 -0.1 -0.2
9 2 2.9 11.1 -0.1 -0.1
9 3 -1.4 9.8 0.1 0.2
9 4 -1.1 -6.4 0.0 0.1
9 5 -13.3 1.5 0.3 0.0
9 6 1.1 4.2 0.0 0.0
9 7 8.9 0.1 -0.1 0.0
9 8 -9.3 -1.4 -0.1 0.0
9 9 -11.9 9.9 0.1 0.0
10 0 -1.9 0.0 0.0 0.0
10 1 -6.2 3.4 0.0 0.0
10 2 0.3 -0.4 0.0 0.0
10 3 0.6 -0.8 0.0 0.0
10 4 1.7 -0.7 0.0 0.0
10 5 -0.5 0.3 0.0 0.0
10 6 0.2 -1.7 0.0 0.0
10 7 1.7 1.3 0.0 0.0
10 8 -0.2 0.2 0.0 0.0
10 9 0.3 0.4 0.0 0.0
10 10 -1.2 0.0 0.0 0.0
11 0 3.1 0.0 0.0 0.0
11 1 -1.5 -0.9 0.0 0.0
11 2 -2.3 -2.2 0.0 0.0
11 3 2.1 -0.3 0.0 0.0
11 4 -0.9 0.3 0.0 0.0
11 5 0.6 0.9 0.0 0.0
11 6 0.5 0.2 0.0 0.0
11 7 -0.4 0.1 0.0 0.0
11 8 0.1 0.5 0.0 0.0
11 9 -0.1 0.4 0.0 0.0
11 10 -0.3 -0.2 0.0 0.0
11 11 0.2 -0.2 0.0 0.0
12 0 -2.0 0.0 0.0 0.0
12 1 -0.1 0.6 0.0 0.0
12 2 0.5 -0.6 0.0 0.0
12 3 1.3 0.1 0.0 0.0
12 4 -0.5 0.3 0.0 0.0
12 5 0.4 -0.1 0.0 0.0
12 6 0.1 0.1 0.0 0.0
12 7 0.0 -0.1 0.0 0.0
12 8 0.2 0.2 0.0 0.0
12 9 -0.3 0.0 0.0 0.0
12 10 -0.1 0.0 0.0 0.0
12 11 0.0 0.0 0.0 0.0
12 12 0.1 0.0 0.0 0.0
999999999999999999999999999999999999999999999999
999999999999999999999999999999999999999999999999
@@ -3,12 +3,10 @@ package com.sap.sailing.declination;
import java.io.IOException;
import java.text.ParseException;
import com.sap.sailing.declination.impl.BGSImporter;
import com.sap.sailing.declination.impl.DeclinationImporter;
import com.sap.sailing.declination.impl.DeclinationServiceImpl;
import com.sap.sailing.declination.impl.WMMCalculatorDeclinationService;
import com.sap.sailing.domain.common.Mile;
import com.sap.sailing.domain.common.Position;
import com.sap.sailing.domain.common.impl.CentralAngleDistance;
import com.sap.sse.common.Distance;
import com.sap.sse.common.TimePoint;
@@ -17,7 +15,7 @@ public interface DeclinationService {
* A default implementation with a spatial default precision of 60 {@link Mile#METERS_PER_GEOGRAPHICAL_MILE
* nautical miles} which equals the length of an arc with one degree on a meridian.
*/
DeclinationService INSTANCE = new DeclinationServiceImpl(new CentralAngleDistance(1./180.*Math.PI), new BGSImporter());
DeclinationService INSTANCE = new WMMCalculatorDeclinationService();
/**
* Obtains declination information with the default precision of this declination service in time and space
@@ -0,0 +1,39 @@
package com.sap.sailing.declination.impl;
import com.sap.sailing.declination.Declination;
import com.sap.sailing.domain.common.Position;
import com.sap.sse.common.Bearing;
import com.sap.sse.common.TimePoint;
public abstract class AbstractDeclinationRecord implements Declination {
private static final long serialVersionUID = 2701783831406402231L;
private final Position position;
private final TimePoint timePoint;
private final Bearing bearing;
public AbstractDeclinationRecord(Position position, TimePoint timePoint, Bearing bearing) {
this.position = position;
this.timePoint = timePoint;
this.bearing = bearing;
}
@Override
public Position getPosition() {
return position;
}
@Override
public TimePoint getTimePoint() {
return timePoint;
}
@Override
public Bearing getBearing() {
return bearing;
}
@Override
public String toString() {
return ""+getTimePoint()+"@"+getPosition()+": "+getBearing();
}
}
@@ -67,7 +67,8 @@ import com.sap.sse.common.impl.MillisecondsTimePoint;
public class BGSImporter extends DeclinationImporter {
// private static final String URL_PATTERN = "http://geomag.bgs.ac.uk/web_service/GMModels/bggm/2015/?latitude=%f&longitude=%f&altitude=0&date=%d-%d-%d&format=xml";
// private static final String URL_PATTERN = "http://geomag.bgs.ac.uk/web_service/GMModels/wmm/2020/?latitude=%f&longitude=%f&altitude=0&date=%d-%d-%d&format=xml";
private static final String URL_PATTERN = "http://geomag.bgs.ac.uk/web_service/GMModels/igrf/13/?latitude=%f&longitude=%f&altitude=0&date=%d-%d-%d&format=xml";
private static final String URL_PATTERN_PRE_2025 = "http://geomag.bgs.ac.uk/web_service/GMModels/igrf/13/?latitude=%f&longitude=%f&altitude=0&date=%d-%d-%d&format=xml";
private static final String URL_PATTERN_POST_2025 = "https://geomag.bgs.ac.uk/web_service/GMModels/wmm/2025?latitude=%f&longitude=%f&altitude=0&date=%d-%d-%d&format=xml";
private static class XmlElementHandler extends DefaultHandler {
private Date dateAsDecimalYear;
@@ -156,7 +157,8 @@ public class BGSImporter extends DeclinationImporter {
throws IOException, ParserConfigurationException, SAXException {
final Calendar cal = new GregorianCalendar(TimeZone.getTimeZone("UTC"));
cal.setTime(timePoint.asDate());
final URL url = new URL(String.format(URL_PATTERN, position.getLatDeg(), position.getLngDeg(), cal.get(Calendar.YEAR), cal.get(Calendar.MONTH)+1, cal.get(Calendar.DAY_OF_MONTH)));
final int year = cal.get(Calendar.YEAR);
final URL url = new URL(String.format(year >= 2025 ? URL_PATTERN_POST_2025 : URL_PATTERN_PRE_2025, position.getLatDeg(), position.getLngDeg(), year, cal.get(Calendar.MONTH)+1, cal.get(Calendar.DAY_OF_MONTH)));
return getDeclinationFromXml(url.openStream());
}
}
@@ -0,0 +1,41 @@
package com.sap.sailing.declination.impl;
import com.sap.sailing.declination.Declination;
import com.sap.sailing.domain.common.Position;
import com.sap.sse.common.Bearing;
import com.sap.sse.common.TimePoint;
import com.sap.sse.common.Util;
/**
* A declination record that does not know about any annual change because it has been computed precisely for
* the time point and position requested. Therefore, no {@link #getBearingCorrectedTo(TimePoint) correction} to
* any time point other than the one {@link #getTimePoint() requested} can be performed. An exception will
* be thrown if that happens.
*
* @author Axel Uhl (d043530)
*
*/
public class DeclinationRecordForExactTimePoint extends AbstractDeclinationRecord implements Declination {
private static final long serialVersionUID = -94512743120385233L;
public DeclinationRecordForExactTimePoint(Position position, TimePoint timePoint, Bearing bearing) {
super(position, timePoint, bearing);
}
@Override
public Bearing getAnnualChange() {
return null;
}
@Override
public Bearing getBearingCorrectedTo(TimePoint timePoint) {
if (Util.equalsWithNull(timePoint, getTimePoint())) {
return getBearing();
} else {
throw new IllegalArgumentException("Declination computed precisely for " + getTimePoint()
+ " cannot be corrected to any other time point " + timePoint
+ " because we lack knowledge of the annual change here.");
}
}
}
@@ -6,35 +6,15 @@ import com.sap.sse.common.Bearing;
import com.sap.sse.common.TimePoint;
import com.sap.sse.common.impl.DegreeBearingImpl;
public class DeclinationRecordImpl implements Declination {
public class DeclinationRecordImpl extends AbstractDeclinationRecord implements Declination {
private static final long serialVersionUID = 6918630656182340186L;
private final Position position;
private final TimePoint timePoint;
private final Bearing bearing;
private final Bearing annualChange;
public DeclinationRecordImpl(Position position, TimePoint timePoint, Bearing bearing, Bearing annualChange) {
super();
this.position = position;
this.timePoint = timePoint;
this.bearing = bearing;
super(position, timePoint, bearing);
this.annualChange = annualChange;
}
@Override
public Position getPosition() {
return position;
}
@Override
public TimePoint getTimePoint() {
return timePoint;
}
@Override
public Bearing getBearing() {
return bearing;
}
@Override
public Bearing getAnnualChange() {
return annualChange;
@@ -50,5 +30,4 @@ public class DeclinationRecordImpl implements Declination {
public String toString() {
return ""+getTimePoint()+"@"+getPosition()+": "+getBearing()+", "+getAnnualChange()+"/year";
}
}
@@ -14,7 +14,7 @@ import com.sap.sailing.domain.common.quadtree.QuadTree;
import com.sap.sse.common.Distance;
import com.sap.sse.common.TimePoint;
public class DeclinationServiceImpl implements DeclinationService {
public class DeclinationServiceImplWithStore implements DeclinationService {
private final Distance defaultMaxDistance;
private final DeclinationStore persistentStore;
private final DeclinationImporter declinationImporter;
@@ -29,7 +29,7 @@ public class DeclinationServiceImpl implements DeclinationService {
*/
private final Map<Integer, QuadTree<Declination>> importerCache;
public DeclinationServiceImpl(Distance defaultMaxDistance, DeclinationImporter declinationImporter) {
public DeclinationServiceImplWithStore(Distance defaultMaxDistance, DeclinationImporter declinationImporter) {
this.declinationImporter = declinationImporter;
this.defaultMaxDistance = defaultMaxDistance;
this.yearStore = new HashMap<Integer, QuadTree<Declination>>();
@@ -207,7 +207,7 @@ public class DeclinationStore {
existingDeclinationRecord = storedDeclinations.get(point);
}
if (existingDeclinationRecord == null
|| DeclinationServiceImpl.timeAndSpaceDistance(existingDeclinationRecord.getPosition().getDistance(point),
|| DeclinationServiceImplWithStore.timeAndSpaceDistance(existingDeclinationRecord.getPosition().getDistance(point),
timePoint, existingDeclinationRecord.getTimePoint()) > 0.1) {
// less than ~6 nautical miles and/or ~.6 months off
fetchAndAppendDeclination(timePoint, point, importer, out);
@@ -2,6 +2,8 @@ package com.sap.sailing.declination.impl;
import java.io.BufferedReader;
import java.io.IOException;
import java.text.ParseException;
import java.text.SimpleDateFormat;
/* License Statement from the NOAA
* The WMM source code is in the public domain and not licensed or
@@ -14,24 +16,115 @@ import java.io.IOException;
import java.util.GregorianCalendar;
import com.sap.sse.common.Duration;
import com.sap.sse.common.Named;
import com.sap.sse.common.TimePoint;
/**
* <p>
* Class to calculate magnetic declination, magnetic field strength, inclination etc. for any point on the earth.
* </p>
* <p>
* Adapted from the geomagc software and World Magnetic Model of the NOAA Satellite and Information Service, National
* Geophysical Data Center
* Geophysical Data Center. Caching removed to make stateless, except for the reading of the coefficients file.
* Results are now returned as instances of the inner class {@link Result}.
* </p>
* http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
* <p>
* © Deep Pradhan, 2017
* </p>
*/
class Geomagnetism {
class Geomagnetism implements Named {
private static final long serialVersionUID = -2814152697634383730L;
class Result {
/** Geomagnetic declination (decimal degrees) [opposite of variation, positive Eastward/negative Westward] */
private final double declination;
/** Geomagnetic inclination/dip angle (degrees) [positive downward] */
private final double inclination;
/** Geomagnetic field intensity/strength (nano Teslas) */
private final double intensity;
/** Geomagnetic horizontal field intensity/strength (nano Teslas) */
private final double bh;
/** Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
private final double bz;
/** Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
private final double bx;
/** Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
private final double by;
public Result(double declination, double inclination, double intensity, double bh, double bz, double bx, double by) {
this.declination = declination;
this.inclination = inclination;
this.intensity = intensity;
this.bh = bh;
this.bz = bz;
this.bx = bx;
this.by = by;
}
/** @return Geomagnetic declination (degrees) [opposite of variation, positive Eastward/negative Westward] */
double getDeclination() {
return declination;
}
/** @return Geomagnetic inclination/dip angle (degrees) [positive downward] */
double getInclination() {
return inclination;
}
/** @return Geomagnetic field intensity/strength (nano Teslas) */
double getIntensity() {
return intensity;
}
/** @return Geomagnetic horizontal field intensity/strength (nano Teslas) */
double getHorizontalIntensity() {
return bh;
}
/** @return Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
double getVerticalIntensity() {
return bz;
}
/** @return Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
double getNorthIntensity() {
return bx;
}
/** @return Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
double getEastIntensity() {
return by;
}
String getModelName() {
return getName();
}
TimePoint getModelIssueTimePoint() {
return getIssueTimePoint();
}
}
private final static SimpleDateFormat simpleDateFormat = new SimpleDateFormat("MM/dd/yyyy");
private final TimePoint startOfValidity;
/**
* Initialise the instance without calculations
*
* @param cofReader
* a reader for a file in .COF format, such as WMM2025.COF; expects to find trailing lines with at least
* "999999" in them
*/
Geomagnetism(BufferedReader r) throws IOException {
Geomagnetism(BufferedReader cofReader) throws IOException, ParseException {
// Initialize constants
maxord = MAX_DEG;
sp[0] = 0;
@@ -39,12 +132,16 @@ class Geomagnetism {
dp[0][0] = 0;
c[0][0] = 0;
cd[0][0] = 0;
final String headerLine = r.readLine();
epoch = Double.parseDouble(headerLine.trim().split("\\s+")[0]);
final String headerLine = cofReader.readLine();
final String[] headerFields = headerLine.trim().split("\\s+");
epoch = Double.parseDouble(headerFields[0]);
startOfValidity = TimePoint.of(new GregorianCalendar((int) epoch, 1, 1, 0, 0, 0).getTimeInMillis()).plus(Duration.ONE_HOUR.times(365.0*24.0*(epoch-((int) epoch))));
name = headerFields[1];
issueTimePoint = TimePoint.of(simpleDateFormat.parse(headerFields[2]));
String[] tokens;
double gnm, hnm, dgnm, dhnm;
String line;
while (!(line=r.readLine()).startsWith("999999")) {
while (!(line=cofReader.readLine()).startsWith("999999")) {
tokens = line.trim().split("\\s+");
final int n = Integer.parseInt(tokens[0]);
final int m = Integer.parseInt(tokens[1]);
@@ -84,7 +181,6 @@ class Geomagnetism {
}
k[1][1] = 0;
fm[0] = 0;
otime = oalt = olat = olon = -1000;
}
/**
@@ -99,7 +195,7 @@ class Geomagnetism {
* @param calendar
* Calendar for date of calculation
*/
void calculate(double longitude, double latitude, double altitude, GregorianCalendar calendar) {
Result calculate(double longitude, double latitude, double altitude, GregorianCalendar calendar) {
double rlon = Math.toRadians(longitude), rlat = Math.toRadians(latitude),
altitudeKm = Double.isNaN(altitude) ? 0 : altitude / 1000,
yearFraction = calendar.get(GregorianCalendar.YEAR)
@@ -111,54 +207,46 @@ class Geomagnetism {
sp[1] = srlon;
cp[1] = crlon;
// Convert from geodetic coords. to spherical coords.
if (altitudeKm != oalt || latitude != olat) {
double q = Math.sqrt(a2 - c2 * srlat2), q1 = altitudeKm * q,
q2 = ((q1 + a2) / (q1 + b2)) * ((q1 + a2) / (q1 + b2)),
r2 = ((altitudeKm * altitudeKm) + 2 * q1 + (a4 - c4 * srlat2) / (q * q));
ct = srlat / Math.sqrt(q2 * crlat2 + srlat2);
st = Math.sqrt(1 - (ct * ct));
r = Math.sqrt(r2);
d = Math.sqrt(a2 * crlat2 + b2 * srlat2);
ca = (altitudeKm + d) / r;
sa = c2 * crlat * srlat / (r * d);
}
if (longitude != olon) {
for (int m = 2; m <= maxord; m++) {
sp[m] = sp[1] * cp[m - 1] + cp[1] * sp[m - 1];
cp[m] = cp[1] * cp[m - 1] - sp[1] * sp[m - 1];
}
double q = Math.sqrt(a2 - c2 * srlat2), q1 = altitudeKm * q,
q2 = ((q1 + a2) / (q1 + b2)) * ((q1 + a2) / (q1 + b2)),
r2 = ((altitudeKm * altitudeKm) + 2 * q1 + (a4 - c4 * srlat2) / (q * q));
ct = srlat / Math.sqrt(q2 * crlat2 + srlat2);
st = Math.sqrt(1 - (ct * ct));
r = Math.sqrt(r2);
d = Math.sqrt(a2 * crlat2 + b2 * srlat2);
ca = (altitudeKm + d) / r;
sa = c2 * crlat * srlat / (r * d);
for (int m = 2; m <= maxord; m++) {
sp[m] = sp[1] * cp[m - 1] + cp[1] * sp[m - 1];
cp[m] = cp[1] * cp[m - 1] - sp[1] * sp[m - 1];
}
double aor = IAU66_RADIUS / r, ar = aor * aor, br = 0, bt = 0, bp = 0, bpp = 0, par, parp, temp1, temp2;
for (int n = 1; n <= maxord; n++) {
ar = ar * aor;
for (int m = 0, d3 = 1, d4 = (n + m + d3) / d3; d4 > 0; d4--, m += d3) {
// Compute unnormalized associated legendre polynomials and derivatives via recursion relations
if (altitudeKm != oalt || latitude != olat) {
if (n == m) {
snorm[n + m * 13] = st * snorm[n - 1 + (m - 1) * 13];
dp[m][n] = st * dp[m - 1][n - 1] + ct * snorm[n - 1 + (m - 1) * 13];
if (n == m) {
snorm[n + m * 13] = st * snorm[n - 1 + (m - 1) * 13];
dp[m][n] = st * dp[m - 1][n - 1] + ct * snorm[n - 1 + (m - 1) * 13];
}
if (n == 1 && m == 0) {
snorm[n + m * 13] = ct * snorm[n - 1 + m * 13];
dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13];
}
if (n > 1 && n != m) {
if (m > n - 2) {
snorm[n - 2 + m * 13] = 0;
}
if (n == 1 && m == 0) {
snorm[n + m * 13] = ct * snorm[n - 1 + m * 13];
dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13];
}
if (n > 1 && n != m) {
if (m > n - 2) {
snorm[n - 2 + m * 13] = 0;
}
if (m > n - 2) {
dp[m][n - 2] = 0;
}
snorm[n + m * 13] = ct * snorm[n - 1 + m * 13] - k[m][n] * snorm[n - 2 + m * 13];
dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13] - k[m][n] * dp[m][n - 2];
if (m > n - 2) {
dp[m][n - 2] = 0;
}
snorm[n + m * 13] = ct * snorm[n - 1 + m * 13] - k[m][n] * snorm[n - 2 + m * 13];
dp[m][n] = ct * dp[m][n - 1] - st * snorm[n - 1 + m * 13] - k[m][n] * dp[m][n - 2];
}
// Time adjust the gauss coefficients
if (yearFraction != otime) {
tc[m][n] = c[m][n] + dt * cd[m][n];
if (m != 0) {
tc[n][m - 1] = c[n][m - 1] + dt * cd[n][m - 1];
}
tc[m][n] = c[m][n] + dt * cd[m][n];
if (m != 0) {
tc[n][m - 1] = c[n][m - 1] + dt * cd[n][m - 1];
}
// Accumulate terms of the spherical harmonic expansions
par = ar * snorm[n + m * 13];
@@ -193,19 +281,16 @@ class Geomagnetism {
// bx must be the east-west field component
// by must be the north-south field component
// bz must be the vertical field component.
bx = -bt * ca - br * sa;
by = bp;
bz = bt * sa - br * ca;
final double bx = -bt * ca - br * sa;
final double by = bp;
final double bz = bt * sa - br * ca;
// Compute declination (dec), inclination (dip) and total intensity (ti)
bh = Math.sqrt((bx * bx) + (by * by));
intensity = Math.sqrt((bh * bh) + (bz * bz));
final double bh = Math.sqrt((bx * bx) + (by * by));
final double intensity = Math.sqrt((bh * bh) + (bz * bz));
// Calculate the declination.
declination = Math.toDegrees(Math.atan2(by, bx));
inclination = Math.toDegrees(Math.atan2(bz, bh));
otime = yearFraction;
oalt = altitudeKm;
olat = latitude;
olon = longitude;
final double declination = Math.toDegrees(Math.atan2(by, bx));
final double inclination = Math.toDegrees(Math.atan2(bz, bh));
return new Result(declination, inclination, intensity, bh, bz, bx, by);
}
/**
@@ -218,8 +303,8 @@ class Geomagnetism {
* @param altitude
* Altitude in metres (with respect to WGS-1984 ellipsoid)
*/
void calculate(double longitude, double latitude, double altitude) {
calculate(longitude, latitude, altitude, new GregorianCalendar());
Result calculate(double longitude, double latitude, double altitude) {
return calculate(longitude, latitude, altitude, new GregorianCalendar());
}
/**
@@ -230,45 +315,30 @@ class Geomagnetism {
* @param latitude
* Latitude in decimal degrees
*/
void calculate(double longitude, double latitude) {
calculate(longitude, latitude, 0);
Result calculate(double longitude, double latitude) {
return calculate(longitude, latitude, 0);
}
/** @return The date in years, for the start of the valid time of the fit coefficients */
public double getEpoch() {
return epoch;
}
/** @return Geomagnetic declination (degrees) [opposite of variation, positive Eastward/negative Westward] */
double getDeclination() {
return declination;
public TimePoint getStartOfValidity() {
return startOfValidity;
}
/** @return Geomagnetic inclination/dip angle (degrees) [positive downward] */
double getInclination() {
return inclination;
public String getName() {
return name;
}
/** @return Geomagnetic field intensity/strength (nano Teslas) */
double getIntensity() {
return intensity;
}
/** @return Geomagnetic horizontal field intensity/strength (nano Teslas) */
double getHorizontalIntensity() {
return bh;
}
/** @return Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
double getVerticalIntensity() {
return bz;
}
/** @return Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
double getNorthIntensity() {
return bx;
}
/** @return Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
double getEastIntensity() {
return by;
public TimePoint getIssueTimePoint() {
return issueTimePoint;
}
private final String name;
private final TimePoint issueTimePoint;
/** Mean radius of IAU-66 ellipsoid, in km */
private static final double IAU66_RADIUS = 6371.2;
@@ -281,30 +351,6 @@ class Geomagnetism {
/** The maximum number of degrees of the spherical harmonic model */
private static final int MAX_DEG = 12;
/** Geomagnetic declination (decimal degrees) [opposite of variation, positive Eastward/negative Westward] */
private double declination = 0;
/** Geomagnetic inclination/dip angle (degrees) [positive downward] */
private double inclination = 0;
/** Geomagnetic field intensity/strength (nano Teslas) */
private double intensity = 0;
/** Geomagnetic horizontal field intensity/strength (nano Teslas) */
private double bh;
/** Geomagnetic vertical field intensity/strength (nano Teslas) [positive downward] */
private double bz;
/** Geomagnetic North South (northerly component) field intensity/strength (nano Tesla) */
private double bx;
/** Geomagnetic East West (easterly component) field intensity/strength (nano Teslas) */
private double by;
/** The maximum order of spherical harmonic model */
private int maxord;
/** The Gauss coefficients of main geomagnetic model (nt) */
private double c[][] = new double[300][300];
@@ -328,17 +374,14 @@ class Geomagnetism {
private double fn[] = new double[13];
private double fm[] = new double[13];
/** The maximum order of spherical harmonic model */
private int maxord;
/** The associated Legendre polynomials for m = 1 (unnormalized) */
private double pp[] = new double[13];
private double k[][] = new double[13][13];
/**
* The variables otime (old time), oalt (old altitude), olat (old latitude), olon (old longitude), are used to store
* the values used from the previous calculation to save on calculation time if some inputs don't change
*/
private double otime, oalt, olat, olon;
/** The date in years, for the start of the valid time of the fit coefficients */
private double epoch;
@@ -0,0 +1,50 @@
package com.sap.sailing.declination.impl;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.text.ParseException;
import java.util.GregorianCalendar;
import java.util.TreeMap;
import com.sap.sailing.declination.Declination;
import com.sap.sailing.declination.DeclinationService;
import com.sap.sailing.declination.impl.Geomagnetism.Result;
import com.sap.sailing.domain.common.Position;
import com.sap.sse.common.Distance;
import com.sap.sse.common.TimePoint;
import com.sap.sse.common.impl.DegreeBearingImpl;
public class WMMCalculatorDeclinationService implements DeclinationService {
/**
* The magnetic models, in ascending order by their {@link Geomagnetism#getIssueTimePoint() issue time point}
*/
private final TreeMap<TimePoint, Geomagnetism> worldMagneticModelsByIssueTimePoint;
public WMMCalculatorDeclinationService() {
final String[] modelNames = new String[] { "/WMM2010.COF", "/WMM2015.COF", "/WMM2020.COF", "/WMM2025.COF" };
worldMagneticModelsByIssueTimePoint = new TreeMap<>();
try {
for (final String modelFileName : modelNames) {
final Geomagnetism model = new Geomagnetism(new BufferedReader(new InputStreamReader(getClass().getResourceAsStream(modelFileName))));
worldMagneticModelsByIssueTimePoint.put(model.getStartOfValidity(), model);
}
} catch (IOException | ParseException e) {
throw new RuntimeException(e);
}
}
@Override
public Declination getDeclination(TimePoint timePoint, Position position, long timeoutForOnlineFetchInMilliseconds) throws IOException, ParseException {
final GregorianCalendar calendar = new GregorianCalendar();
calendar.setTimeInMillis(timePoint.asMillis());
final Result wmmResult = worldMagneticModelsByIssueTimePoint.floorEntry(timePoint).getValue().calculate(position.getLngDeg(), position.getLatDeg(), /* altitude */ 0, calendar);
return new DeclinationRecordForExactTimePoint(position, timePoint, new DegreeBearingImpl(wmmResult.getDeclination()));
}
@Override
public Declination getDeclination(TimePoint timePoint, Position position, Distance maxDistance,
long timeoutForOnlineFetchInMilliseconds) throws IOException, ParseException {
return getDeclination(timePoint, position, timeoutForOnlineFetchInMilliseconds);
}
}