Matemáticas para calcular la longitud terrestre directamente bajo el sol con el tiempo

Estoy tratando de calcular la longitud en la tierra donde es mediodía en algún momento. (Es decir, la longitud que es coplanar con el plano definido por el sol y el eje de la tierra).

Aquí está mi código Python:

from math import sin, cos, tan, atan, pi
def sun_longitude(when):
    """Given time in ms since 1/1/1970, return the longitude the sun is over at that moment"""
    # https://en.wikipedia.org/wiki/Position_of_the_Sun
    jdn = 2440587.5 + when / (1000.0 * 3600 * 24)
    n = jdn - 2451545.0 # 1/1/2000
    L = (280.460 + 0.9856474 * n) % 360.0
    g = (357.528 + 0.9856004 * n) % 360.0
    degtorad = 2.0 * pi / 360.0
    lambda_ = L + 1.915 * sin(g * degtorad) + 0.020 * sin(2 * g * degtorad)

    # https://en.wikipedia.org/wiki/Axial_tilt#Short_term
    T = n / (365 * 100.0) # julian centuries from 2000
    # ε = 23° 26′ 21.406″ − 46.836769″ T − 0.0001831″ T2
    epsilon = 23.4392794 - 0.780612817 * T - 5.0861e-8 * T * T;
    alpha = atan(cos(epsilon * degtorad) * tan(lambda_ * degtorad)) / degtorad

    # https://en.wikipedia.org/wiki/Sidereal_time
    # https://en.wikipedia.org/wiki/Hour_angle
    GMST = (18.697374558 + 24.06570982441908 * n) % 24.0
    print "jdn = {jdn}, n = {n}, L = {L}, g = {g}, lambda_ = {lambda_}, T = {T}, epsilon = {epsilon}, alpha = {alpha}, GMST = {GMST}".format(**locals())
    LHA = GMST * 360/24 + alpha;

    return LHA

from datetime import datetime
def sun_long_from_str(whenstr):
    when = datetime.strptime(whenstr, "%Y-%m-%d %H:%M")
    secs = (when - datetime(1970, 1, 1)).total_seconds()
    return sun_longitude(secs * 1000.0)

Soy consciente de que no estoy corrigiendo el cuadrante correcto al calcular alfa, pero tengo errores en otros lugares; por ejemplo, para el mediodía GMT del 1 de enero de 2000, obtengo 279° en lugar de cerca de 0, que es lo que esperaría.

Esto no me está dando respuestas correctas, y no sé cómo depurarlo. ¿Alguien puede encontrar mis errores o indicarme algún código de muestra razonable o un ejemplo elaborado para esto?

Porque, una vez que obtenga el algoritmo correcto, lo volveré a implementar para un dispositivo integrado, no puedo simplemente usar un paquete de implementación y no he encontrado una biblioteca que sea lo suficientemente sencilla para entender cómo puedo implementar esto. pregunta especifica

Gracias por un par de correcciones de errores. Ahora, cuando lo ejecuto, obtengo lo siguiente para estos ejemplos:

2000-01-01 12:00:
  jdn = 2451545.0, n = 0.0, L = 280.46, g = 357.528, lambda_ = 280.375680197, T = 0.0, epsilon = 23.4392794, alpha = -78.7141369122, GMST = 18.697374558
  result: 201.74648145775257

2000-03-20 07:35:
  jdn = 2451623.81597, n = 78.815972222, L = 358.144758099, g = 75.2090537484, lambda_ = 360.006175505, T = 0.00215934170471, epsilon = 23.4375937902, alpha = 0.00566598789588, GMST = 19.4596915825
  291.9010397253072

Como puede ver, la primera consulta (mediodía del 1 de enero de 2000) ahora tiene un número de día juliano correcto. Al mediodía GMT, esperaría que el sol estuviera por encima de una longitud cercana a 0, por lo que 201 es incorrecto.

La segunda consulta es la hora del equinoccio de primavera en 2000, que esperaba que pusiera a cero la parte sideral del cálculo, pero no es así.

Intenté consultar la interfaz web de NASA HORIZONS para las efemérides del sol el 1/1/2000 al mediodía GMT, tanto para las ubicaciones "Geocéntricas" como para las de Greenwich, pero no sé cómo evaluar el resultado y compararlo con mi trabajo. encima.

Geocéntrico:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00     18 45 09.36 -23 01 59.7 -26.78 -10.59 0.98332762653520  -0.0127281   0.0000 /?   0.0000

Greenwich:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00 *m  18 45 09.36 -23 02 08.3 -26.78 -10.59 0.98331613178086  -0.0166655   0.0000 /?   0.0000

Gracias de nuevo por cualquier ayuda que pueda proporcionar.

[Editado para corregir el cálculo del día juliano y el uso de grados frente a radianes para alfa y para agregar ejemplos]