Add magField2Elements function.

README.md

@@ -66,7 +66,10 @@ void loop() {
 If you have a position in latitude, longitude, and height, 
 you can convert it to geocentric cartesian coordinates 
 with `geodetic2ecef`. Note that `geodetic2ecef` uses 
-single precision floats, so it will only be accurate to about 1 meter:
+single precision floats, so it will only be accurate to about 1 meter.
+You can also convert the magnetic field to 
+the [seven magnetic elements](https://www.ncei.noaa.gov/products/world-magnetic-model) 
+in units of nanoTesla and degrees.
 ~~~cpp
 #include "geomag.hpp"
 void setup() {
@@ -81,14 +84,24 @@ void loop() {
   float lat = val + 43.0f; // latitude in degrees
   float lon = val + 75.0f; // longitude in degrees
   float height = val + 305; // height above WGS84 ellipsoid in meters
-  geomag::Vector in = geomag::geodetic2ecef(lat,lon,height);
-  geomag::Vector out = geomag::GeoMag(2022.5,in,geomag::WMM2020);
-  Serial.print(out.x*1E9);
-  Serial.println(" nT x");
-  Serial.print(out.y*1E9);
-  Serial.println(" nT y");
-  Serial.print(out.z*1E9);
-  Serial.println(" nT z");
+  geomag::Vector position = geomag::geodetic2ecef(lat,lon,height);
+  geomag::Vector mag_field = geomag::GeoMag(2022.5,position,geomag::WMM2020);
+  geomag::Elements out = geomag::magField2Elements(mag_field, lat, lon);
+  Serial.print(out.north);
+  Serial.println(" nT north");
+  Serial.print(out.east);
+  Serial.println(" nT east");
+  Serial.print(out.down);
+  Serial.println(" nT down");
+  Serial.print(out.horizontal);
+  Serial.println(" nT horizontal");
+  Serial.print(out.total);
+  Serial.println(" nT total");
+  Serial.print(out.inclination);
+  Serial.println(" deg inclination");
+  Serial.print(out.declination);
+  Serial.println(" deg declination");
+  Serial.println();
   delay(2000);
 }
 ~~~

geomag.hpp

@@ -81,6 +81,50 @@ typedef struct {
 } Vector;
 
 
+typedef struct {
+    float north;// local north magnetic field (nT)
+    float east;// local east magnetic field (nT)
+    float down;// local down magnetic field (nT)
+    float horizontal;// local horizontal magnetic field intensity (nT)
+    float total;// local total magnetic field intensity (nT)
+    float inclination;// also called the dip angle, 
+    // the angle measured from the horizontal plane to the 
+    // magnetic field vector; a downward field is positive (deg)
+    float declination;// also called the magnetic variation,
+    // the angle between true north and the horizontal component 
+    // of the field, a eastward magnetic field of true North is positive (deg)
+} Elements;
+
+/** Return a struct containing the 7 magnetic elements.
+See https://www.geomag.nrcan.gc.ca/mag_fld/comp-en.php and
+https://www.ngdc.noaa.gov/geomag/icons/faqelems.gif for more info.
+ INPUT:
+    mag_field_itrs: local magnetic field in the itrs coordinate system (T)
+    lat: latitude in degrees, -90 at the south pole, 90 at the north pole.
+    lon: longitude in degrees.
+**/
+inline Elements magField2Elements(Vector mag_field_itrs, float lat, float lon){
+    float x = mag_field_itrs.x*1E9f;
+    float y = mag_field_itrs.y*1E9f;
+    float z = mag_field_itrs.z*1E9f;
+    float phi = lat*((float)(M_PI/180.0));
+    float lam = lon*((float)(M_PI/180.0));
+    float sphi = sinf(phi);
+    float cphi = cosf(phi);
+    float slam = sinf(lam);
+    float clam = cosf(lam);
+    float x1 = clam*x + slam*y;
+    float north = -sphi*x1 + cphi*z;
+    float east = -slam*x + clam*y;
+    float down = -cphi*x1 + -sphi*z;
+    float horizontal = sqrtf(north*north + east*east);
+    float total = sqrtf(horizontal*horizontal + down*down);
+    float inclination = atan2f(down, horizontal)*((float)(180.0/M_PI));
+    float declination = atan2f(east, north)*((float)(180.0/M_PI));
+    return {north, east, down, horizontal, total, inclination, declination};
+}
+
+
 /** Return the position in International Terrestrial Reference System coordinates, units meters.
 Using the WGS 84 ellipsoid and the algorithm from https://geographiclib.sourceforge.io/
  INPUT: