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()

Inspired by your question, the next release of Skyfield will include a function that computes the position angle for you. Here's a preview of the documentation that shows how to use the function: https://github.com/skyfielders/python-skyfield/blob/2c457ecbc18f6965447b2c184f2ac125a052ba80/skyfield/documentation/examples.rst#at-what-angle-is-the-sun-to-the-crescent-moon — Brandon Rhodes, Dec 05 '19 at 23:37
Oh, you are the developer of Skyfield. How very nice to meet you. Also thank you for this new addition. I really really like Skyfield. — SirDagen, Dec 07 '19 at 10:38

SirDagen · Answer 1 · 2019-12-14T06:52:05.810

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)

If it helps simplify your code, note that `alt.radians` and `azi.radians` will give you radians directly and let you avoid the calls to `math.radians`. (But maybe you prefer values that are meaningful if you need to print them out during debugging?) — Brandon Rhodes, Dec 05 '19 at 13:57
Thank you very much for pointing this out, Brandon. Surely it makes more sense to work in rad. (My program is still in an experimental state and it is currently easier for me to work in degrees.) — SirDagen, Dec 05 '19 at 18:11

Calculate the tilt of the moon crescent

1 Answers1