There is more to coding mathematical expressions than that of translating them into C language. The email excerpts below prove that point.

 

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.