Saturn/Neptune positions

[george7378]

Hi everyone,

I have been looking through the algorithm on this page:

http://stjarnhimlen.se/comp/ppcomp.html

http://stjarnhimlen.se/comp/tutorial.html

...and I decided to make a program to calculate the positions of the planets. I've finished it (it's available here: http://gkastro.0sites.org/Ephemcalc.exe) and it seems to give pretty accurate positions for most of the planets when compared to calsky. However, I've noticed that the positions for Saturn and Neptune seem to be pretty far off what Calsky says, so I was wondering if it anyone could explain why? Is there a general inaccuracy when calculating Neptune and Saturn's positions, or is it likely ot be just me?

I haven't included the peturbations in the program (but there aren't any described for Neptune anyway).

Also, the program sometimes outputs negative right ascension values, which doesn't seem to be possible. Can anyone tell me what to do if I get negative R.A. values (do I need to add/subtract a set value until RA falls within a certain boundary)?

Thanks!

 

[RisingFury]

Ra is defined from 0 to 24 hours, so if you get a negative value, add 24 hours to it.

This should do the trick:

Ra = Ra + 24;
Ra = Ra - int(Ra / 24) * 24;

As for the perturbations... I'm not sure if Saturn suffers from any, but I know that Neptune does. Before the discovery of Pluto, it was believed that an Earth mass planet is causing them. IIRC they're caused by the Kupier belt, but I'm not sure...

---------- Post added at 08:17 ---------- Previous post was at 08:12 ----------

The website describes perturbations for Jupiter, Saturn and Uranus... try including those to see if the position for Saturn get corrected.

 

[george7378]

OK - I've corrected the RA problem and it works fine now (thanks for the tip!). Now that I have corrected RA, I can see that the error for Saturn and Neptune only occurs in the Declination output - RA is actually very accurate (+/- 30 seconds). I'll do the peturbations and see if it makes a difference for Saturn, and then I'll try and find something for Neptune.

EDIT: Maybe it's my iteration technique for finding the eccentric anomaly. At the moment, I'm just iterating 20 times (no matter what the difference between the results), like so:

             for (int it=0; it<=1; it++){

	          it = 0;

iterstart:

	          it++;

	          double E = E0 - ( E0 - e*(180/PI) * sin(DEGRAD*(E0)) - M ) / ( 1 - e * cos(DEGRAD*(E0)) );
	          E0 = E;

	          if (it < 20)
		      goto iterstart;

	          else goto cont;

	          }

cont:

//Rest of code.

I'd expect it to converge after 20 iterations no matter what the inputs, but still - I'm no expert.

 

[RisingFury]

OK - I've corrected the RA problem and it works fine now (thanks for the tip!). Now that I have corrected RA, I can see that the error for Saturn and Neptune only occurs in the Declination output - RA is actually very accurate (+/- 30 seconds). I'll do the peturbations and see if it makes a difference for Saturn, and then I'll try and find something for Neptune.

EDIT: Maybe it's my iteration technique for finding the eccentric anomaly. At the moment, I'm just iterating 20 times (no matter what the difference between the results), like so:

             for (int it=0; it<=1; it++){

              it = 0;

iterstart:

              it++;

              double E = E0 - ( E0 - e*(180/PI) * sin(DEGRAD*(E0)) - M ) / ( 1 - e * cos(DEGRAD*(E0)) );
              E0 = E;

              if (it < 20)
              goto iterstart;

              else goto cont;

              }

cont:

//Rest of code.

I'd expect it to converge after 20 iterations no matter what the inputs, but still - I'm no expert.

Ouch... Given the syntax I'm assuming you're working in C / C++. If so, you should drop the iterstart and cont tags. You can write the loop like this:

for (int I = 0; I < 20; I++)
{
   // Code goes here...
}

This will execute the loop 20 times, with the value for I being 0 on the first loop and 19 on the last. If you wanted 1 to 20 you could write like this:

for (int I = 1; I <= 20; I++)
{
   // Code goes here...
}

Your loop can then be written like this:

for (int I = 1; I <= 20; I++)
{
   double E = E0 - (E0 - e*(180 / PI) * sin(DEGRAD*(E0)) - M) / (1 - e * cos(DEGRAD*(E0));
   E0 = E;
}

You can interpret this loop as such: Start with I = 1, continue the loop until I equals 20, increase the value of I after each iteration has finished.

It's generally bad programming practice to use goto statements in your code... such code is difficult to read.

If you want to continue the loop until E has converged enough, you could try something like this:

double Margin = 0; // Define a reasonable value here, value must be greater then 0.

while (true)
{
   double E = E0 - (E0 - e*(180 / PI) * sin(DEGRAD*(E0)) - M) / (1 - e * cos(DEGRAD*(E0));

   if (abs(E - E0) < Margin)
   {
      E0 = E;

      break;
   }

   E0 = E;
}

In this case the loop is set to true, it will continue until it's broken from the inside. The condition that breaks the loop is the if statement. You can interpret it as such: If the difference between the previous value and the current value is low enough, break the loop. You can enter a reasonable value for the margin yourself. The abs() stands for absolute, so you get rid of any potential negative sign. You're only interested in the difference, not it's sign...

 

[george7378]

Thanks for your help - I've added the last bit of code to do the iteration until it converges (which is better practice than just doing it 20 times and hoping that it converges, I suppose).

I also found that the problem with the declination values is to do with the way I'm converting the decimal declinations to degrees, minutes and seconds. I think it's fixed now, and it's giving pretty nice data:

http://gkastro.0sites.org/Ephemcalc.exe

Thanks again!