skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 211 forks source link

The Pole and the circumpolar #663

Closed lucaberti closed 2 years ago

lucaberti commented 2 years ago

Hello, is there a function to know if a target is circumpolar at a given location? I tried to search in the API references and in the code, but I did not find it. I also tried to write a function to do it myself, but I ran in some problems. The assumtion is that an object becomes circumpolar at a certain Latitude when it has a greater declination than 90° - the alt of the pole on that location. Now, in order to compute the pole alt, I create a Star object like this: thepole = Star(ra_hours=(00, 00, 00), dec_degrees=(90, 00 , 00), epoch=2451545.0) and I compute its alt like this: thepole_alt, az, distance = (ephemeris['earth']+topos).at(time).observe(thepole).apparent().altaz() Checking the result against some other planetaries software, I get that the pole is a little off. For my case, I get 46° 06' 52.3" in Skyfield and 46° 12' 30ish" in the planetaries (the location the same in all cases).

Then I have this dumb problem: I'm not able to subtract two Angle objects. I tried to convert it in float numbers, but I have to admit that I did not catch how the Angle class works. Thank you for any help!

brandon-rhodes commented 2 years ago

How are you choosing time? A circumpolar star's altitude varies all night. As it circles the pole its altitude will be at a minimum when below the pole and at a maximum about 12 hours later when it is above it.

I think you want to compare the object's declination to the latitude of the topos. If the topos is at 30° N, for example, then any object with a declination of 90° - 30° = 60° or higher will be circumpolar, if I'm remembering the definition correctly?

lucaberti commented 2 years ago

Time is the time of observation. I do not calculate altitude of the target, but the altitute of the pole. By definition it should not change during the day.

Does the 90°-L account for pole account for precession and nutation?

brandon-rhodes commented 2 years ago

Does the 90°-L account for pole account for precession and nutation?

No, but I would expect the effects of precession and nutation to be small compared to the inherent error in the concept of “circumpolar” itself, since of course true circumpolar-ness depends on obstacles on the horizon and on the temperature and humidity raising objects at the horizon by some particular number of degrees. If you could describe your problem in a bit more detail, we could work out what kind of precision you need and expect from the calculation.

What do you mean by “in the planetaries”? Is that the name of software, or maybe of book with astronomy coordinates? And, what source code are you using in each case to compute the angle of the pole above the horizon? I would need to see the source code in both cases to determine what difference in the choice of coordinate system or Earth model might be making the small difference in altitude.

Note that determining the sky position of the pole is only half of the problem — since the pole wanders around to different RA/dec positions over the months and years, you will then need to compute the angular separate between the true pole and your object. Simply using the declination won't work, since that doesn't account for polar motion. (Or maybe you're ignoring that? It's difficult to guess which effects you want to include without seeing your code, and learning more about the problem you're trying to solve.)

I should have mentioned in my first answer, but hit "Comment" too soon: to subtract two angles in degrees, I would try a1.degrees - a2.degrees and Python should return a floating point number of degrees.

lucaberti commented 2 years ago

Sorry, for "in the planetaris" I ment softwares like Stellarium, C2A and Cartes du Ciel. But you are right: the pole changes RA/dec, so my idea is, at best, a silly propposal. Sorry. Nevertheless this is the code:

thepole = Star(ra_hours=(00, 00, 00), dec_degrees=(90, 00 , 00), epoch=2451545.0)
thepole_alt, az, distance = (ephemeris['earth']+topos).at(time).observe(thepole).apparent().altaz()

So, probably the best course of action would be to search for the nearest antimeridian pass and check if alt is above ore below ground. Or simply just apply the 90°-L rule ;)

brandon-rhodes commented 2 years ago
thepole = Star(ra_hours=(00, 00, 00), dec_degrees=(90, 00 , 00), epoch=2451545.0)
thepole_alt, az, distance = (ephemeris['earth']+topos).at(time).observe(thepole).apparent().altaz()

If I understand this code correctly, it computes the position of the pole back on 2000 January 1 and then asks, on a current date, how high that point is in the sky. I suspect that this is not what you want — the pole has now precessed past that position, and while it might be of historical interest to be able to point up and say "See right there? That's where the pole was way back in the year 200!", it will probably not do the work of helping to identify current-day circumpolar objects.

So, probably the best course of action would be to search for the nearest antimeridian pass and check if alt is above ore below ground.

That sounds like a good plan! Though I'll admit it will be a bit expensive. You could possibly cache the result: for, say, each longitude 0–360°, build a series of Star objects at different declinations until the expensive antimeridian test reveals whether it's circumpolar, and zero in on the exact declination where that answer changes. You will then have a lookup table for the crucial declination at which something ceases to (or begins to be) circumpolar, across the whole sky of 360 meridians.

Then you could check particular objects later against the declination for the nearest line of longitude.

Or simply just apply the 90°-L rule ;)

My guess is that the 90°-minus-latitude rule will give you an answer within the precision of other factors like atmospheric refraction and horizon obstacles — but as I don't know your application, that's only a guess. :)

brandon-rhodes commented 2 years ago

Since it sounds like you have a few good ways forward here, I'm going to go ahead and close this issue, but feel free to comment further if you still have questions.