# Calculating the sun’s position in ECI

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

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.

• Unfortunately I am using c# not python. Dec 19, 2018 at 17:24
• Then, i don't know how to help you ... sorry Dec 19, 2018 at 20:51

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)

• If anyone has any more precise values for how to calculate the ecliptic longitude in radians, I would appreciate an edit. Jul 29, 2020 at 16:14
• Also the formula for the obliquity of the ecliptic plane only has one significant digit, so any help there would also be appreciated! Jul 29, 2020 at 17:52
• Alright I just corrected the formula for obliquity of the ecliptic plane, which was also actually in degrees, and made it have two significant digits (47 arcseconds per century of change). Jul 29, 2020 at 18:02
• Apparently a better number is 46.8093 arcseconds per century starting at J2000.0 according to en.wikipedia.org/wiki/Axial_tilt#Short_term, which is good for the next few thousand years? Updating above... Jul 29, 2020 at 18:19
• Also apparently the starting point for Earth's tilt at J2000.0 is 23.439291... degrees according to the same article...updating formula again. Jul 29, 2020 at 18:23
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;
}