geoidh.m

GPS TOOLBOX包含以下内容： 1、GPS相关常量和转换因子； 2、角度变换； 3、坐标系转换： &#61656 点变换； &#61656 矩阵变换； &#61656 向量变换

%                             geoidh.m
%  Scope:   This MATLAB macro computes the WGS-84 geoid height correction for
%           the current location specified by latitude and longitude.
%  Usage:   xoutput = geoidh(lat,lon,tgeoid)
%  Description of parameters:
%           lat     - input, latitude in radians,  - pi/2 to pi/2
%           lon     - input, longitude in radians,  0 to 2*pi (automatically
%                     adjustable  -2*pi to 4*pi)
%           tgeoid  - input, geoid heights table (in meters) for a 19 by 36
%                     Latitude/Longitude grid; the latitudes go northward 
%                     from -90 to 90 degrees, and longitudes go eastward 
%                     from 0 to 350 degrees
%           xoutput - output, geoid height correction in meters
%  Reference:
%           [1] Anonym., Department of Defense World Geodetic System 1984.
%               DMA Technical Report 8350.2, Second Edition, September 1, 1991.
%  Last update:  04/23/00
%  Copyright (C) 1997-00 by LL Consulting. All Rights Reserved. 

function   xoutput = geoidh(lat,lon,tgeoid)

% Locate the grid for the specified location (latitude/longitude)
  
lat10 = (lat * (18.0 / pi)) + 10.;
phi1 = floor(lat10);
y = lat10 - phi1;
phi2 = phi1 + 1;
if (phi2 > 19)
  phi2 = 19;            %  adjustment for North pole only
end

lon10 = lon * (18.0 / pi) + 1.;
if (lon10 < 0)
  lon10 = lon10 + 36.0; % adjustment for negative longitude (up to -360 degrees)
end
lambda1 = floor(lon10);
x = lon10 - lambda1;
lambda2 = lambda1 + 1;
if (lambda1 > 36)
  lambda1 = lambda1 - 36;
end
if (lambda2 > 36)
  lambda2 = lambda2 - 36;
end

% Determine the interpolation weighting factors; xp = relative longitude 
% location, 0 <= xp < 1, and yp = relative latitude location, 0 <= yp < 1 

onemx = 1 - x;
onemy = 1 - y;

xp = x;
yp = y;
xpyp = xp*yp;
gweight(1) = (xpyp * xpyp * (9 - (6 * (xp + yp)) + 4 * xpyp));

xp = onemx;
yp = y;
xpyp = xp*yp;
gweight(2) = (xpyp * xpyp * (9 - (6 * (xp + yp)) + 4 * xpyp));

xp = onemx;
yp = onemy;
xpyp = xp*yp;
gweight(3) = (xpyp * xpyp * (9 - (6 * (xp + yp)) + 4 * xpyp));

xp = x;
yp = onemy;
xpyp = xp*yp;
gweight(4) = (xpyp * xpyp * (9 - (6 * (xp + yp)) + 4 * xpyp));

% Determine the interpolated value for geoid height

xoutput =  gweight(1) * tgeoid(phi2,lambda2) + ...
           gweight(2) * tgeoid(phi2,lambda1) + ...
           gweight(3) * tgeoid(phi1,lambda1) + ...
           gweight(4) * tgeoid(phi1,lambda2);