skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.39k stars 209 forks source link

Non-physical gravitational deflection corrections for Solar System bodies #734

Open andreww5au opened 2 years ago

andreww5au commented 2 years ago

When correcting for the gravitational deflection of light by Solar System bodies, deflections are calculated assuming each body is a point mass, up to the point where the position of the distant object and the centre of the Sun differ by 1 arcsecond. For most objects, this is OK, but for the Sun, this assumption produces huge (half-degree) deflections for few-arcsecond separations that would only exist in reality if the Sun had collapsed to a black hole.
It would be nice if the deflectors dictionary in relativity.py included the object diameter as well as the relative mass, so the cutoff could be changed to the apparent angular radius of the deflector, instead of a constant 1 arcsecond.

brandon-rhodes commented 2 years ago

Yes, that's a good point—and is an approach I have thought of trying for another similar artifact that Skyfield has inherited from the USNO’s NOVAS Fortran code on which all of Skyfield's astrometry is based. Here's a test that verifies the current behavior:

https://github.com/skyfielders/python-skyfield/blob/master/skyfield/tests/test_earth_deflection.py#L12

And a notebook in which I did the investigation behind that test:

https://github.com/skyfielders/astronomy-notebooks/blob/master/Skyfield-Notes/Fixing-earth-deflection.ipynb

A Skyfield user has suggested that the abrupt discontinuity at the point that an object is 0.8 of the way behind Earth's limb might be resolved by instead considering the amount of deflection to start decreasing once the object is behind the Earth, by the proportion of the Earth's mass that's still between the object and the Earth's center.

I wonder if we could do the same thing in the case you've encountered: once the object is behind the limb of the deflector, could we reduce the amount of deflection by the fraction of (in this case) the Sun's mass that still sitting on the "inside" of the path of the deflected light? That would eliminated the discontinuity, and the deflection would drop smoothly to zero for an object directly behind the Sun.

I don't have a good physical sense of how gravitational deflection really works for light that passes through the deflector's mass, though, so let me know if there's a better model that could be applied.

If we could come up with a solution that eliminates discontinuity, it's a big help to folks using solvers and iterator to catch particular circumstances. What do you think?

andreww5au commented 2 years ago

I don't think using the fraction of the Sun's mass inside the impact parameter would work even for a purely Newtonian calculation - for a stationary object inside the Sun, gravitational forces from all the mass in the Sun outside that point cancels out, but that's not true for the time a moving body spends approaching the Sun, then departing. For gravitational lensing, it's much more complicated - the formula you're using now is for a point mass, and doing it properly for rays passing through the Sun would mean knowing the Sun's density as a function of radius, collapsing it to a lensing plane perpendicular to the observer, and integrating that deflection function over the surface density in the two-dimensional lensing plane. I don't think there's anything wrong with a sharp discontinuity at the limb of the Sun, because there is a real, physical discontinuity there - the Sun isn't really transparent to anything, apart from neutrinos, and we can't measure accurate directions for them. Earth deflection might be worth that hassle, but I doubt it...

brandon-rhodes commented 1 year ago

I don't think there's anything wrong with a sharp discontinuity at the limb of the Sun

The problem is that solvers and position-finding routines that happen to cross that boundary won't know its there and will have a big problem with any discontinuity. I think we need an approach to deflection that provides a position curve without first-order discontinuity, if solvers are to work well that happen to reach those boundaries?