C++ Question Strange maths results and right ascension problems

[asbjos]

So while making a MFD, I need to calculate the right ascension of a body with the Earth (or any other object for that sake) as reference.
I get an angle, which turns out to be wrong (I assume the RA in Info screen [Ctrl+I] is correct), but when I try to convert the angle from radians to degrees-minutes-secons, I get seconds raging from 0 to 3600. Of course, I assume I did a mistake, but can't find any. So I load the same code in Matlab, but there it works as intended.

I ran this code in Orbiter (language is C++):

for (int n = 1; n <= 101; n++)
{
	float Angle = -4 * PI - 8 * PI / 101 + n * 8 * PI / 101;

	int Degrees = floor(abs(Angle) * DEG);
	int Minutes = floor((abs(Angle) * DEG - Degrees) * 60);
	float Seconds = (abs(Angle) * DEG - Degrees - Minutes / 60) * 3600;

	char Sign [256];
	if (Angle < 0) // Add the sign of the angle
		sprintf(Sign, "-");
	else
		sprintf(Sign, "");

	sprintf(LogMsg, "%.3f = %s%id %im %.2fs", Angle, Sign, Degrees, Minutes, Seconds);
	oapiWriteLog(LogMsg);
}

Then in Matlab, I ran this:

for n = 1 : 101
	Angle = -4 * pi - 8 * pi / 101 + n * 8 * pi / 101;

    DEG = 180 / pi;
	Degrees = floor(abs(Angle) * DEG);
	Minutes = floor((abs(Angle) * DEG - Degrees) * 60);
	Seconds = (abs(Angle) * DEG - Degrees - Minutes / 60) * 3600;

	if (Angle < 0) % Add the sign of the angle
		Sign = '-';
	else
		Sign = '';
    end

	LogMsg = sprintf('%.3f = %s%id %im %.2fs', Angle, Sign, Degrees, Minutes, Seconds);
	disp(LogMsg);
end

I got these results (C++ in Orbiter on the left, Matlab on the right:

C++					|	Matlab
-12.566 = -720d 0m 0.07s		|	-12.566 = -720d 0m 0.00s
-12.318 = -705d 44m 2673.27s		|	-12.318 = -705d 44m 33.27s
-12.069 = -691d 29m 1746.46s		|	-12.069 = -691d 29m 6.53s
-11.820 = -677d 13m 819.85s		|	-11.820 = -677d 13m 39.80s
-11.571 = -662d 58m 3493.04s		|	-11.571 = -662d 58m 13.07s
-11.322 = -648d 42m 2566.43s		|	-11.322 = -648d 42m 46.34s
-11.073 = -634d 27m 1639.63s		|	-11.073 = -634d 27m 19.60s
-10.824 = -620d 11m 712.82s		|	-10.824 = -620d 11m 52.87s
-10.576 = -605d 56m 3386.21s		|	-10.576 = -605d 56m 26.14s
-10.327 = -591d 40m 2459.41s		|	-10.327 = -591d 40m 59.41s
-10.078 = -577d 25m 1532.60s		|	-10.078 = -577d 25m 32.67s
-9.829 = -563d 10m 605.99s		|	-9.829 = -563d 10m 5.94s
-9.580 = -548d 54m 3279.19s		|	-9.580 = -548d 54m 39.21s
-9.331 = -534d 39m 2352.38s		|	-9.331 = -534d 39m 12.48s
-9.083 = -520d 23m 1425.77s		|	-9.083 = -520d 23m 45.74s
-8.834 = -506d 8m 498.96s		|	-8.834 = -506d 8m 19.01s
-8.585 = -491d 52m 3172.35s		|	-8.585 = -491d 52m 52.28s
-8.336 = -477d 37m 2245.55s		|	-8.336 = -477d 37m 25.54s
-8.087 = -463d 21m 1318.74s		|	-8.087 = -463d 21m 58.81s
-7.838 = -449d 6m 392.03s		|	-7.838 = -449d 6m 32.08s
-7.590 = -434d 51m 3065.33s		|	-7.590 = -434d 51m 5.35s
-7.341 = -420d 35m 2138.62s		|	-7.341 = -420d 35m 38.61s
-7.092 = -406d 20m 1211.91s		|	-7.092 = -406d 20m 11.88s
-6.843 = -392d 4m 285.10s		|	-6.843 = -392d 4m 45.15s
-6.594 = -377d 49m 2958.40s		|	-6.594 = -377d 49m 18.42s
-6.345 = -363d 33m 2031.69s		|	-6.345 = -363d 33m 51.68s
-6.097 = -349d 18m 1104.98s		|	-6.097 = -349d 18m 24.95s
-5.848 = -335d 2m 178.17s		|	-5.848 = -335d 2m 58.22s
-5.599 = -320d 47m 2851.47s		|	-5.599 = -320d 47m 31.49s
-5.350 = -306d 32m 1924.76s		|	-5.350 = -306d 32m 4.75s
-5.101 = -292d 16m 998.05s		|	-5.101 = -292d 16m 38.02s
-4.852 = -278d 1m 71.24s		|	-4.852 = -278d 1m 11.29s
-4.604 = -263d 45m 2744.54s		|	-4.604 = -263d 45m 44.55s
-4.355 = -249d 30m 1817.83s		|	-4.355 = -249d 30m 17.82s
-4.106 = -235d 14m 891.12s		|	-4.106 = -235d 14m 51.09s
-3.857 = -220d 59m 3564.36s		|	-3.857 = -220d 59m 24.36s
-3.608 = -206d 43m 2637.61s		|	-3.608 = -206d 43m 57.62s
-3.359 = -192d 28m 1710.90s		|	-3.359 = -192d 28m 30.89s
-3.110 = -178d 13m 784.14s		|	-3.110 = -178d 13m 4.16s
-2.862 = -163d 57m 3457.43s		|	-2.862 = -163d 57m 37.43s
-2.613 = -149d 42m 2530.68s		|	-2.613 = -149d 42m 10.69s
-2.364 = -135d 26m 1603.97s		|	-2.364 = -135d 26m 43.96s
-2.115 = -121d 11m 677.21s		|	-2.115 = -121d 11m 17.23s
-1.866 = -106d 55m 3350.51s		|	-1.866 = -106d 55m 50.50s
-1.617 = -92d 40m 2423.77s		|	-1.617 = -92d 40m 23.76s
-1.369 = -78d 24m 1497.04s		|	-1.369 = -78d 24m 57.03s
 -1.120 = -64d 9m 570.31s		|	-1.120 = -64d 9m 30.30s
 -0.871 = -49d 54m 3243.56s		|	-0.871 = -49d 54m 3.56s
-0.622 = -35d 38m 2316.83s		|	-0.622 = -35d 38m 36.83s
-0.373 = -21d 23m 1390.10s		|	-0.373 = -21d 23m 10.10s
 -0.124 = -7d 7m 463.37s		|	-0.124 = -7d 7m 43.37s
 0.124 = 7d 7m 463.37s			|	0.124 = 7d 7m 43.37s
 0.373 = 21d 23m 1390.10s		|	0.373 = 21d 23m 10.10s
 0.622 = 35d 38m 2316.83s		|	0.622 = 35d 38m 36.83s
 0.871 = 49d 54m 3243.56s		|	0.871 = 49d 54m 3.56s
 1.120 = 64d 9m 570.31s			|	1.120 = 64d 9m 30.30s
 1.369 = 78d 24m 1497.04s		|	1.369 = 78d 24m 57.03s
 1.617 = 92d 40m 2423.77s		|	1.617 = 92d 40m 23.76s
 1.866 = 106d 55m 3350.51s		|	1.866 = 106d 55m 50.50s
 2.115 = 121d 11m 677.21s		|	2.115 = 121d 11m 17.23s
 2.364 = 135d 26m 1603.97s		|	2.364 = 135d 26m 43.96s
 2.613 = 149d 42m 2530.68s		|	2.613 = 149d 42m 10.69s
 2.862 = 163d 57m 3457.43s		|	2.862 = 163d 57m 37.43s
 3.110 = 178d 13m 784.14s		|	3.110 = 178d 13m 4.16s
 3.359 = 192d 28m 1710.90s		|	3.359 = 192d 28m 30.89s
 3.608 = 206d 43m 2637.61s		|	3.608 = 206d 43m 57.62s
 3.857 = 220d 59m 3564.36s		|	3.857 = 220d 59m 24.36s
 4.106 = 235d 14m 891.12s		|	4.106 = 235d 14m 51.09s
 4.355 = 249d 30m 1817.83s		|	4.355 = 249d 30m 17.82s
 4.604 = 263d 45m 2744.54s		|	4.604 = 263d 45m 44.55s
 4.852 = 278d 1m 71.24s			|	4.852 = 278d 1m 11.29s
 5.101 = 292d 16m 998.05s		|	5.101 = 292d 16m 38.02s
 5.350 = 306d 32m 1924.76s		|	5.350 = 306d 32m 4.75s
 5.599 = 320d 47m 2851.47s		|	5.599 = 320d 47m 31.49s
 5.848 = 335d 2m 178.17s		|	5.848 = 335d 2m 58.22s
 6.097 = 349d 18m 1104.98s		|	6.097 = 349d 18m 24.95s
 6.345 = 363d 33m 2031.69s		|	6.345 = 363d 33m 51.68s
 6.594 = 377d 49m 2958.40s		|	6.594 = 377d 49m 18.42s
 6.843 = 392d 4m 285.10s		|	6.843 = 392d 4m 45.15s
 7.092 = 406d 20m 1211.91s		|	7.092 = 406d 20m 11.88s
 7.341 = 420d 35m 2138.62s		|	7.341 = 420d 35m 38.61s
 7.590 = 434d 51m 3065.33s		|	7.590 = 434d 51m 5.35s
 7.838 = 449d 6m 392.03s		|	7.838 = 449d 6m 32.08s
 8.087 = 463d 21m 1318.74s		|	8.087 = 463d 21m 58.81s
 8.336 = 477d 37m 2245.55s		|	8.336 = 477d 37m 25.54s
 8.585 = 491d 52m 3172.35s		|	8.585 = 491d 52m 52.28s
 8.834 = 506d 8m 498.96s		|	8.834 = 506d 8m 19.01s
 9.083 = 520d 23m 1425.77s		|	9.083 = 520d 23m 45.74s
 9.331 = 534d 39m 2352.38s		|	9.331 = 534d 39m 12.48s
 9.580 = 548d 54m 3279.19s		|	9.580 = 548d 54m 39.21s
 9.829 = 563d 10m 605.99s		|	9.829 = 563d 10m 5.94s
 10.078 = 577d 25m 1532.60s		|	10.078 = 577d 25m 32.67s
 10.327 = 591d 40m 2459.41s		|	10.327 = 591d 40m 59.41s
 10.576 = 605d 56m 3386.21s		|	10.576 = 605d 56m 26.14s
 10.824 = 620d 11m 712.82s		|	10.824 = 620d 11m 52.87s
 11.073 = 634d 27m 1639.63s		|	11.073 = 634d 27m 19.60s
 11.322 = 648d 42m 2566.43s		|	11.322 = 648d 42m 46.34s
 11.571 = 662d 58m 3493.04s		|	11.571 = 662d 58m 13.07s
 11.820 = 677d 13m 819.85s		|	11.820 = 677d 13m 39.80s
 12.069 = 691d 29m 1746.46s		|	12.069 = 691d 29m 6.53s
 12.318 = 705d 44m 2673.27s		|	12.318 = 705d 44m 33.27s

It's basically the same code, but the result is different. Why?




Also, bonus question:
How do I get the right ascension? I know that in an equatorial coordinate system where the x-axis points at RA = 0 and the y-axis points to celestial north, the right ascension is given by [math]RA = atan2 (Vector.y, Vector.x)[/math], where [math]Vector[/math] is obtained from

oapiGetRelativePos(TargetObject, ReferenceObject, &Vector);

But then the vector is oriented with x-axis out of 0N 0E (not far away from the west coast of central Africa). So what angle do I have to multiply with to orient x-axis towards RA = 0?

 

[PolliMatrix]

I think the problem is that "Minutes / 60" is evaluated as an integer - and rounded down to zero - because both "Minutes" and "60" are integers and this part is calculated before the first calculation that requires converting the number to a float (the subtraction).

The solution would be to use "60.0" or "60." instead of "60", so it gets cast to a float at the division.

 

[jedidia]

Bingo. I thought it must have something to do with the integers the moment I saw it, but didn't have the time to look closely at it.

 

[jangofett287]

This is c++, so you want "60.0f" Also float angle = (blah) needs the float treatment as well.

 

[martins]

Using "Minutes/60.0f" should work, but as an alternative, in your equation for seconds, couldn't you just pull a factor 60 into the parentheses to avoid the problem?

float Seconds = ((abs(Angle) * DEG - Degrees)*60 - Minutes) * 60;

Note that using explicit floating point constants ("60.0f" instead of "60") and explicit conversion to float ("(float)Minutes" instead of "Minutes") where appropriate is cleaner, but is really only required where the other operand isn't already a float. Otherwise an implicit conversion to float is performed.

Incidently, the #defined constant "DEG" expands to a double, so the entire right-hand side will expand to double, before being downgraded again to float by the assignment operator. This may prompt a compiler warning. You can avoid that by defining Seconds as double, or explicit conversion of the right-hand side:

float Seconds = (float)(((abs(Angle) * DEG - Degrees)*60 - Minutes) * 60);

 

[BrianJ]

Also, bonus question:
How do I get the right ascension? I know that in an equatorial coordinate system where the x-axis points at RA = 0 and the y-axis points to celestial north, the right ascension is given by [math]RA = atan2 (Vector.y, Vector.x)[/math], where [math]Vector[/math] is obtained from

oapiGetRelativePos(TargetObject, ReferenceObject, &Vector);

But then the vector is oriented with x-axis out of 0N 0E (not far away from the west coast of central Africa). So what angle do I have to multiply with to orient x-axis towards RA = 0?

Once you have got your reference vector using

oapiGetRelativePos(TargetObject, ReferenceObject, &Vector);

you need to rotate the vector by Earth Obliquity (~23.4 deg) around the X-axis, then you can use
[math]RA = atan2 (Vector.y, Vector.x)[/math]to find the RA.
Any help?

 

[Col_Klonk]

Kerbal Space Program suffers from a problem which I'm sure is Floating Point Runaway (that Kracken oscillation), which is probably similar to this, or due to them not using significant digits with their Float (Real4 or REAL8) numbers.

You'll probably find that Matlab uses a different FP rounding mode (this is different to Integer Floor or Ceiling) to C/C++. After a lot of calculations in the equation you get resolution loss, or disparate results. As mentioned use floats for everything except the last result... It's best to keep the results unassigned to any variable (ie stored in the FPU/SSE unit). I'm not sure whether you can do this in C++ or C, as soon as you assign it to a variable, you reduce the results from 80 bit (or 64) to 32 bit (Floats). I haven't tried this but I think you can do this by making the the calculation one long single line in C++, instead of 3.. if that is possible?

Wrt to significant digits... decide on the resolution you desire, and beyond that run an average of random numbers beyond the least significant digit. Add this average to the result before Floor'ing or Ceiling it.
I'm not sure whether this can be done in C or C++, but you can also control the rounding mode of the FPU or SSE calculations... Flipping between two rounding modes halfway through calculations, should help 'eliminate' that tiny creeping error.
A slight performance loss.. but at least accuracy would be more consistent.
With calculations done on the GPU.. I think the situation is the same.. I'll check on this