Saemundsson formula dropped value

[izzatzubir]

While working using altaz(), I went through your earthlib.py file. It seems at line 148, you dropped an extra decimal.

The correct formula (from Meuss's Astronomical Algorithm, p107) is P/1010 * 283/(temp +273). However, on line 148, you used 0.28 instead of 0.283.

Can you please correct this? Or perhaps you have other reasons for this?

Thanks.

[brandon-rhodes]

The formula that Skyfield uses is based on the official United States Naval Observatory code for computing refraction:

https://github.com/brandon-rhodes/python-novas/blob/5128fdd4deeceec9dbecf936b751d2ef2aee4fd1/Cdist/novas.c#L9195

As you can see, it uses 0.28 instead of 0.283.

Because Skyfield is designed to give the same results as the Naval Observatory, it will need to keep using the constant 0.28 there.

Does Meeus explain where he himself got the number from? If we could find what paper or source the number 0.28 comes from, and where Meeus got 0.283 comes from, I would be happy to update the documentation of the refraction routine to explain how that constant was determined and where it comes from — what does Meeus say?

[izzatzubir]

From page 106-107, Meuss did not mention where he get the 283 from (converted to 0.283). However, this can be inferred from his statement regarding the assumed temperature of 10 degrees Celcius which then be converted to 283 Kelvin.

Thanks.

[brandon-rhodes]

I just tried changing the constant to 0.283 to see the result, and the result is that 43 tests fail, which is around 7% of all Skyfield tests. As an example, here's the first failure to show up on my screen:

  skyfield/tests/test_against_novas.py line 2159 in
  test_jupiter_barycenter_topocentric_date0
    compare(alt.degrees, 49.42056980196601, 0.0005 * arcsecond)
  skyfield/tests/test_against_novas.py line 25 in
  compare(<SpiceKernel 'de405.bsp'>)
    assert max(abs(value - expected_value)) <= epsilon
AssertionError: 0.00015050063548471826
      is not <= 1.3888888888888888e-07

Alas, this shows that a case where Skyfield agreed with the Naval Observatory library to within a half-milli-arcsecond, adjusting the constant creates a difference between the two libraries of about 1000 times that tolerance.

There is one other data source that I sometimes test Skyfield against: the NASA Horizons system.

https://ssd.jpl.nasa.gov/horizons/app.html#/

It's very tedious to test against, because instead of being scriptable like a local library, the HORIZONS system is most often manipulated using a complicated multi-step web form. But I note that HORIZONS does seem to include a table setting named:

I don't have time to try setting up a HORIZONS circumstance where an object would be right at the horizon and thus experience significant refraction, but if you have a few minutes to step through the interface and create such a circumstance, I think that a demonstration that Skyfield can't come up with the same numbers as HORIZONS — if it turns out they use 0.283 — would be a very good argument in favor of at least setting up Skyfield where users could choose to have 0.283 used instead of 0.28. All I would need is the full text output from a HORIZONS run (include all the 'small print' down at the bottom, which is important because it gives definitions), showing a computed altitude with refraction turned 'on' that is different from the refraction Skyfield produces. If you decide to tackle it, let me know what you discover!

[izzatzubir]

Hello brandon. I have checked the values, and the value you inputted is Correct!

0.28 actually comes from 283/1010 (10 degree Celcius divided by 1010 millibar) ~ 0.28019802.

Therefore the value is correct. We should close this issue.

[brandon-rhodes]

Thanks so much for tracking down the true value! I'll go ahead and close, then.