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-18-2019, 03:49 PM   #16
dgatsoulis
ele2png user
 
dgatsoulis's Avatar
Default

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.
dgatsoulis is offline   Reply With Quote
Thanked by:
Old 09-18-2019, 09:07 PM   #17
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
Originally Posted by martins View Post
 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.
JMW is offline   Reply With Quote
Old 09-19-2019, 04:03 PM   #18
martins
Orbiter Founder
Default

Quote:
Originally Posted by JMW View Post
 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?

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.
martins is offline   Reply With Quote
Old 09-20-2019, 01:56 PM   #19
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
Originally Posted by martins View Post
 and I do get the same results with both codes.. So this factor 0.11 is a mystery to me.
Now I'm confused
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?

Quote:
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....
Thanks for your patience.

Screenshots attached
Attached Thumbnails
Difference.png   No Difference.png  
JMW is offline   Reply With Quote
Old 09-20-2019, 03:14 PM   #20
martins
Orbiter Founder
Default

Quote:
Originally Posted by JMW View Post
 Do you mean you are getting the same results as me (a difference) using the two codes?
No.
Quote:
Or are you getting the same result with both codes?
Yes.

Quote:
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
Quote:
Can you explain a bit more please how I do this check....
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.
martins is offline   Reply With Quote
Old 09-20-2019, 07:15 PM   #21
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
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.
(Kourou --> Overberg - Best I could do)
Yours: 8811 km great circle distances
Mine: 9255 km (Only 5% diff now, not 11% ? )

True ?: 10,007.54339801029 km (but this is straight line distance isn't it ? ^ ^ ^ )
JMW is offline   Reply With Quote
Old 09-20-2019, 08:54 PM   #22
martins
Orbiter Founder
Default

For the great circle distance between Kourou (Long=-52.73, Lat=5.2) and Overberg (Long=+20.32, Lat=-34.58) I get dis=8.8052e3 with both our codes, which also agrees with the MFD display. This doesn't seem to match either of your results, but at least corresponds approximately to 8811. Maybe the vessel wasn't parked exactly on the base's reference point?

9255 is however significantly different. Can you explain in detail how you got to that value?
martins is offline   Reply With Quote
Old 09-20-2019, 09:09 PM   #23
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
9255 is however significantly different. Can you explain in detail how you got to that value?
It's just by compiling with "my code".......
Am I on the same page here ?

Quote:
I get dis=8.8052e3 with both our codes, which also agrees with the MFD display.
Don't understand how this happens.....

Hope I'm not leading you a merry dance here..........
JMW is offline   Reply With Quote
Old 09-20-2019, 10:00 PM   #24
martins
Orbiter Founder
Default

We'll get there eventually
Quote:
Originally Posted by JMW View Post
 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;
So, step by step, this is what I get (leaving out dis = (dis)*DEG;, which I hope we agree is incorrect):

Code:
latitude = 5.2
longitude = -52.73
tgtlat = -34.58
tgtlng = +20.32
dis = sin((latitude)*RAD) * sin((tgtlat)*RAD) + cos((latitude)*RAD) * cos((tgtlat)*RAD) * cos((longitude-tgtlng)*RAD) = 0.1876
dis = acos(dis) = 1.3821
dis = dis*6371 = 8.8052e3
whereas you ... ?
martins is offline   Reply With Quote
Old 09-20-2019, 11:19 PM   #25
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
(leaving out dis = (dis)*DEG;, which I hope we agree is incorrect):
Well, only way I can get significant reading is by including dis = (dis)*DEG;

I then get 9248.97, otherwise it's a silly 165

Sorry, but I've gotta hit the sack now but this is doing my head in, so we need to get to the bottom of it.
Thanks for all so far...
JMW is offline   Reply With Quote
Old 09-21-2019, 05:20 PM   #26
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Quote:
Well, only way I can get significant reading is by including dis = (dis)*DEG;
Ok, I'm awake now and have remedied my mistake.

Checks out to 8805.18 in my environment (don't have Matlab)
Attached Thumbnails
Check.png  
JMW is offline   Reply With Quote
Old 09-21-2019, 07:32 PM   #27
martins
Orbiter Founder
Default

Excellent. Sounds like this is all fixed now.
martins is offline   Reply With Quote
Old 09-21-2019, 07:47 PM   #28
JMW
Aspiring Addon Developer
 
JMW's Avatar
Default

Thank you so much for your patient help
JMW is offline   Reply With Quote
Thanked by:
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 07:59 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.