ClosestPointDistanceToClosestPointOnSegment fails for some pairs on globe... see video and code.. any ideas ?

[SkybuckFlying]

See picture in second comment for problem.

See here for example pairs fail, near end of video:

https://www.youtube.com/live/-EsixpO9mBk?feature=share

// engineer who provided it to you for assistance.
// Define some earth constants
const MAJA = (6378137.0); // semi major axis of ref ellipsoid
const FLAT = (1.0/298.2572235630); // flattening coef of ref ellipsoid.
// These are derived values from the above WGS-84 values
const ESQR = (FLAT * (2.0-FLAT)); // 1st eccentricity squared
const OMES = (1.0 - ESQR); // 1 minus eccentricity squared
const EFOR = (ESQR * ESQR); // Sqr of the 1st eccentricity squared
const ASQR = (MAJA * MAJA); // semi major axis squared := nlmaja**

procedure ConvertECEFToLTP
(
const ParaEcef : TEcef;
var ParaLtp : TLtp
);
var
a0,a1,a2,a3,a4,b0,b1,b2,b3 : double;
b,c0,opqk,qk,qkc,qks,f,fprm : double;
k : integer;
ytemp, xtemp : double;
begin
// Convert from ParaEcef (XYZ) to LTP (ParaLtp.Lat/ParaLtp.Lon/ParaLtp.Alt)

//-------------------------------------------
// b := (0*0 + 1*1) / semi_major_axis_squared
// c := (z*z) / semi_major_axis_squared
//-------------------------------------------
b :=
	(
		(ParaEcef.x * ParaEcef.x) +
		(ParaEcef.y * ParaEcef.y)
	) / ASQR;
c0 :=
	(ParaEcef.z * ParaEcef.z) / ASQR;

//-------------------------------------------
// a0 := c * one_minus_eccentricity_sqr
// a1 := 2 * a0
// a2 := a0 + b - first_eccentricity_to_fourth
// a3 := -2.0 * first_eccentricity_to_fourth
// a4 := - first_eccentricity_to_fourth
//-------------------------------------------
a0 := OMES * c0;
a1 := 2.0 * a0;
a2 := a0 + b - EFOR;
a3 := -2.0 * EFOR;
a4 := - EFOR;

//-------------------------------------------
// b0 := 4 * a0, b1 := 3 * a1
// b2 := 2 * a2, b3 := a3
//-------------------------------------------
b0 := 4.0 * a0;
b1 := 3.0 * a1;
b2 := 2.0 * a2;
b3 := a3;

//-------------------------------------------
// Compute First Eccentricity Squared
//-------------------------------------------
qk := ESQR;

// for (k := 1; k <:= 3; k++)
for k := 1 to 3 do
begin
qks := qk * qk;
qkc := qk * qks;
f :=
(a0 * qks * qks) +
(a1 * qkc) +
(a2 * qks) +
(a3 * qk) + a4;
fprm :=
(b0 * qkc) +
(b1 * qks) +
(b2 * qk) + b3;
qk := qk - (f / fprm);
end;

//-------------------------------------------
// Compute latitude, longitude, altitude
//-------------------------------------------
opqk := 1.0 + qk;
if
(
	(ParaEcef.x = 0.0) and
	(ParaEcef.y = 0.0)
) then
// on the earth's axis
begin
	// We are sitting EXACTLY on the earth's axis
	// Probably at the center or on or near one of the poles
	ParaLtp.Lon := 0.0; // as good as any other value
	if (ParaEcef.z >= 0.0) then
	begin
		ParaLtp.Lat := PI / 2; // ParaLtp.Alt above north pole
	end	else
	begin
		ParaLtp.Lat := -PI / 2; // ParaLtp.Alt above south pole
	end;
end else
begin
	ytemp := opqk * ParaEcef.z;
	xtemp := sqrt
	(
		( ParaEcef.x * ParaEcef.x ) +
		( ParaEcef.y * ParaEcef.y )
	);
	ParaLtp.Lat := arctan2( ytemp, xtemp );
	ParaLtp.Lon := arctan2( ParaEcef.y, ParaEcef.x );
end;
ParaLtp.Alt :=
	(
		1.0 - OMES / ESQR * qk
	) * MAJA *
	sqrt
	(
		(
			b / (opqk * opqk)
		) + c0
	);

end; // ConvertECEFToLTP()

// written by Skybuck Flying ! ;) :) =D
function ClosestPointDistanceToClosestPointOnSegment
(
ParaPointX, ParaPointY, ParaPointZ : double;
ParaStartX, ParaStartY, ParaStartZ : double;
ParaEndX, ParaEndY, ParaEndZ : double;
ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres).
) : double;
var
vPoint,
vSegmentPoint1,
vSegmentPoint2,
vClosestPointOnSegment : TLatLonNvectorSpherical;
begin
vSegmentPoint1.FromXYZ2( ParaStartX, ParaStartY, ParaStartZ, ParaRadius );
vSegmentPoint2.FromXYZ2( ParaEndX, ParaEndY, ParaEndZ, ParaRadius );
vPoint.FromXYZ2( ParaPointX, ParaPointY, ParaPointZ, ParaRadius );

vClosestPointOnSegment := nearestPointOnSegment
(
	vPoint,
	vSegmentPoint1,
	vSegmentPoint2,
	ParaRadius
);

result := distanceTo
(
	vPoint, vClosestPointOnSegment, ParaRadius
);

end;

function ClosestPointDistanceToClosestPointOnSegmentV2
(
ParaPointX, ParaPointY, ParaPointZ : double;
ParaStartX, ParaStartY, ParaStartZ : double;
ParaEndX, ParaEndY, ParaEndZ : double;
ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres).
) : double;
var
vPoint,
vSegmentPoint1,
vSegmentPoint2,
vClosestPointOnSegment : TLatLonNvectorSpherical;
vEcef : Tecef;
vLtp : Tltp;
begin
// convert start x,y,z to lat lng
vEcef.X := ParaStartX;
vEcef.Y := ParaStartY;
vEcef.Z := ParaStartZ;

ConvertECEFToLTP
(
	vEcef,
	vLtp
);

vSegmentPoint1.mLat := vLtp.Lat;
vSegmentPoint1.mLon := vLtp.Lon;

// convert end x,y,z to lat lang
vEcef.X := ParaEndX;
vEcef.Y := ParaEndY;
vEcef.Z := ParaEndZ;

ConvertECEFToLTP
(
	vEcef,
	vLtp
);

vSegmentPoint2.mLat := vLtp.Lat;
vSegmentPoint2.mLon := vLtp.Lon;

// ParaRadius ??

// convert point x,y,z to lat lng
vEcef.X := ParaPointX;
vEcef.Y := ParaPointY;
vEcef.Z := ParaPointZ;

ConvertECEFToLTP
(
	vEcef,
	vLtp
);

vPoint.mLat := vLtp.Lat;
vPoint.mLon := vLtp.Lon;

vClosestPointOnSegment := nearestPointOnSegment
(
	vPoint,
	vSegmentPoint1,
	vSegmentPoint2,
	ParaRadius
);

result := distanceTo
(
	vPoint, vClosestPointOnSegment, ParaRadius
);

end;

Any idea what could be wrong ? If you need more code let me know...

[SkybuckFlying]