SDK Question Distance calc overrun ?

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
Hi People.

Quick question, probably not quick answer.....

On a distance calculation I'm getting a value of "-1.#IND" intermittently.

It's messing up an autopilot.

Is it an overrun of some kind.

If so, what do I need to do....?

Code used:
Code:
dis = sin((latitude)*RAD) * sin((tgtlat)*RAD) + cos((latitude)*RAD) * cos((tgtlat)*RAD) * cos((longitude-tgtlng)*RAD);
  dis = acos(dis);
  dis = (dis)*DEG;
  dis = dis * 6371;
 
Last edited:

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
36,879
Reaction score
1,537
Points
203
Location
Langendernbach
The difference in longitude should be absolute:
Code:
fabs(longitude-tgtlng)
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,407
Reaction score
341
Points
83
Website
orbit.medphys.ucl.ac.uk
Hi People.

Quick question, probably not quick answer.....

On a distance calculation I'm getting a value of "-1.#IND" intermittently.

It's messing up an autopilot.

Is it an overrun of some kind.

If so, what do I need to do....?

Code used:
Code:
dis = sin((latitude)*RAD) * sin((tgtlat)*RAD) + cos((latitude)*RAD) * cos((tgtlat)*RAD) * cos((longitude-tgtlng)*RAD);
  dis = acos(dis);
  dis = (dis)*DEG;
  dis = dis * 6371;

The only function that could potentially throw an error I can see here is the acos. Did you check that the input is valid (range [-1,1])? Even if your calculation of the argument is mathematically correct, there is always the possibility of rounding errors.
The difference in longitude should be absolute:
Code:
fabs(longitude-tgtlng)
Would that make a difference here, given that cos(-x) = cos(x)?
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
36,879
Reaction score
1,537
Points
203
Location
Langendernbach
Would that make a difference here, given that cos(-x) = cos(x)?


Not really, but the other terms should not exceed 2.0 in any case and not 1.0 if the shortest great circle distance is searched. Was just the first difference to the reference I noticed
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,407
Reaction score
341
Points
83
Website
orbit.medphys.ucl.ac.uk
Also, for the cases where it doesn't fail, does it give you the expected result? I am a bit concerned about those lines:

Code:
  dis = (dis)*DEG;
  dis = dis * 6371;
Is that making the approximation sin(x) ~ x (i.e. small distances)? Even then, converting the central angle distance to degrees rings alarm bells here ...
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
Gentlemen, sorry for the delay, have been checking a few things.

Urwumpe: I've tried 'fabs(longitude-tgtlng)' and it sent the calcs haywire.

martins:
By 'the input' do you mean the result of line 1 ?
That does seem to be within -1 && 1 always.

Yes, when it doesn't read "-1.#IND" it gives a perfect distance.

It's occurring whilst hovering directly over a pad, so very small distances are involved....

I got that formula from a website and also from a post in the forum, but can't remember whose.
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
36,879
Reaction score
1,537
Points
203
Location
Langendernbach
Urwumpe: I've tried 'fabs(longitude-tgtlng)' and it sent the calcs haywire.


That shouldn't happen at all, since it, as martins pointed out, should be just a superfluous term if cos is properly implemented.



"-1.#IND" is the windows version of telling you "NaN" or "not a number", which occurs, when the result of an floating point operation is undefined.

acos(1) and acos(-1) should be no NaN, but acos(1.0000000000000000000000000000000000000000000000000000000000000000000000000001) would be one.



Your formula is also mentioned in Wikipedia for great circle distance and in the math book I used at university, but always with an absolute definition of difference in latitude.
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
Urwumpe: Apologies, I didn't implement fabs correctly. Is fine now.

Marijn: I'll take a look there later, (RL taking over for now)

Sounds like the small distance might be the culprit.

Thanks all.

---------- Post added at 07:48 PM ---------- Previous post was at 07:30 PM ----------

martins:
What did you use for distance calc in VOR/VTOL mfd ? (if it's not classified:))

---------- Post added 17-09-19 at 05:48 PM ---------- Previous post was 16-09-19 at 07:48 PM ----------

martins: Can it be "open source" ? ^ ^ ^ ^ ^ :bighug:
I'm curious because there is about 0.22 km difference (mine greater)
at 2 km, reducing to 0.11 at 1 km, 0.05 at 500 mtrs, down to less than 0.09 at zero.
 

Urwumpe

Not funny anymore
Addon Developer
Donator
Joined
Feb 6, 2008
Messages
36,879
Reaction score
1,537
Points
203
Location
Langendernbach
I'm curious because there is about 0.22 km difference (mine greater)
at 2 km, reducing to 0.11 at 1 km, 0.05 at 500 mtrs, down to less than 0.09 at zero.


Which radius did you use and what is the true radius?
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
6371 km.........

---------- Post added at 06:42 PM ---------- Previous post was at 06:37 PM ----------

True... depends where you are ?
The radius of Earth at the equator is 3,963 miles (6,378 kilometers), according to NASA's Goddard Space Flight Center. However, Earth is not quite a sphere. The planet's rotation causes it to bulge at the equator. Earth's polar radius is 3,950 miles (6,356 km) — a difference of 13 miles (22 km).
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,407
Reaction score
341
Points
83
Website
orbit.medphys.ucl.ac.uk
What did you use for distance calc in VOR/VTOL mfd ? (if it's not classified:))

martins: Can it be "open source" ? ^ ^ ^ ^ ^ :bighug:

Let me just quickly declassify this for you. Here is the relevant VOR/VTOL MFD code:

Code:
void Orthodome (double lng1, double lat1, double lng2, double lat2,
				double &dist, double &dir)
{
	double A = lng2-lng1;
	double dlng = fabs(A);
	double dlat = fabs(lat2-lat1);
	if (dlat < 1e-14) {
		dist = dlng;
		dir = (lng2 > lng1 ? PI05 : 3*PI05);
		return;
	} else if (dlng < 1e-14) {
		dist = dlat;
		dir = (lat2 > lat1 ? 0.0 : PI);
	} else {
		double sinA  = sin(A),    cosA  = cos(A);
		double slat1 = sin(lat1), clat1 = cos(lat1);
		double slat2 = sin(lat2), clat2 = cos(lat2);
		double cosa  = slat2*slat1 + clat2*clat1*cosA;
		dist = acos (cosa);
		dir = atan2 (sinA*clat2, clat1*slat2 - slat1*clat2*cosA);
		if (dir < 0.0) dir += Pi2;     // 0 <= dir < 2pi
	}
}

Orthodome (sp->lng, sp->lat, tlng, tlat, adist, bdir);
bdir -= sp->dir;
if      (bdir <  0.0) bdir += Pi2;
else if (bdir >= Pi2) bdir -= Pi2;
dist = adist * sp->ref->Size();

Let me know if you spot a mistake.
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
A "mistake" ?
No radius !
We're now in sync :lol: :thumbup:
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
Dunno.

I ended up adding this

dist = dist *6371;
deleting
sp->adist * ref->Size();
 
Last edited:

dgatsoulis

ele2png user
Joined
Dec 2, 2009
Messages
1,846
Reaction score
142
Points
78
Location
Sparta
IIRC Earth's radius in Orbiter is 6371.01 km. Check Earth.cfg to make sure, or simply get the size of the reference body as Martin did.
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
So, was the discrepancy just down to a difference in radius? It sounded quite significant.

As I used the same old "dist = dist *6371;" and we ended in sync, I guess difference is due to something in the calcs.
Way above my head to spot though.

Thanks Doctor for your interest, appreciated.
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,407
Reaction score
341
Points
83
Website
orbit.medphys.ucl.ac.uk
As I used the same old "dist = dist *6371;" and we ended in sync, I guess difference is due to something in the calcs.

You mean we are now in sync because you replaced your own calculation with mine? :blink:

I am flattered, but as long as the 11% difference you reported isn't explained, there is no guarantee that my code is any better than yours, so there is still work to be done.

As far as I can tell, both our codes are functionally identical (except for the spurious line "dis = (dis)*DEG;" of course), and I do get the same results with both codes. So this factor 0.11 is a mystery to me.

Can you check which of the two results is the correct one, if any? (should be trivial for a suitable set of input values, e.g. (0,0) and (0,90)). Then at least we know where to look for the error.
 

JMW

Aspiring Addon Developer
Joined
Aug 5, 2008
Messages
580
Reaction score
40
Points
43
Location
Happy Wherever
and I do get the same results with both codes.. So this factor 0.11 is a mystery to me.

Now I'm confused :lol:
Just to clarify:
Do you mean you are getting the same results as me (a difference) using the two codes?
Or are you getting the same result with both codes?

When I use your code (plus substitution of "dis = dis *6371;" for "sp->ref->Size();" ) the readout agrees with mfd.
When I use my original code the 11% diff occurs.
So difference is definitely in the code isn't it?

Can you check which of the two results is the correct one, if any? (should be trivial for a suitable set of input values, e.g. (0,0) and (0,90))

Can you explain a bit more please how I do this check....:shifty:
Thanks for your patience.

Screenshots attached
 

Attachments

  • Difference.png
    Difference.png
    315.6 KB · Views: 10
  • No Difference.png
    No Difference.png
    367.3 KB · Views: 9

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,407
Reaction score
341
Points
83
Website
orbit.medphys.ucl.ac.uk
Do you mean you are getting the same results as me (a difference) using the two codes?
No.
Or are you getting the same result with both codes?
Yes.

When I use your code (plus substitution of "dis = dis *6371;" for "sp->ref->Size();" ) the readout agrees with mfd.
When I use my original code the 11% diff occurs.
So difference is definitely in the code isn't it?
Not as far as I can tell. I quickly copied both our codes into Matlab to compare, and I got the same results for some arbitrary position pairs (not an exhaustive proof, but then, I can't really see how they could be different - It really is the same formula in both cases).

Here is the Matlab code (jmw is your code, ms is mine. {mylat, mylng} and {tgtlat, tgtlng} are two random positions. dis1 and dis2 are the great circle distances with your and my code.) As result I get dis1 = dis2 = 4.6836e+03.
Code:
function compare_distances
RAD = pi/180.0;
DEG = 1/RAD;
sz = 6371;

mylat = 23;
mylng = 48;
tgtlat = 55;
tgtlng = 11;

dis1 = jmw(mylat, mylng, tgtlat, tgtlng);
dis2 = ms(mylng*RAD, mylat*RAD, tgtlng*RAD, tgtlat*RAD);

function dis = jmw(lat0, lng0, lat1, lng1)
  dis = sin((lat0)*RAD) * sin((lat1)*RAD) + cos((lat0)*RAD) * cos((lat1)*RAD) * cos((lng0-lng1)*RAD);
  dis = acos(dis);
  %dis = (dis)*DEG;
  dis = dis * sz;
end

function dist = ms(lng1, lat1, lng2, lat2)
	A = lng2-lng1;
	dlng = abs(A);
	dlat = abs(lat2-lat1);
	if dlat < 1e-14
		dist = dlng;
        if lng2 > lng1
            dir = pi/2;
        else
            dir = 3*pi/2;
        end
    elseif dlng < 1e-14
		dist = dlat;
        if lat2 > lat1
            dir = 0;
        else
            dir = pi;
        end
    else
		sinA  = sin(A),    cosA  = cos(A);
		slat1 = sin(lat1), clat1 = cos(lat1);
		slat2 = sin(lat2), clat2 = cos(lat2);
		cosa  = slat2*slat1 + clat2*clat1*cosA;
		dist = acos (cosa);
		dir = atan2 (sinA*clat2, clat1*slat2 - slat1*clat2*cosA);
		if dir < 0.0
            dir = dir + pi*2;     % 0 <= dir < 2pi
        end
    end
    dist = dist*sz;
end

end

Can you explain a bit more please how I do this check....:shifty:
Thanks for your patience.
If you pick two points 90 degrees apart, you know the distance: dis = 2*pi*rad/4. So if the two codes give you different results for this case, at most one can be correct, and you can validate directly with the known truth.
 
Top