¿Cómo calcular el ángulo exacto a la luna?

Para un cálculo único, la solución más sencilla sería utilizar un software de planetario, como stellarium.

Para una solución de programación, un paquete de python como pyephem es una forma eficiente de calcular la posición de la luna en múltiples momentos, de una manera que se puede importar a una hoja de cálculo.

El movimiento aparente de la luna es más complejo de lo que piensas. La luna se mueve unos 15 grados por hora debido a la rotación de la Tierra. Pero además de este movimiento bastante simple, está la órbita real de la Luna alrededor de la Tierra, que es elíptica (por lo que este movimiento no es uniforme, la Luna se mueve más rápido cuando está cerca de la Tierra) y la posición e inclinación de la órbita es perturbado por el sol. Las perturbaciones son bastante regulares y se pueden calcular, pero estos efectos (la rotación de la tierra, la órbita de la luna, la excentricidad de la órbita, la perturbación de la órbita) se combinan para hacer que el cálculo real sea más complejo que "una fórmula" que puede caber fácilmente en una hoja de cálculo, por lo que se recomienda un paquete como pyephem.

En esta respuesta , menciono una fuente donde se puede encontrar una fórmula útil para el comportamiento aproximado de los objetos celestes.

Algoritmos astronómicos de Jean Meeus

No tengo una copia personalmente, así que no puedo verificar si hay exactamente lo que necesita allí o no, pero este pdf de un pequeño subconjunto de páginas de la cuarta edición incluye la tabla de contenido, que incluye:

30. Position of the moon....................... 147
31. Illuminated Fraction of the Moon’s Disk.... 155
32. Phases of the Moon......................... 159

así que hay una buena probabilidad de que en esas siete u ocho páginas del capítulo 30 haya justo lo que está buscando.

Creo que los "ángulos 3D" que menciona serían solo altitud (o ángulo de elevación) y acimut. El primero es el ángulo por encima (o por debajo) del horizonte, y el segundo es como la "dirección de la brújula", excepto que comienza en el norte verdadero, no en el norte magnético, por lo que es solo 2D.

Si se siente cómodo usando Python, uno de los lenguajes informáticos más fáciles, existe un paquete llamado Skyfield . Creo que debería poder ejecutarlo en una Raspberry Pi, por ejemplo, si su creación es electromecánica. El paquete es muy preciso, por lo que generalmente necesita una conexión a Internet para "verificar cosas", lo que a veces es un poco frustrante.

También puede calcular una tabla para su ubicación, luego guardarla y usarla para ejecutar una pequeña computadora o controlador, si está usando uno.

GREEN  = Moon (being made of green Cheese of course)
YELLOW = Sun

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos

halfpi, pi, twopi = [f*np.pi for f in 0.5, 1, 2]
degs, rads        = 180/pi, pi/180

load   = Loader('~/Documents/fishing/SkyData')
data   = load('de421.bsp')
ts     = load.timescale()

de421  = load('de421.bsp')
earth  = de421['earth']
sun    = de421['sun']
moon   = de421['moon']

Haifa   = earth + Topos(32.82, 34.98)
minutes = np.arange(60*24)
times   = ts.utc(2018, 9, 24, 0, minutes)

salt, saz, sd = Haifa.at(times).observe(sun).apparent().altaz()
malt, maz, md = Haifa.at(times).observe(moon).apparent().altaz()

saltr, sazr, maltr, mazr = [x.radians for x in (salt, saz, malt, maz)]

sr, mr = [1. - x/halfpi for x in (saltr, maltr)]

sx, sy = [sr * f(sazr) for f in (np.sin, np.cos)]
mx, my = [mr * f(mazr) for f in (np.sin, np.cos)]

theta = np.linspace(0, twopi, 201)
xc, yc = [f(theta) for f in (np.cos, np.sin)]

if True:
    plt.figure()
    plt.plot(sx, sy, '-y')
    plt.plot(mx, my, '-g')
    plt.plot(xc, yc, '-k', linewidth=2)
    plt.plot([0], [0], 'ok')
    plt.xlim(-1.8, 1.8)
    plt.ylim(-1.8, 1.8)
    plt.text(-0.15, 0.1, 'Zenith', fontsize=16)
    plt.text(-0.02, 1.04, 'N', fontsize=16)
    #plt.text(-1.6, 1.0, 'North is up', fontsize=16)
    plt.show()


saltd, sazd, maltd, mazd = [x.degrees for x in (salt, saz, malt, maz)]
mbreaks = np.where(np.abs(mazd[1:]-mazd[:-1] < 0))
sbreaks = np.where(np.abs(sazd[1:]-sazd[:-1] < 0))

saltd, sazd, maltd, mazd = [x[:-1] for x in (saltd, sazd, maltd, mazd)]
sazd[sbreaks] = np.nan
mazd[mbreaks] = np.nan

if True:
    plt.figure()
    plt.plot(sazd, saltd, '-y')
    plt.plot(mazd, maltd, '-g')
    plt.plot([0, 360], [0, 0], '-k')
    plt.xlabel('azimuth (degs)', fontsize=18)
    plt.ylabel('altitude/elevation (degs)', fontsize=18)
    plt.xlim(0, 360)
    plt.show()

if True:
    plt.figure()
    hours = minutes[:-1]/60.

    plt.subplot(2, 1, 1)

    plt.plot(hours, sazd, '-y')
    plt.plot(hours, mazd, '-g')
    # plt.xlabel('UTC hours', fontsize=16)
    plt.ylabel('azimuth (degs)', fontsize=16)

    plt.subplot(2, 1, 2)

    plt.plot(hours, saltd, '-y')
    plt.plot(hours, maltd, '-g')
    plt.xlabel('UTC hours', fontsize=16)
    plt.ylabel('altitude/elevation (degs)', fontsize=16)

    plt.suptitle("Haifa 24/9/2018", fontsize=16)

    plt.show()