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.