Calculate the tilt of the moon crescent

Question

Is it possible to calculate the tilt of the moon crescent with Skyfield? Or directly? I know the crescent rocks back and forth and I would like to display it correctly.

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

Okay, thanks to stackoverflows recommendations I found this answer here which might be what I am looking for

Calculating moon face rotation as a function of Earth coordinates

I am currently testing this. At least for today it looks very accurate if I look out of the window. Way to go!

[edit:] I tested it for a week now and udpated the code below. Now everything seems to be on point. Thanks 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)