Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Apparent planetary magnitudes #210

Closed
ghost opened this issue Oct 11, 2018 · 16 comments
Closed

Apparent planetary magnitudes #210

ghost opened this issue Oct 11, 2018 · 16 comments
Assignees

Comments

@ghost
Copy link

ghost commented Oct 11, 2018

Unless I missed it completely, skyfield doesn't appear to compute planetary magnitude.

@brandon-rhodes
Copy link
Member

I'll check the Explanatory Supplement when I get the chance and see what it recommends.

@brandon-rhodes brandon-rhodes self-assigned this Nov 11, 2018
@brandon-rhodes
Copy link
Member

The Supplement does have some tables and formulae, so I'll keep this open until I have time to tackle magnitudes!

@brandon-rhodes
Copy link
Member

Oh, this appears to be dated even more recently, and comes with code — well, kind of; they say it's on SourceForge!?

https://arxiv.org/pdf/1808.01973.pdf

@brandon-rhodes
Copy link
Member

The 2,035 lines of Fortran seem to mostly consist of I/O to read in its own test cases. I'll see if I can extract the actual formulae.

@ghost
Copy link
Author

ghost commented Nov 30, 2018

Brandon,

My apologies. I actually did some research on this myself to answer https://astronomy.stackexchange.com/questions/27918/what-is-the-planet-relative-to-earth-that-shows-the-greatest-change-in-apparent/27954#27954, and should have linked to there.

Basically, illumination (magnitude is essentially the logarithm of illumination) varies exactly as you'd expect: as the square of distance from the Sun, as the square of the distance from the Earth, and linearly as the albedo of the planet.

The only real bugaboo is the phase and something called the "phase curve" (https://en.wikipedia.org/wiki/Phase_curve_(astronomy)). Planets don't reflect like mirrors, but they don't reflect like perfect Lambertian surfaces (https://en.wikipedia.org/wiki/Lambertian_reflectance) either.

Instead, the phase curve (illumination based on phase) is computed from the magnitude itself, at least for the major planets. I believe there are generic formulas for asteroids, comets, etc, but, for the major planets, it's basically just curve fitting.

In the FORTRAN code you're looking at, I believe this is ph_ang_factor or something.

@Bernmeister
Copy link

@brandon-rhodes
Copy link
Member

@Bernmeister — Thanks for pointing out that Stack Exchange question! I have just jumped in to provide an answer based on actual planetary observations, from this paper and source code:

https://arxiv.org/pdf/1808.01973.pdf

https://sourceforge.net/projects/planetary-magnitudes/

These are cited as well at the top of the planetary magnitude module that I started working on last year but have not yet had the time to finish:

https://github.com/skyfielders/python-skyfield/blob/master/skyfield/magnitudelib.py

I should really have added a link to that file to this issue back when I started working on it, but I was probably thinking "I'm almost done! I'll wait and tell them when I'm finished tomorrow!" and so I failed to update this issue with info about the approach or the file where work is ongoing.

Hopefully I'll soon find time to return to it; we’ll see! In the meantime, feel free to try out any of the routines that are already written, and let me know if they work for you.

@brandon-rhodes
Copy link
Member

I have just dug back in to the module to see where things stand. First, I went ahead and added some code to take a planetary position, generate the parameters needed by the magnitude routines, and return the routine’s result. So you can now do this:

from skyfield import api
from skyfield.magnitudelib import planetary_magnitude

ts = api.load.timescale()
t = ts.utc(2020, 7, 27)
eph = api.load('de421.bsp')
p = eph['earth'].at(t).observe(eph['mercury'])
m = planetary_magnitude(p)
print(m)

But several steps will be necessary to complete the work and get it ready to add to the documentation.

  1. Some of the routines work when t is a vector (like Mercury) but some don’t (like Jupiter).
  2. Much more difficult will be planets like Uranus and Mars, that need to know not only their own position and those of the Earth and Sun, but what part of their surface each of those bodies is facing.
  3. And to know where their surface is facing, Skyfield’s support for planetary reference frames will need to be further expanded. Discussion on that happens to have taken place just earlier today! See: Central meridian longitude of planets and Jovian satellites? #417 (comment)

But given a busy schedule over the next few weeks, I’m not sure how soon I’ll be able to sit down and tackle this — and I'm unsure what the API will look like that lets users provide the extra orientation data that the models need.

Maybe I should just go with simpler and less accurate models that don't need that extra data?

But then folks will probably worry that it doesn't match HORIZONS. Hmm.

@brandon-rhodes
Copy link
Member

Documentation for the prototype, which was just released in Skyfield 1.26, is here:

https://rhodesmill.org/skyfield/api.html#skyfield.magnitudelib.planetary_magnitude

@reza-ghazi
Copy link

Hello Brandon,

Speaking of apparent magnitude, as you might know, from 2021 (next year), Astronomical Almanac will use the formulas and algorithm of (Mallama, Hilton 2018) the document 1808.01973.pdf as you mentioned above. So it would be amazing if you focus on these formulas, and please do not think about low accuracy :-).

Anyhow, what about Lunar apparent magnitude?
I didn't even find any empirical formula for that on the web. Do you have any idea?

Thank you,

@brandon-rhodes
Copy link
Member

So it would be amazing if you focus on these formulas, and please do not think about low accuracy

I agree that we should focus on the modern formulas, yes! If you the any interest, I would welcome a pull request that adds Mars, Saturn, or Neptune to magnitudelib.py alongside the planets already there. Here, for our future reference (and for anyone else that might jump in and want to take a look at this issue), is the FORTRAN file associated with the above PDF:

Ap_Mag_V3.f90.txt

Have we searched, by the way, for whether anyone else has already translated this to Python? The FORTRAN is actually very easy to read compared to some FORTRAN code that I've seen before, but it still imposes the cost of someone transcribing it into Python.

Anyhow, what about Lunar apparent magnitude?

I don’t remember seeing that tabulated in an almanac, and I don't see any mention of it having just flipped through the Explanatory Supplement. I do remember seeing the Moon's magnitude in one of the diagrams in Guy Ottewell’s annual astronomy calendar, but I don't know what formula he was using.

@reza-ghazi
Copy link

I added to my schedule to see what I can do with the Fortran file.

Please see the following image that I generated from the Horizon web interface. There is an APmag column there, which is the object's apparent magnitude (here the Moon).

moon_ephemeris_20_11_01_sample

Thank you,

@brandon-rhodes
Copy link
Member

Does HORIZONS have a bibliography, I wonder, that might suggest the name of the Moon magnitude model they're using?

@reza-ghazi
Copy link

Hello,

Not actually, and I am still working to find about it from other resources.

Best,

@tmoliter
Copy link

tmoliter commented Jan 1, 2021

This would be a super super cool feature to have!!

@brandon-rhodes
Copy link
Member

With the recent release of Skyfield 1.40, magnitudes are now supported for all of the major planets. Feel free to comment further here if you run into any snags, but otherwise I think we can declare this issue complete!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

4 participants