brandon-rhodes / pyephem

Scientific-grade astronomy routines for Python
Other
781 stars 121 forks source link

Satellite elong is not computed #79

Closed dmopalmer closed 8 years ago

dmopalmer commented 9 years ago

When a satellite position is computed, the elong (angular separation from the Sun) is not computed and is left with a value of 0.

(The elong is directly related to phase angle, hence brightness.)

albertw commented 9 years ago

It doesn't look as if libastro makes any effort to compute the elongation for earth satellites.

The elongation function is in circum.c for solar system bodies, and adding that logic explicitly into earthsat.c does give you an elongation. I haven't checked to see if the elongation reported is actually correct mind you!

$ git diff
diff --git a/libastro-3.7.6/earthsat.c b/libastro-3.7.6/earthsat.c
index f3c661a..a751d29 100644
--- a/libastro-3.7.6/earthsat.c
+++ b/libastro-3.7.6/earthsat.c
@@ -66,6 +67,8 @@ static void GetBearings (double SatX, double SatY, double SatZ,
 static int Eclipsed (double SatX, double SatY, double SatZ,
     double SatRadius, double CrntTime);
 static void InitOrbitRoutines (double EpochDay, int AtEod);
+static void elongation (double lam, double bet, double lsn, double *el);

 #ifdef USE_ORBIT_PROPAGATOR 
 static void GetPrecession (double SemiMajorAxis, double Eccentricity,
@@ -124,6 +127,13 @@ static double SunEpochTime,SunInclination,SunRAAN,SunEccentricity,
 static double SinPenumbra,CosPenumbra;

+static void
+elongation (double lam, double bet, double lsn, double *el)                                                            
+{
+        *el = acos(cos(bet)*cos(lam-lsn));
+        if (lam>lsn+PI || (lam>lsn-PI && lam<lsn)) *el = - *el;
+}
+
 /* given a Now and an Obj with info about an earth satellite in the es_* fields
  * fill in the s_* sky fields describing the satellite.
  * as usual, we compute the geocentric ra/dec precessed to np->n_epoch and
@@ -146,6 +156,9 @@ obj_earthsat (Now *np, Obj *op)
        double dtmp;
        double CrntTime;
        double ra, dec;
+        double lsn, rsn;        /* true geoc lng of sun, dist from sn to earth*/
+        double lam, bet;        /* geocentric ecliptic long and lat */
+        double el;              /* elongation */

 #ifdef ESAT_TRACE
        printf ("\n");
@@ -212,6 +225,13 @@ obj_earthsat (Now *np, Obj *op)

        op->s_eclipsed = Eclipsed(SatX,SatY,SatZ,Radius,CrntTime);

+        sunpos (mjed, &lsn, &rsn, NULL);
+        eq_ecl (mjed, ra, dec, &bet, &lam);
+        elongation (lam, bet, lsn, &el);
+        el = raddeg(el);
+        op->s_elong = (float)el;
+
 #ifdef ESAT_TRACE
        printf ("CrntTime = %g\n", CrntTime);
        printf ("SatX = %g\n", SatX);
brandon-rhodes commented 8 years ago

To keep PyEphem manageable, I have generally been avoiding changes to the underlying libastro library because then it is far harder to track changes in the original library when a new version comes out. For the moment, since this is not a feature of libastro, I am going to close this issue and not plan on trying to add it to PyEphem. I will instead spend that time adding an elongation function to Skyfield, so that going forward we are not stuck with the limitations of this C library!