1
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Unfortunately I am using c# not python. $\endgroup$ – Nick Brown Dec 19 '18 at 17:24
  • $\begingroup$ Then, i don't know how to help you ... sorry $\endgroup$ – Yoann A. Dec 19 '18 at 20:51
2
$\begingroup$

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)

$\endgroup$
  • $\begingroup$ If anyone has any more precise values for how to calculate the ecliptic longitude in radians, I would appreciate an edit. $\endgroup$ – brethvoice Jul 29 at 16:14
  • $\begingroup$ Also the formula for the obliquity of the ecliptic plane only has one significant digit, so any help there would also be appreciated! $\endgroup$ – brethvoice Jul 29 at 17:52
  • $\begingroup$ 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). $\endgroup$ – brethvoice Jul 29 at 18:02
  • $\begingroup$ 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... $\endgroup$ – brethvoice Jul 29 at 18:19
  • $\begingroup$ Also apparently the starting point for Earth's tilt at J2000.0 is 23.439291... degrees according to the same article...updating formula again. $\endgroup$ – brethvoice Jul 29 at 18:23
0
$\begingroup$
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;
        }
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.