Unable to detect planetary occultation by the moon

[meiphoo7-Mae]

As I already told in #980 I'm still trying to calculate several astronomical events with Skyfield.
Now I want to see if it is possible to find a planetary occultation by the moon. I've searched the internet and this question (some mention Skyfield) shows up several times, but apparently no one tried it. So I did.

I use this function (based upon what I found in your source code):

def separation_t1_t2(target1, target2):
   def separation(t):
      e = eph['earth'].at(t)
      tg1 = e.observe(target1).apparent()
      tg2 = e.observe(target2).apparent()
      return tg1.separation_from(tg2).degrees

   separation.step_days = 1/24
   return separation

And feed that to find_minima to figure out the smallest value returned by the function.

When I try this with the objects moon and saturn barycenter I get a minimum returned for August 21 2024 which is correct because at that date there was an occultation of Saturn (at least in the part of the world where I live). However the result is not entirely correct because the smallest separation found is larger than half the width of the moon, so that means no occultation. According to Stellarium there was an occultation at my location. Also the found time is a bit off because it is a little less than a hour to early.
I assume that the minimum distance between the objects occurs during the occultation but apparently I'm wrong.
Do you see a problem with this approach?
I hope this is not the same thinking error as in #980...

[brandon-rhodes]

I think you are right that an occultation is considered to be at its maximum when the angular separation is smallest. Could you share a short script that prints an answer, and share the date and time that you wish it was printing out instead? Thanks!

[meiphoo7-Mae]

I thought the snippet I provided was good enough, but I understand the problem. Here is a working script:

from skyfield.api import N, S, E, W, load, wgs84
from skyfield.searchlib import find_minima

#from imports import *
#from common_data import *

ts = load.timescale()
tnow =ts.now ()
eph = load('de421.bsp')

location = wgs84.latlon(52.092876 * N, 5.104480 * E)
observer = eph['Earth'] + location


t0 = ts.utc(2024, 8, 1)
t1 = ts.utc(2024, 8, 31)

moon = eph['moon']
saturn = eph['saturn barycenter']

def separation_t1_t2(target1, target2):
   def separation(t):
      e = eph['earth'].at(t)
      tg1 = e.observe(target1).apparent()
      tg2 = e.observe(target2).apparent()
      return tg1.separation_from(tg2).degrees

   separation.step_days = 1/24
   return separation

f = separation_t1_t2(moon, saturn)
t, y = find_minima(t0, t1, f)
for ti, yi in zip(t, y):
   text = f'{ti.utc_iso()} Smallest angular distance: {yi}'
print (text)

This gives:
2024-08-21T02:41:38Z Smallest angular distance: 0.40688657224346086

I don't know the exact time or distance that it is supposed to be because I compared it with in-the-sky.org and that site seems to be down at the moment I write this.

This source code is not optimized because I adapted it from a much larger project that I created.

[brandon-rhodes]

Thanks for the full script! It helps me see exactly what the arguments to each routine are.

It took me several readings to realize that you carefully build the observer but then instead use the geocenter for your angle:

e = eph['earth'].at(t)

Try doing:

e = observer.at(t)

Hopefully that will drop the angle low enough that Saturn goes behind the Moon!

[meiphoo7-Mae]

Wow, thank you.
That was a really quick response...

I found something interesting: The site is up again and the result of the original version of the script (geocenter) is exactly the same as: https://in-the-sky.org/news.php?id=20240821_15_100&latitude=52.092876&longitude=5.104480&timezone=%2B02%3A00.
Even the distance of about 24 arcminutes is right.
I've searched on internet and an appulse is indeed calculated for the location of the observer (apparent).
The adapted version with your correction seems to verify with Stellarium, but the exact time and minimum distance is difficult to detect with Stellarium.

But I gets weirder:
I've tried another possible occultation on July 24.
Now my original version produces: 2024-07-24T20:28:25Z Smallest angular distance: 0.34716988131943555, which verifies with https://in-the-sky.org/news.php?id=20240724_15_100&latitude=52.092876&longitude=5.104480&timezone=%2B02%3A00.
The adapted version with your correction produces: 2024-07-24T20:08:08Z Smallest angular distance: 0.6143410093807207, so that means no occultation which is confirmed by Stellarium.

It looks like in-the-sky.org uses an 'alternative' method to calculate the moment of appulse.
Anyway, my goal was to detect a planetary occultation by the moon and with your correction that is possible.
Once again: thank you!

[brandon-rhodes]

I'm glad we got your script working, and successfully detecting occultations! I hope you continue to enjoy using Skyfield.

[meiphoo7-Mae]

I certainly will. Skyfield is the best software for astronomical calculations that I've encountered so far: Only the mathematics are provided, no predefined layout stuff, so complete freedom of the presentation of the results. It's not very difficult to devise new calculations (if the theory behind those calculations is completely understood ;-) ).
So, yes, I hope to use Skyfield for many more calculations.