skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

Calculation of whether the moon can be seen with the eye at the beginning of the lunar month #889

Closed sgbmzm closed 10 months ago

sgbmzm commented 11 months ago

Previously I asked here if there is a possibility of such a feature in skyfield. Now I publish the function I made.

`

#This function returns information on whether it is possible to see the moon with the eye on a certain date and time
#This is useful for the first moon in a lunar month or the last moon in a lunar month
#This function makes use of the research of Dr. Roy Hoffman

def Hoffmann_moon_normal_vision_factor_at(any_time_datetime):
    any_time = load.timescale().from_datetime(any_time_datetime)
    sun_topocentric = (eph['earth'] + location).at(any_time).observe(eph['sun']).apparent()
    moon_topocentric = (eph['earth'] + location).at(any_time).observe(eph['moon']).apparent()
    sun_alt, sun_az, sun_dist = sun_topocentric.altaz()
    moon_alt, moon_az, moon_dist = moon_topocentric.altaz()
    moon_percent = 100.0 * moon_topocentric.fraction_illuminated(eph['sun'])            
    moon_equatorial_radius_km = 1738.1 
    moon_apparent_diameter = Angle(radians=np.arcsin(moon_equatorial_radius_km / moon_dist.km) * 2.0)
    moon_thickness = moon_apparent_diameter.arcminutes() * (moon_percent / 100)
    # Hoffman's vision criterion. The formula is: the height difference between the sun and the moon, plus 7.37, times the square root of the thickness of the crescent 
    moon_vision_factor = abs(moon_alt.degrees-sun_alt.degrees) + (7.37 * (moon_thickness ** (0.5))) 
    # Normalization of the Hoffman vision coefficient so that zero means no vision in the eye at all and one means that there is certain vision in the eye and between zero and one is the range of uncertainty
    moon_normal_vision_factor = (moon_vision_factor - 12.2) / 5.2 

    if moon_normal_vision_factor <= -0.7:
        print(f"moon_normal_vision_factor is: {moon_normal_vision_factor:2.01f}, The moon does not! will be seen, even with binoculars")
    elif moon_normal_vision_factor > -0.7 and moon_normal_vision_factor < 0:
        print(f"moon_normal_vision_factor is: {moon_normal_vision_factor:2.01f}, The moon may be visible through binoculars")
    elif moon_normal_vision_factor >= 0 and moon_normal_vision_factor < 0.5:
        print(f"moon_normal_vision_factor is: {moon_normal_vision_factor:2.01f}, The moon may be visible to the eye, with great difficulty")
    elif moon_normal_vision_factor >= 0.5 and moon_normal_vision_factor < 1:
        print(f"moon_normal_vision_factor is: {moon_normal_vision_factor:2.01f}, The moon will probably be visible, there may be difficulty")
    elif moon_normal_vision_factor >= 1:
        print(f"moon_normal_vision_factor is: {moon_normal_vision_factor:2.01f}, The moon will be easily visible to the eye, when there are no clouds")

    return moon_normal_vision_factor #or return moon_normal_vision_factor > 0

from skyfield.api import N, S, E, W, wgs84, load, Angle
from pytz import timezone
import datetime as dt
import numpy as np

location = wgs84.latlon(31 * N, 35 * E)

eph = load('de440.bsp')

israel = timezone('israel')
time = dt.datetime(2023, 8, 18, 14, 53, 59, 731882)
local_time = israel.localize(time)
Hoffmann_moon_normal_vision_factor_at(local_time)

`

brandon-rhodes commented 10 months ago

Thank you for sharing this code, hopefully it will help anyone else who is interested in making the same calculation! I'm going to close the issue since there is no bug in Skyfield to fix, but, happily, that won't prevent this issue from coming up in searches for folks who Google for how to perform this operation with Skyfield.