Geoid

Recently I was testing Matlab2012a geoid functions of which there are two: geoidheight (in Aerospace Toolbox) and egm96geoid (In Mapping Toolbox). I wrote the following script to check if there is a difference between these two function results for the same inputs (takes about half an hour to run):

clc; clear all;
[NVec,refvec] = egm96geoid(1, [-89.75  89.75],[-180 180]);
i = 0;
for lat_deg = -89:89
    fprintf(1, 'lat_deg = %1.1d\n', lat_deg);
    for lon_deg = 0:360
        i = i+1;
        N1_m(i) = geoidheight(lat_deg, lon_deg); %#ok<*SAGROW>

        N2Linear_m(i) = ltln2val(NVec,refvec,lat_deg,lon_deg,'bilinear');
        N2BiCubic_m(i) = ltln2val(NVec,refvec,lat_deg,lon_deg,'bicubic'); %to check if interpolation method makes any difference
    end
end
diff1 = N1_m-N2Linear_m;
diff2 = N2BiCubic_m-N2Linear_m;
subplot(2,1,1)
plot(diff1)
ylabel('diff1')
grid on
subplot(2,1,2)
plot(diff2)
ylabel('diff2')
grid on

As you can see from the plot, the difference between the two functions (diff1) can be as high as 7m! I checked the results with the zip files from NASA EGM96 and found that geoidheight is within 0.01m. This means that there is something wrong with egm96geoid or ltln2val.

Conclusion: Use geoidheight to get geoid offsets. Do not assume that Matlab's functions are correct and ABT (Always be testing). This error seems to be fixed in Matlab 2014 and later versions.