function [deltaX, deltaY, deltaZ] = ecefOffsetFormula( ...
    spheroid, phi1, lambda1, h1, phi2, lambda2, h2, inDegrees)
%Cartesian ECEF offset between geodetic positions
%
%   The outputs are the same as in ecefOffset. The inputs are also the
%   same, except that the angleUnit string is replaced by logical scalar,
%   inDegrees, which is true when phi1, lambda1, phi2, and lambda2 are in
%   degrees and false when they are in radians.

% Copyright 2019-2020 The MathWorks, Inc.

%#codegen
    
    if ~isobject(spheroid) || ~isscalar(spheroid)
        error(message('geodesy:spheroid:ExpectedScalarSpheroid',upper('spheroid')))
    end
    
    coder.ignoreCatch
    try %#ok<EMTC>
        a = spheroid.a;
        f = spheroid.Flattening;
    catch e
        if ~isa(spheroid,'map.geodesy.Spheroid')
            error(message('geodesy:spheroid:ExpectedScalarSpheroid',upper('spheroid')))
        else
            throwAsCaller(e)
        end
    end
    
    if inDegrees
        s1 = sind(phi1);
        c1 = cosd(phi1);
        
        s2 = sind(phi2);
        c2 = cosd(phi2);
        
        p1 = c1 .* cosd(lambda1);
        p2 = c2 .* cosd(lambda2);
        
        q1 = c1 .* sind(lambda1);
        q2 = c2 .* sind(lambda2);
    else
        s1 = sin(phi1);
        c1 = cos(phi1);
        
        s2 = sin(phi2);
        c2 = cos(phi2);
        
        p1 = c1 .* cos(lambda1);
        p2 = c2 .* cos(lambda2);
        
        q1 = c1 .* sin(lambda1);
        q2 = c2 .* sin(lambda2);
    end
    
    if f == 0
        % Spherical case
        deltaX = a * (p2 - p1) + (h2 .* p2 - h1 .* p1);
        deltaY = a * (q2 - q1) + (h2 .* q2 - h1 .* q1);
        deltaZ = a * (s2 - s1) + (h2 .* s2 - h1 .* s1);
    else
        % Oblate spheroid case
        e2 = f * (2 - f);
        w1 = 1 ./ sqrt(1 - e2 * s1.^2);
        w2 = 1 ./ sqrt(1 - e2 * s2.^2);
        deltaX =            a * (p2 .* w2 - p1 .* w1) + (h2 .* p2 - h1 .* p1);
        deltaY =            a * (q2 .* w2 - q1 .* w1) + (h2 .* q2 - h1 .* q1);
        deltaZ = (1 - e2) * a * (s2 .* w2 - s1 .* w1) + (h2 .* s2 - h1 .* s1);
    end
    
    % Implementation Notes from Rob Comer
    % -----------------------------------
    % In order to minimize numerical round off for points that are
    % closely-spaced relative to the dimensions of the reference
    % ellipsoid, use a sequence of computations that avoids
    % computing a small difference by subtracting two large
    % numbers.  Rather than convert each point to ECEF, then
    % subtract, note that each ECEF coordinate in the following
    % formulas for x, y, and z:
    %
    %     N  = a ./ sqrt(1 - e2 * sin(phi).^2);
    %
    %     x = (N + h) .* cos(phi) .* cos(lambda);
    %     y = (N + h) .* cos(phi) .* sin(lambda);
    %     z = (N*(1 - e2) + h) .* sin(phi);
    %
    % contains a term that is a multiple of a and a term that is a
    % multiple of h.  For example,
    %
    %  x = a * cos(phi) * cos(lambda) * w + h * cos(phi) * cos(lambda);
    %
    % where w = 1 ./ sqrt(1 - e2 * sin(phi).^2).  Constructing x this
    % way for both x1 and x2, then taking the difference and factoring
    % out a gives:
    %
    %  dx = a * (cos(phi2) * cos(lambda2) * w1 - cos(phi1) * cos(lambda1) * w1)
    %       + (h2 * cos(phi2) * cos(lambda2) - h1 * cos(phi1) * cos(lambda1))
end
