Approximating Magnetic Declination Using a Lookup Table

[scottfromscott]

I have a Baby Orangutan robot controller with an AVR 328B microcontroller. I have an HM55B compass module attached and wanted to be able to adjust the compass readings for true north. There is software at the NOAA web site to do this, but I was unable to figure out how to make it work with my controller. I decided to use the NOAA software to calculate magnetic declinations for points around the globe, then write the software to approximate declination values in between those points. The following C functions accomplish this within a few degrees of accuracy...

Code: [Select]

#include <stdio.h>

/*
Adjust electronic compass readings from magnetic north to true north using a subset of data generated by the NOAA WMM70.exe software.
Approximate magnetic declination for latitudes -60..60 and all longitudes for July 15, 2012 (the midpoint of the the period covered by
the model.) Points at 10 degree intervals of a (latitude, longitude) grid are plotted for declination with the WMM2010 software.
Declinations are stored in the source program in a literal string. Bilinear interpolation is then used to approximate the declinations
within these 10 degree grids. This method allows microcontrollers with minimal memory resources to correct readings from a compass sensor
for true north. Test readings should be within a few degrees of the WMM2010 model for most inhabited locations.

If you just need to use the existing lookup table here for use in your program, you will need only the following:
sc() function to decode declination values from encoded ascii to integer format
declination() function to retrieve values from the lookup table dec_tbl[]
Call declination with your latitude and longitude. approximation for declination is returned.

To create your own lookup table:

STEP 1: download and install geomag70.exe from [url=http://www.ngdc.noaa.gov/geomag/WMM/soft.shtml]www.ngdc.noaa.gov/geomag/WMM/soft.shtml[/url]

STEP 2: read the readme.txt file
eg. run the program at the command prompt for help: geomag70 h
then modify fprintf() format string in the build_wmm_input_file() for your required values
run the build_wmm_input_file() function to create geomag70_input.txt
the file contents will look something like this:
 
2012.5 D M100 -60 -180
2012.5 D M100 -60 -170
:
2012.5 D M100 60 170
2012.5 D M100 60 180
 
STEP 3: To create the raw output file run: geomage70.exe  WMM2010.COF f geomag70_input.txt geomag70_output.txt
the file contents will look something like this:
 
Date Coord-System Altitude Latitude Longitude D_deg D_min I_deg I_min H_nT X_nT Y_nT Z_nT F_nT dD_min dI_min dH_nT dX_nT dY_nT dZ_nT dF_nT
2012.5 D M100 -60 -180    46d 37m   -77d 45m   13243.3   9095.8   9625.6 -60953.3  62375.4     6.6       1.2         10.5    -11.2     25.0     55.3    -51.8
2012.5 D M100 -60 -170    45d 31m   -75d 45m   14906.9  10446.2  10634.5 -58725.1  60587.5     5.8       1.4          8.1    -12.3     23.3     66.1    -62.1
:
2012.5 D M100 60 170    -1d 28m    70d 21m   18169.1  18163.1   -467.1  50891.8  54037.8    -2.2       1.9         -6.5     -6.8    -11.6     68.4     62.3
2012.5 D M100 60 180     3d 58m    70d 23m   18041.0  17997.9   1247.3  50595.9  53716.1    -6.4       1.8         -5.0     -2.7    -33.6     69.0     63.3
 
STEP 4: To build the lookup table run build_lookup_table()
this will extract the declination values 99d 99m from the 6th and 7th columns in the geomag70_output.txt file for each lat,lon
to create the file wmm_ss.txt containing encoded declinations in ascii text format
the file contents will look something like this:

}}|{zxwuqm .. Z]`a`]XPICABEIOT

STEP 5: using a text editor, copy the ascii data into the declination() function source code to create the dec_tbl[] char array
to make your lookup table

STEP 6: test the declination() function by calling it with various latitudes and longitudes
compare the approximated declination the geomag_70.exe calculated declination

[email protected]
*/

void build_wmm_input_file()
{
/*
Build input file for processing by the geomag70.exe program available from the NOAA web site above
The input file has the following format
where columns are date, coord_system, meters_above_sea_level, lat, lon
2012.5 is midpoint of WMM2010 model's range
100 meters above sea level is is the average elevation of populated areas

2012.5 D M100 -60 -180
2012.5 D M100 -60 -170
:
2012.5 D M100 60 170
2012.5 D M100 60 180
*/

int lat_index, lon_index;

FILE *fp;

fp = fopen("geomag70_input.txt", "w");

for (lat_index=-60;lat_index<=60;lat_index += 10)
{
for(lon_index=-180;lon_index<=180;lon_index += 10)
{
fprintf(fp,"2012.5 D M100 %hd %hd\n",lat_index,lon_index);
}
}
fclose(fp);
}


char sc(int x, char c)
{
/*
a 'short char' for representing numbers in range -45..45 in displayable character strings
-45 maps to '#'(ascii 35) 0 maps to 'P' (ascii 80) and 45 maps to '}' (ascii 125)
'\' (ascii 92) is mapped to '!' (ascii 126) to avoid unintended escape sequences
store declinations in a string literal in the source program
values for char c are 'e'ncode or 'd'ecode
declinations out of range are truncated to +|- limits (very few qualify)
*/

if (c=='e') return (x  < 0)? ((x < -45)? 35 : x + 80) : ((x > 45)? 125 :(x == 12)? 126 : x + 80); // short int -> short char
else return (x < 80)? ((x < 35)? -45 : x - 80) :((x == 126)? 12 : x - 80); // short char -> short int
}

void build_lookup_table()
{
/*
for each lat,lon, read declination value from geomag70_output.txt file
write declination value formatted as encoded ascii text to the wmm_ss.txt file
*/

FILE *inf, *outf;

float d, m;
int dec, i, j, k, lat_index, lon_index;
char line[160], dec_deg[4], dec_min[4];

inf = fopen("geomag70_output.txt", "r");
outf = fopen("wmm_ss.txt", "w");

fgets ( line, sizeof(line), inf );

k = 0;

for (lat_index=-60;lat_index<=60;lat_index += 10)
{
for(lon_index=-180;lon_index<=180;lon_index += 10)
{
// get declination from file

fgets ( line, sizeof(line), inf );

// dig out the degrees and minutes

i = 0;
for(j=1;j<=5;j++)
{
while(line[i]!=' ') i++;
i++;
}
i+=2;
j = 0;
while (line[i]!='d') dec_deg[j++] = line[i++];
dec_deg[j]='\0';
i+=2;
j=0;
while (line[i]!='m') dec_min[j++] = line[i++];
dec_min[j]='\0';
d = atoi(dec_deg);
m = atoi(dec_min);
m=(d<0)? -m : m;

dec = round(d+m/60);

//output to clean file...

fputc(sc(dec,'e'),outf);

//debug...
k++;
printf("%i: %i %i %i \n", k, lat_index, lon_index, dec);
}
}
fclose(inf);
fclose(outf);
}

int declination(int lon, int lat)
{
/* lookup declination values from 10 x 10 degree grid and return approximate declination for (lon, lat) */

char dec_tbl[482] = \
"}}|{zxwuqmgaZTOJFB<6/'#########7fx}}}nonnnnnmlhb[TLGCA>;60)#" \
"######0AQ]ejmnffgffffgfc^VME?<:9973-($#$(0:DMTZ_befaaaaaaaa`" \
"^XPG?:8679;;830038?FLQTX[^_a]]]^^]]~~YTLC<878:>CFFB?>?CGMPRT" \
"VX[~][[[[[[ZZYVQIA<9:<@EJMNLIGGILOQRSUWYZ[YYYYYYYYWTNGA=<>BF" \
"JMPQPNLLLNPQQRSUWYYYYYYYYYXVRLFA>>AFJMOQRRPONNOOPOPQSVXYXYYZ" \
"[[ZYVQKEA?ADHLOQRSSRPPOOONMMNPSVXWY[~]]~ZVPIC@@BEIMPQSTTSRRQ" \
"POMKJKMPTWUY~^__^[VOGA??BEILPRTUVVUUSQOKIGHJMQUTY]`bba]VMD>;" \
"<?CGKORTWYZZYWTOJFDDGKOTTY^befd_UH=757:?DIMRVZ]`a`]XPICABEIO" \
"T";

float dec00, dec01, dec11, dec10, decmin, decmax;
float lonmin,latmin;

/* set base point (latmin, lonmin) of grid */

lonmin = (lon != 180)? (int)(lon / 10) * 10: 170;
if (lat >= 60) latmin = 50;
else if (lat < -60) latmin = -60;
else latmin = (int)(lat / 10) * 10;

dec00 = sc(dec_tbl[(int)((60 + (latmin +  0)) * 3.7 + (180 + (lonmin +  0)) / 10)],'d');
dec01 = sc(dec_tbl[(int)((60 + (latmin + 10)) * 3.7 + (180 + (lonmin +  0)) / 10)],'d');
dec11 = sc(dec_tbl[(int)((60 + (latmin + 10)) * 3.7 + (180 + (lonmin + 10)) / 10)],'d');
dec10 = sc(dec_tbl[(int)((60 + (latmin +  0)) * 3.7 + (180 + (lonmin + 10)) / 10)],'d');

decmin = (float)(lon - lonmin) / 10 * (dec10 - dec00) + dec00;
decmax = (float)(lon - lonmin) / 10 * (dec11 - dec01) + dec01;

return round((float)(lat - latmin) / 10 * (decmax - decmin) + decmin); // by bilinear interpolation
}

int main()
{
/*
for lon -175, lat -55, this program returns a declination of 36...

> ch grab60.c
> Lon: -175, Lat: -55, Dec: 36

geomag70.exe returns 36d 50m...

> geomag70 WMM2010.COF 2012.5 D M100 -55 -175

Geomag v7.0 - Jan 25, 2010

Model: WMM2010
Latitude: -55.00 deg
Longitude: -175.00 deg
Altitude: 100.00 meters
Date of Interest:  2012.50
-------------------------------------------------------------------------------
   Date          D           I           H        X        Y        Z        F
   (yr)      (deg min)   (deg min)     (nT)     (nT)     (nT)     (nT)     (nT)
2012.50    36d  50m   -73d  54m   16595.5  13283.4   9948.0 -57468.9  59817.1
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Date         dD          dI           dH       dX       dY       dZ       dF
   (yr)      (min/yr)    (min/yr)    (nT/yr)  (nT/yr)  (nT/yr)  (nT/yr)  (nT/yr)
2012.50       5.3        0.8         -1.4    -16.4     19.4     55.1    -53.3
-------------------------------------------------------------------------------
*/

int lat = -55, lon = -175;

//build_wmm_input_file();
//build_lookup_table();

printf("Lon: %i, Lat: %i, Dec: %i\n", lon, lat, declination(lon, lat));

return 0;
}

[scottfromscott]

October 1. 2010
The attached file contains version 3 of the 'Declination Table Development Kit' and replaces the code of the original posting above.
Changes made were:

• Menu to assist building a custom table or to lookup declinations from entered latitudes and longitudes
  
• improved table storage allowing declinations from -128 to +127 using one byte each