Formulas to determine the illuminated phase and orientation of the Moon

Question

The phase of the Moon (the bright limb) looks different depending on where on Earth the observer is located and the time of the observation. I've located this video and the accompanying spreadsheet that supposedly handles this sort of calculation, however I am somewhat confused on what formulas are actually used. I would greatly appreciate if someone with a bit more experience could decipher the spreadsheet a bit and perhaps write down or link me to the formulas being used.

I'm most interested in determining the illuminated phase (i.e the size of the bright limb) and the apparent orientation of the illuminated portion relative to the horizon. My goal is to understand how this type of calculation is done, in order to re-implement it in code for an application to use - displaying the typical crescent moon icon for night-time, but at the correct angle for a specific location and time.

PS! First post here, so I hope I have used the correct terminology for things.

EDIT: Eventhough the author of the spreadsheet left an answer to the topic, it doesn't feel like it really answers the question. Not to sound clueless and needy, but the question talks about formulas, which I know can be derived from the spreadsheet, however as I have no idea how to approach this, then the answer I was hoping for would help a bit more.

Answer 1

I created the spreadsheet. It's a bit complicated, to say the least. Initially I was inputting the formulas by hand from the book I cited in the video, but the process was rather error-prone so I ditched those tables and simply scanned the pages electronically and used the numbers from the scans to construct the tables of the Fourier series used to compute the moon's position. Some of the tables are therefore redundant and aren't used. You may be able to use a much simpler method to compute the position of the moon though; for my purposes I needed arcsecond accuracy. If all you're interested in is computing the approximate phase angle and apparent orientation you can get away with a much more simplified calculation for the moon's position that is accurate to within a few tenths of a degree. In that case the main thing you need from my spreadsheet are the formulas for computing the phase angle and apparent orientation. You'll need to calculate the approximate location of the sun as well in order to determine these things. The spreadsheet does this as well, using a simplified method described by Jean Meeus in "Astronomical Formulae for Calculators." He also gives a simplified method for computing the position of the moon. Either "Astronomical Formulae for Calculators" or "Astronomical Algorithms" (both by Jean Meeus) would be excellent references for your project. I can walk you through the formulas you need if you need help with that part.

Answer 2

I was struggling for a while to find a ready-made formula to calculate Moon orientation at my location. Eventually, after several bits and pieces from the internet, I now have a formula for it (for my Northern Hemisphere, not sure whether it works for Southern Hemisphere). I've written a more detailed blog describing 3 solutions you can get a Moon phase image file with the right moon orientation/tilt as seen from your location/sky here Moon Phase orientation as seen from your location/sky

The key value you need to calculate is the Moon Zenith Angle ((Moon Phase Angle (or called Position Angle of the Moon's Bright Limb in Jean Meeus's book) minus Moon Parallactic Angle)). You also need a have a correct/current Moon Phase with North-up orientation (vertical Moon Phase with Rabbit ears at 90 degree clockwise from North Starting Moon Position) to start with (I can get the current Moon Phase png file from CumulusMX software, it always produces the current Moon Phase with North-up orientation. I can also DIY the current Moon Phase, though a bit complicated. Let me know if you want to know how.). From there, you may follow the orientation logic from my python code below based on the above assumption. I've verified its orientation with timeanddate moon phase website, Stellarium, and Skysafari.Click this link for the current oriented Moon at my Lat 14d N, 100.5E: current oriented Moon Phase

#def moonorientation():
# compute moon zenith angle, alitude, and azimuth at current time
# compute moon orientation as seen from my location by rotating moon.png...
# ... which is a moon phase image created by CumulusMX, and then save the rotated moon phase...
# ... in /var/www/html/moonphase_CumulusMX_rotated.png
from skyfield.api import N, E, load, wgs84
from skyfield.trigonometry import position_angle_of
from skyfield.framelib import ecliptic_frame
import datetime
# get current time in utc
now = datetime.datetime.utcnow()
ts = load.timescale()
t = ts.utc(now.year, now.month, now.day, now.hour, now.minute, now.second)

eph = load('de421.bsp')
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']
# find Moon phase (0-180 = new to waxing to full moon, >180-0 is waning moon)
e = earth.at(t)
_, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
_, mlon, _ = e.observe(moon).apparent().frame_latlon(ecliptic_frame)
phase = (mlon.degrees - slon.degrees) % 360.0
print("moon phase (180deg=full)",phase)
# find moon zenith angle, need observer location
# set my location
myLat = 13.71732
myLng = 100.5907
europacafe = earth + wgs84.latlon(myLat * N, myLng * E)
b = europacafe.at(t)
m = b.observe(moon).apparent()
s = b.observe(sun).apparent()
#print(position_angle_of(m.altaz(), s.altaz()).degrees)
moonangle = position_angle_of(m.altaz(), s.altaz()).degrees
print("zenith moon angle: ", moonangle)
# compute moon phase bright limb degree to be oriented from what seen at north pole
if phase > 180: #waning moon
    moonorient = moonangle - 90
else: #waxing moon
    moonorient = moonangle + 90
# as rotate function in html is limited to +180 and -180, so need adjustment
if moonorient > 180:
   moonorient = moonorient - 360
elif moonorient < -180:
    moonorient = moonorient + 360
moonorient = moonorient * -1


print("moon phase rotation angle, from upright phase, at my location: ", moonorient)  
# from the above angle degree, rotate clockwise if it is positive degree
#return moonorient
#print(m.altaz())

# rotate moon phase (py rotate +deg = anti-clockwise)
from PIL import Image 
img = Image.open(r"/home/bthoven/CumulusMX/web/moonES.png")
#rotate_img= img.rotate(moonorient, resample=Image.BICUBIC, expand=True)
rotate_img= img.rotate(moonorient * -1, resample=Image.BICUBIC)
#rotate_img.show() 
rotate_img.save('/var/www/html/moonphase_CumulusMX_rotated.png') 

Answer 3

The Moon's or a planet's bright limb orientation can be obtained by calculating the Position Angle of the Moon's Bright Limb. For detailed explanations and formula, you can see the following source:

Jean Meeus, Astronomical Algorithms, Secon Edition, 1998, PP 345-347