Calculating the sun’s position in ECI

Question

I’m trying to calculate visible ISS passes and I need to calculate if the ISS is behind the earth by “drawing” a line between the ISS and the sun and seeing if it crosses earth. I have the ECI position of the ISS, but I don’t have the position of the sun. How do I calculate the suns ECI position given a date and time?

Edit: I need to use c# for this so I can’t use python libraries

Answer 1

If you know how to use python, you can do something like :

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

It will return a position in geocentric coordinates for the given time.

Answer 2

The answer below is a combination of the first answer with a cross-check from Wikipedia (for the obliquity of the ecliptic plane specifically) and here (the formula for ecliptic longitude of the Sun in the first answer uses 0.918994643 to multiply sin(2 * g * pi) in the final term instead of 0.020, so I used the factor below, but I am not sure which is correct). Also, the unit vector's direction does not give the distance, so a complete answer needs to scale the unit vector as follows:

First, calculate the Julian offset from J2000 as a decimal number of days d. Then, do the following steps, in order:

1. Calculate parameters
2. Calculate unit directional vector
3. Calculate distance to sun and scale the unit vector

NOTE: All trigonometric function calls (sine, cosine) used below are the versions which accept radians, not degrees! Almost all of the literature in this field uses degrees rather than radians because they are more intuitive for readers of the publications. Also, the reason for dividing by 36525 is to obtain a decimal number of centuries, given that the Julian year is 365.25 days exactly.

1. Calculate parameters

L = 280.4606184 + [(36000.77005361 / 36525) * d] a.k.a. mean longitude, in degrees

g = 357.5277233 + ([35999.05034 / 36525] * d) a.k.a. mean anomaly, in degrees

p = L + [1.914666471 * sin(g * pi / 180)] + [0.918994643 * sin(2*g * pi / 180)] a.k.a. ecliptic longitude lambda, in degrees

q = 23.43929 - ((46.8093/3600) * (d / 36525)) a.k.a. obliquity of ecliptic plane epsilon, in degrees

2. Calculate unit directional vector in ECI coordinates

u_x = cos(p * pi / 180)

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

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

3. Calculate distance to sun and scale the unit vector

a = 1.000140612 - [0.016708617 * cos(g * pi / 180)] - [0.000139589 * cos(2 * g * pi / 180)] a.k.a distance from Earth's center to Sun's center in astronomical units (AU)

m = a * 149597870700 a.k.a. center-to-center distance from Earth to Sun in meters

v = m * u_v (or a * u_v) is distance to sun in meters (or in AU)

Answer 3

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;
        }