The relative error of a product is the sum of
the relative errors of the factors. When forced to work with simple
precision (4 bytes), each value carries a relative error of +/-0.0000001.
This doesn't seem like much, but let's consider what happens when every
number in the following expression carries such an error:

latitude_table[i]= 180.0 / PI * acos(sqrt((1
- cos(PI / 30.0)) / (1 - cos(2 * PI / (i + 1)))));

The errors of the cos(x) and the sin(x)
functions, are equal to the error in the argument dx multiplied by the
function values of sin(x) or cos(x) respectively, this meaning that they are
a smaller than the error in the angle expressed in radians, but for the square roots,
the error is greater than that of the argument when it is less than one (as
the case goes). Same goes for the acos(x).

Let us neglect the multiplying effect of the functions for the moment an just count the products. Functions though, are solved by successive approximations algorithm and/or tables, meaning that, at least, they introduce as mush error as if squaring the argument.

Even integer constants like 1 or 2 can not be
expressed exactly when promoted to floating point numbers, every single
product must be counted (division counts just as a product when propagating
errors).

The expression (1 - cos(PI / 30.0)) could be a
great contributor to the error.

cos(PI / 30.0) =
0.99452189536827333692269194498057

but the error of evaluating it in single
precision is +/-0.0000003

(1 - cos(PI / 30.0)) = -0.005478 +/-0.0000003

Relative error is 0.00006, this is 600 times the
contribution of a single product.

The error contribution of expression (1 - cos(2
* PI / (i + 1)), will be less for all i <59. For 25 degrees latitude (i =
53), the contribution to the error will be of 0.000044.

The 11 product operations will add to the total
too and that will bring it up to .0001 in a best case scenario. Actual
values have been seen farther off than that. In geodesic terms, this error
at 25 degrees is about
0.006 degrees or one sixth of a nautical mile (~1000 ft), this
accounts pretty well for the observed glitch.

If instead of making the table by calculating
the expression for each integer from 1 - 60, the table could go
hard coded
with values calculated using double precision in EXCEL. This way we can
bring down the above error to the tens of feet range (50 ft). Better yet, if
instead of one table, we split it in 2; a table of shorts for the degrees
and one of doubles for the fraction of degrees, then were are down to the
feet range (2 feet).

***************************************************************************************************

Before plunging into coding all this ECEF (*Earth
Centered Earth Fixed reference system*) stuff ,a task that I regard as a
fun one, I decided to do some math to estimate what kind of an error are
we introducing by considering a constant radius to the earth (the point
of going into ECEF is for dealing with a variable one).

The relative difference in earth radius from the
poles to the equator is 2.4 x 10^{-3}. In the range we run our
intruders, which is 190 miles, the change in earth radius is only 8.4x10^{-5}.
Ignoring such a change may introduce a maximum error in the Latitude or
Longitude of dynamic intruders of 0.15 minutes. This will have no
further implications than a very slight change in course of 0.05 seconds. I
don't think such a correction is necessary and more knowing that it comes
with a heavy toll in CPU usage.

This said, the spherical correction must be
introduced. Current algorithm neglects earth curvature and assumes nearness
to equator. this introduces errors in longitude greater than 10 minutes for
a 190 EW range when latitude goes above 20 degrees.

My proposition is to concentrate on that main
source of error and neglect the ellipsoid effect.