Cálculo de la posición del sol en ECI

Si sabe cómo usar Python, puede hacer algo como:

import astropy as ap
ap.coordinates.get_sun(time)

Devolverá una posición en coordenadas geocéntricas para el tiempo dado.

La respuesta a continuación es una combinación de la primera respuesta con una verificación cruzada de Wikipedia (específicamente para la oblicuidad del plano de la eclíptica) y aquí (la fórmula para la longitud de la eclíptica del Sol en la primera respuesta usa 0.918994643 para multiplicar sin(2 * g* pi) en el término final en lugar de 0,020, por lo que utilicé el siguiente factor, pero no estoy seguro de cuál es el correcto). Además, la dirección del vector unitario no da la distancia, por lo que una respuesta completa debe escalar el vector unitario de la siguiente manera:

Primero, calcule la compensación juliana de J2000 como un número decimal de días d. Luego, realice los siguientes pasos, en orden:

1. Calcular parámetros
2. Calcular vector direccional unitario
3. Calcular la distancia al sol y escalar el vector unitario

NOTA: ¡Todas las llamadas a funciones trigonométricas (seno, coseno) utilizadas a continuación son las versiones que aceptan radianes, no grados! Casi toda la literatura en este campo usa grados en lugar de radianes porque son más intuitivos para los lectores de las publicaciones. Además, la razón de dividir por 36525 es para obtener un número decimal de siglos, dado que el año juliano tiene exactamente 365,25 días.

1. Calcular parámetros

L= 280.4606184 + [(36000.77005361 / 36525) * d] también conocido como longitud media, en grados

g= 357.5277233 + ([35999.05034 / 36525] * d) también conocida como anomalía media, en grados

p= L+ [1.914666471 * sin( g* pi/ 180)] + [0.918994643 * sin(2* g* pi/ 180)] también conocido como longitud eclíptica lambda, en grados

q= 23.43929 - ((46.8093/3600) * ( d/ 36525)) también conocido como oblicuidad del plano de la eclíptica epsilon, en grados

2. Calcular el vector direccional unitario en coordenadas ECI

u_x= cos( p* pi/ 180)

u_y= cos( q* pi/ 180) * pecado( p* pi/ 180)

u_v= pecado( q* pi/180) * pecado( p* pi/180)

3. Calcular la distancia al sol y escalar el vector unitario

a= 1.000140612 - [0.016708617 * cos( g* pi/ 180)] - [0.000139589 * cos(2 * g* pi/ 180)] también conocido como distancia desde el centro de la Tierra al centro del Sol en unidades astronómicas (AU)

m= a* 149597870700 también conocido como distancia de centro a centro de la Tierra al Sol en metros

v= m* u_v(o a* u_v) es la distancia al sol en metros (o en AU)

public static Vector3 GetSunDirection(DateTime time)
        {
            time = time.ToUniversalTime();   
            double JD = 367*time.Year-Math.Floor(7.0*(time.Year+Math.Floor((time.Month+9.0)/12.0))/4.0)+Math.Floor(275.0*time.Month/9.0)+time.Day+1721013.5+time.Hour/24.0+time.Minute/1440.0+time.Second/86400.0;
            double pi = 3.14159265359;
            double UT1 = (JD-2451545)/36525;
            double longMSUN = 280.4606184+36000.77005361*UT1;
            double mSUN = 357.5277233+35999.05034*UT1;
            double ecliptic = longMSUN+1.914666471*Math.Sin(mSUN*pi/180)+0.918994643*Math.Sin(2*mSUN*pi/180);
            double eccen = 23.439291-0.0130042*UT1;

            double x = Math.Cos(ecliptic*pi/180);
            double y = Math.Cos(eccen*pi/180)*Math.Sin(ecliptic*pi/180);
            double z = Math.Sin(eccen*pi/180)*Math.Sin(ecliptic*pi/180);

            return new Vector3(x, y, z);
        }

        public static Vector3 GetSun(DateTime time)
        {
            double sunDistance = 0.989 * 1.496E+8;
            var sunPosition = GetSunDirection(time.ToUniversalTime());

            sunPosition = sunPosition * sunDistance;

            return sunPosition;
        }