function [dn,de,du] = deltablh(b1,l1,h1,b2,l2,h2);
% function [dn,de,du] = deltablh(b1,l1,h1,b2,l2,h2);
% Get distance along surface
% Inputs:
%  b1 = latitude of point 1 in degrees
%  l1 = longitude of point 1 in degrees
%  h1 = height of point 1 in meters
%  b2 = latitude of point 2 in degrees
%  l2 = longitude of point 2 in degrees
%  h2 = height of point 2 in meters
% Outputs (all in meters) :
%  dn = delta north 
%  de = delta east
%  du = delta up
% 
% Created by Will Lentz
% Trimble Navigation Limited
% 
% Modified by Vivek Nadkarni
% December 20, 1999

% Convert to radians
b1r = b1*pi/180;  % latitude
b2r = b2*pi/180;
l1r = l1*pi/180;  % longitude
l2r = l2*pi/180;

% Define WGS84 Ellipsoid Parameters
a = 6378137;
f = 1 / 298.257223563;
e2 = 2*f - f^2;

% Compute Radii of Curvature
W = sqrt( 1 - e2 * sin((b1r+b2r)/2).^2 );
N = a ./ W;
M = a .* ( 1 - e2 ) ./ W.^3;

% Compute Metric Components
dn = M .* (b2r-b1r);
de = N .* cos((b1r+b2r)/2) .* (l2r-l1r);
du = h2-h1;

