Sun or moon, finding the date and time of crossing a specific longitude

Question

the Sun is currently in Virgo ♍, 13°7'9" or longitude 56.0762 degrees how does one calculate when exactly the Sun will reach this position again. Since the moon is doing a spin in 29.x days lets use that as example also, "she" is in Gemini ♊, 1°33'23" or 56.1074 degrees , can we calculate when will it reach Gemini 05° 45' 00" or 56.1263 degrees ( a random coordination) Any math or programming language will be just fine even tho I messed with swisseph ( solcross, mooncross) and skyfield and no result ...

Answer 1

In skyfield, to find the times when something happens, you build a function that takes a Time object, outputs an integer value, and has a step_days attribute, and pass that function to the find_discrete method to find the times when your function transitions from one integer or boolean value to the next during its search period.

In the below script, I create a function ecliptic_longitude_exceeds that returns functions that can be used in the find_discrete method. It's a simple function, that should work for the Sun or the Moon, because we can trust their ecliptic longitudes to be constantly increasing from 0° to 360°. Not as easy to work with planets (they'll backtrack, requiring more code to detect and ignore the 360° to 0° transition) We'll use the year 2023 as our time window.

from typing import Callable

import pytz
import skyfield.searchlib
from skyfield import api, timelib
from skyfield.jpllib import SpiceKernel, ChebyshevPosition
from skyfield.vectorlib import VectorSum
from skyfield.framelib import ecliptic_frame


def ecliptic_longitude_exceeds(
        longitude: float,
        target: str|VectorSum|ChebyshevPosition,
        ephemeris: SpiceKernel) -> Callable[[timelib.Time], int]:
    """
    Creates functions that check whether target ecliptic longitude exceeds
    value at a specified time
    :param longitude: Ecliptic Longitude in decimal degrees
    :param target: Object to be checked
    :param ephemeris: Ephemeris to be used.
    :return: find_discrete-compatible function
    """

    earth = ephemeris['earth']
    target = ephemeris[target] if isinstance(target, str) else target

    def function(time: timelib.Time) -> bool:
        """
        :param time: Time of Observation
        :return: True if ecliptic longitude >= longitude False otherwise
        """
        observation = earth.at(time).observe(target).apparent()
        _, observed_longitude, _ = observation.frame_latlon(ecliptic_frame)
        return observed_longitude.degrees >= longitude

    function.step_days = 60

    return function


def main():

    ephemeris = api.load('de421.bsp')
    ts = api.load.timescale()
    utc = pytz.timezone("UTC")
    start, stop = ts.utc(2023), ts.utc(2024)

    sun_exceeds = ecliptic_longitude_exceeds(
        longitude=56.1074, target="sun", ephemeris=ephemeris  
    )

    times, states = skyfield.searchlib.find_discrete(start, stop, sun_exceeds)

    longitude_times = list(time for time, in_state
                           in zip(times.astimezone(utc), states)
                           if in_state)

    print('\n'.join(str(lt) for lt in longitude_times))


if __name__ == '__main__':
    main()

This script returns:

2023-05-17 06:11:22.543808+00:00

Quickly checking in my desktop copy of Stellarium puts the Sun's Ecliptic longitude at that time as 56°06'28.8", which matches up.

Answer 2

You had the right idea with solcross and mooncross from swisseph. The functions utc_to_jd and jdet_to_utc can help convert from a date+time to a Julian date and back.

This code:

import swisseph as swe
from datetime import datetime, timezone

def jd_from_dt(dt_ut):
    """Converts a datetime in UTC to a Julian date in TT."""
    (jd_tt, jd_ut) = swe.utc_to_jd(dt_ut.year, dt_ut.month, dt_ut.day,
            dt_ut.hour, dt_ut.minute, dt_ut.second)
    return jd_tt

def dt_from_jd(jd_tt):
    """Converts a Julian date in TT to a datetime in UTC."""
    (yr, mo, day, h, m, s) = swe.jdet_to_utc(jd_tt)
    dt_ut = datetime(yr, mo, day, h, m, round(s), tzinfo=timezone.utc)
    return dt_ut

# examples from question

post_jd = jd_from_dt(datetime(2023, 9, 5, 23, tzinfo=timezone.utc))

sun1_lon = 150 + 13 + 7/60 + 9/3600
sun1a_jd = swe.solcross(sun1_lon, post_jd - 30)
sun1b_jd = swe.solcross(sun1_lon, post_jd + 335)
print('Sun at lon={0:.4f}:\n  {1}\n  {2}'.format(
        sun1_lon, dt_from_jd(sun1a_jd), dt_from_jd(sun1b_jd)))

moon1_lon = 60 + 1 + 33/60 + 23/3600
moon1a_jd = swe.mooncross(moon1_lon, post_jd - 3)
moon1b_jd = swe.mooncross(moon1_lon, post_jd + 24)
print('Moon at lon={0:.4f}:\n  {1}\n  {2}'.format(
        moon1_lon, dt_from_jd(moon1a_jd), dt_from_jd(moon1b_jd)))

moon2_lon = 60 + 5 + 45/60
moon2a_jd = swe.mooncross(moon2_lon, post_jd - 3)
moon2b_jd = swe.mooncross(moon2_lon, post_jd + 24)
print('Moon at lon={0:.4f}:\n  {1}\n  {2}'.format(
        moon2_lon, dt_from_jd(moon2a_jd), dt_from_jd(moon2b_jd)))

# example from notovny's answer

jan1_jd = jd_from_dt(datetime(2023, 1, 1, 0, tzinfo=timezone.utc))

sun2_lon = 56.1074
sun2a_jd = swe.solcross(sun2_lon, jan1_jd)
sun2b_jd = swe.solcross(sun2_lon, jan1_jd + 365)
print('Sun at lon={0:.4f}:\n  {1}\n  {2}'.format(
        sun2_lon, dt_from_jd(sun2a_jd), dt_from_jd(sun2b_jd)))

produces this output:

Sun at lon=163.1192:
  2023-09-05 22:54:26+00:00
  2024-09-05 04:39:52+00:00
Moon at lon=61.5564:
  2023-09-05 22:58:23+00:00
  2023-10-03 07:50:07+00:00
Moon at lon=65.7500:
  2023-09-06 06:44:16+00:00
  2023-10-03 15:23:09+00:00
Sun at lon=56.1074:
  2023-05-17 06:11:23+00:00
  2024-05-16 11:57:53+00:00