((lat2-lat1)*BETA*(lat2-lat1)*BETA + (long2-long1)*delta*(long2-long1)*delta)
I'm writing some C code and I am getting some screwy results. Basically lat1, lat2, long1, long2 are all latitudes and longitudes in radians, delta is another precalculated parameter. BETA is a defined constant.
As far as I can tell this code always results in an answer equal to BETA regardless of changing inputs. Additionally when I brake out the two addition terms: (lat2-lat1)*BETA*(lat2-lat1)*BETA and (long2-long1)*delta*(long2-long1)*delta), and compute them separately I get a result of "0". I have been trying to wrap my head around what is wrong without success for a day now. Any suggestions for things to try would be appreciated.
Thanks in advance.
But the big thing here is of course that earth is not flat, so pythagoras is no fun to use except for very small areas - unless the goal is to bild a straight tunnel (which would feel like a downhill slope followed by an uphill slope if really long).
The wikipedia entry for evaluating the great-circle distance is: en.wikipedia.org/.../Great-circle_distance
But the simplification of computing using an averaged earth radius instead of caring about earth being oblongated will still give _way_ better results than the naive pythagoras.
What quality that is needed comes back to the original question: what is the evaluation intended to be used for?
Note that a PC using single-precision arithmetic will not compute in the same way that the Keil code does. When Keil multiply two float values they do it in single precision. A PC will have extended precision for intermediate operations. Most PC compilers then have compilation flags to control if the generated code should follow the IEEE regulations for storing back intermediate results into limited-sized variables or if the compiler should keep the intermediate values in the floating-point registers (which normally are 80-bit large).
But in your case, you will get very large numbers because you do not divide with 2*pi before you square.
earth is not flat, so pythagoras is no fun
It's not quite that bad. As I said, the formula uses Pythagoras on arc lengths instead of straight lines (secants). I.e. roads, not tunnels.
But yes, a less incorrect function would have to go from arcs to secants, use Pythagoras, then convert the resulting secant back into a surface arc.
That's got nothing to do with the problem. Given his input coordinates are in radians, dividing by 2*pi would be correct only if accompanied by multiplying the result by (2*pi)^2 afterwards --- a waste of effort for nothing unless intermediate results were in danger of hitting overflow. Which they aren't.
If one really wanted to reduce number scales, it's the earth's radius in meters that should be factored out:
ecc=ALPHA/BETA # a constant
arc_lat = lat1 - lat2 arc_long = ecc*cos(lat1)*(long1-long2)
result_angle = hypot(arc_lat, arc_long)
result_meters = result_angle*BETA
Thanks all for the help, I believe I have resolved this with your help.
arc_lat = lat1 - lat2; part1 = cos(lat1); part2 = long1-long2; arc_long = ECC*part1*part2; part1 = pow(arc_lat,2); part2 = pow(arc_long,2); return BETA * sqrt( part1+part2 );
Without using pow() I could not get valid results for computing the hypotenuse...
This algorithm is used for short distances. For that reason these assumptions produce what I would call an acceptable error (1 to 2 meters of error for a 500m measurement). While more elegant solutions are possible, it does not mean they are necessary. If I ever did find a need to squeeze more accuracy out I could always use more locally defined radii for the specific area of interest.
Thank you for your suggestions, regardless.
So are you telling us that not only did you observe a difference between pow(x,2) and x*x here, but pow(x,2) was closer to correct? Now that would be seriously worrying. If there's a difference at all, it's pow() that should be less precise than simple multiplication.
What were the actual results for both methods?
What I mean by invalid results was the same supposed cancellation problem we were discussing before. x*x resulted in a 0 answer where pow(x,2) gave me an acceptable result.