Orbiter-Forum  

Go Back   Orbiter-Forum > Orbiter Space Flight Simulator > Orbiter SDK
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

Orbiter SDK Orbiter software developers post your questions and answers about the SDK, the API interface, LUA, meshing, etc.

Reply
 
Thread Tools
Old 09-16-2019, 03:53 PM   #1
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default Distance calc overrun ?

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 by JMW; 09-16-2019 at 03:56 PM. Reason: include code
JMW is offline   Reply With Quote
Old 09-16-2019, 04:18 PM   #2
Urwumpe
Certain Super User
 
Urwumpe's Avatar

Default

The difference in longitude should be absolute:
Code:
fabs(longitude-tgtlng)
Urwumpe is offline   Reply With Quote
Thanked by:
Old 09-16-2019, 04:29 PM   #3
martins
Orbiter Founder
Default

Quote:
Originally Posted by JMW View Post
 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.
Quote:
Originally Posted by Urwumpe View Post
 The difference in longitude should be absolute:
Code:
fabs(longitude-tgtlng)
Would that make a difference here, given that cos(-x) = cos(x)?
martins is offline   Reply With Quote
Thanked by:
Old 09-16-2019, 05:09 PM   #4
Urwumpe
Certain Super User
 
Urwumpe's Avatar

Default

Quote:
Originally Posted by martins View Post
  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
Urwumpe is offline   Reply With Quote
Old 09-16-2019, 05:26 PM   #5
martins
Orbiter Founder
Default

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 ...
martins is offline   Reply With Quote
Old 09-16-2019, 05:40 PM   #6
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

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.
JMW is offline   Reply With Quote
Old 09-16-2019, 05:56 PM   #7
Urwumpe
Certain Super User
 
Urwumpe's Avatar

Default

Quote:
Originally Posted by JMW View Post
 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.0000000000000000000000000000000000000000000 000000000000000000000000000000001) 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.
Urwumpe is offline   Reply With Quote
Thanked by:
Old 09-16-2019, 06:05 PM   #8
Marijn
Orbinaut
 
Marijn's Avatar
Default

Have a look at this site with formulas. Under the head "Spherical Law of Cosines", it says that you might need a different formula for very small distances: http://www.movable-type.co.uk/scripts/latlong.html
Marijn is offline   Reply With Quote
Thanked by:
Old 09-17-2019, 04:48 PM   #9
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

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" ? ^ ^ ^ ^ ^
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.
JMW is offline   Reply With Quote
Old 09-17-2019, 04:54 PM   #10
Urwumpe
Certain Super User
 
Urwumpe's Avatar

Default

Quote:
Originally Posted by JMW View Post
  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?
Urwumpe is offline   Reply With Quote
Old 09-17-2019, 05:42 PM   #11
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

6371 km.........

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

True... depends where you are ?
Quote:
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).
JMW is offline   Reply With Quote
Old 09-17-2019, 08:15 PM   #12
martins
Orbiter Founder
Default

Quote:
Originally Posted by JMW View Post
 What did you use for distance calc in VOR/VTOL mfd ? (if it's not classified)

martins: Can it be "open source" ? ^ ^ ^ ^ ^
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.
martins is offline   Reply With Quote
Thanked by:
Old 09-17-2019, 10:57 PM   #13
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

A "mistake" ?
No radius !
We're now in sync
JMW is offline   Reply With Quote
Old 09-18-2019, 11:04 AM   #14
martins
Orbiter Founder
Default

So, was the discrepancy just down to a difference in radius? It sounded quite significant.
martins is offline   Reply With Quote
Old 09-18-2019, 02:36 PM   #15
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Dunno.

I ended up adding this

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

Last edited by JMW; 09-18-2019 at 04:03 PM.
JMW is offline   Reply With Quote
Reply

  Orbiter-Forum > Orbiter Space Flight Simulator > Orbiter SDK


Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump


All times are GMT. The time now is 01:18 AM.

Quick Links Need Help?


About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2019, vBulletin Solutions Inc.
Copyright 2007 - 2017, Orbiter-Forum.com. All rights reserved.