计算新月的倾斜度

Question

可以用Skyfield计算月牙的倾斜度吗？还是直接？我知道来回的月牙石，我想正确地显示它。

https://www.calsky.com/cs.cgi/Moon/15?obs=4311082095747103&c=

# https://rhodesmill.org/skyfield/toc.html 
from skyfield import api

tmsc = api.load.timescale()
planets = api.load('de421.bsp')

import datetime
import pytz

toposloc = api.Topos('50.9409116 N', '6.9576131 E') # Köln local position on earth
toposabs = planets['earth'] + toposloc # absolute position incl earth 
DTUTC = datetime.datetime.now(datetime.timezone.utc)

# calc moon
astro = toposabs.at(tmsc.utc(DTUTC.year, DTUTC.month, DTUTC.day, DTUTC.hour, DTUTC.minute, DTUTC.second)).observe(planets['moon'])
app = astro.apparent()
alt, azi, distance = app.altaz()
altazmoon=[alt.degrees, azi.degrees] # save for later
_, lonmoon, _ = app.ecliptic_latlon()

Answer 1

好的，感谢 stackoverflows 的建议，我在这里找到了这个答案，这可能就是我正在寻找的答案

Calculating moon face rotation as a function of Earth coordinates

我目前正在测试这个。至少在今天，如果我从窗外看，它看起来非常准确。好样的！

[edit:] 我已经测试了一周并更新了下面的代码。现在一切似乎都对了。感谢 Stackoverflow

import math

def moontilt(altazsun, altazmoon): # mathematical angle where the sun light shines from onto the moon disc (measured from 3 o’clock anticlockwise)
    dLon = math.radians(altazsun[1]-altazmoon[1])
    radaltsun = math.radians(altazsun[0])
    radaltmoon = math.radians(altazmoon[0])
    cosaltsun = math.cos(radaltsun)
    y = math.sin(dLon) * cosaltsun
    x = math.cos(radaltmoon) * math.sin(radaltsun) - math.sin(radaltmoon) * cosaltsun * math.cos(dLon)
    brng = math.atan2(y, x)
    return 90-math.degrees(brng) 

# calc sun
astro = toposabs.at(tmsc.utc(DTUTC.year, DTUTC.month, DTUTC.day, DTUTC.hour, DTUTC.minute, DTUTC.second)).observe(planets['sun'])
app = astro.apparent()
alt, azi, distance = app.altaz()
altazsun=[alt.degrees, azi.degrees] # save for later
_, lonsun, _ = app.ecliptic_latlon()

# get moon phase
moonphase = ((lonmoon.degrees - lonsun.degrees) % 360.0) / 360 #  0: new moon, 0.5: full moon, 1: new moon

# this should be the tilt from the standard displaying symbols in degrees of the crescent I draw with PIL
tilt=-moontilt(altazsun, altazmoon) # PIL.ImageDraw.Draw.chord measures angles increasing clockwise. Thus minus
if moonphase>0.5: tilt+=180 # after full moon the light is expected to shine from the opposite side of the moon (which is already taken into account in the standard moon icons)