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.
No comments:
Post a Comment