Document how to compute the phase angle of a sunlit object (like an Earth satellite)

[mworion]

Hi, I wonder how to calculate the phase angle of a sunlit satellite in a performant way. I learned that for satellites and the precision the observe() method could be avoided. So my basic attempt was:

        sun = ephemeris['sun']
        earth = ephemeris['earth']

        vecObserverSat = (sat - loc).at(tEv)
        vecObserverSun = (sun - (earth + loc)).at(tEv)
        phase = vecObserverSat.separation_from(vecObserverSun)

for a given time tEv and an observer position loc on earth. Does this simple attempt make sense or is it wrong? I just winding up my head where a mistake in thinking might be in there. It seems too simple. Any feedback appreciated. Michel

PS: perhaps when sorting this out a magnitude_satellite feature could be added as well.

[brandon-rhodes]

Here's how phase angle is calculated when doing planetary magnitudes:

https://github.com/skyfielders/python-skyfield/blob/b438e9fb8e1cbc6ea5b224d06f9a7dba6f0fff53/skyfield/magnitudelib.py#L80

I would want an angle meeting at the satellite, so I would set up the positions a little differently. Also, the Sun's position can be treated as the Solar System Barycenter for the purposes of magnitude, I suspect? So I would compute the separation between these two vectors:

sat - loc       # loc → sat vector
sat + earth   # sun → sat vector

Does that give an angle that looks similar to the one you're after?

I would indeed be interested to see an approach for computing the magnitude of a satellite; I'm not familiar with a general formulae for it.

[mworion]

@brandon-rhodes : Thanks for the quick response! For the formulae to get an estimation of the apparent magnitude of a satellite I did some search work:

https://astronomy.stackexchange.com/questions/28744/calculating-the-apparent-magnitude-of-a-satellite/28765#28765 https://www.researchgate.net/publication/268194552_Large_phase_angle_observations_of_GEO_satellites https://amostech.com/TechnicalPapers/2013/POSTER/COGNION.pdf https://apps.dtic.mil/dtic/tr/fulltext/u2/785380.pdf

Unfortunately it could be only an estimation based data from real observations. Secondly as you mentioned in your planetary part as well, the real material, construction, size and reflection differs from satellite to satellite. So you need somehow the "intrinsic" magnitude of the satellite. They are not included in main satellite databases like celestak et.al. Empiric values of -1.3 (sometimes -1.8) mag are referenced in this papers. So we need definitely a parameter for this value to overwrite standard settings if somebody nows this value.

In overall, the work leads finally to the following function:

def calcAppMag(self, phase, satRange, intMag=-1.3)
    term1 = intMag
    term2 = +5.0 * np.log10(satRange / 1000.)
    arg = np.sin(phase) + (np.pi - phase) * np.cos(phase)
    term3 = -2.5 * np.log10(arg)
    appMag = term1 + term2 + term3
    return appMag

Question to you calculation of the phase angle. the first term is used like I had it in mind. It's a geocentric vector. I agree with your opinion that the sun's position can be treated as the Solar System Barycenter. If so I get an error message from angle_between in functions.py line 52:

[functions.py], line: 52 msg: [unsupported operand type(s) for *: 'Barycentric' and 'Barycentric']

Any idea? Michel

[brandon-rhodes]

Oh — you can only use angle_between() with raw NumPy vectors; I should have mentioned that you can ignore that call. Please continue to use separation_from() when you're passing higher level Skyfield position objects.

[mworion]

Ok, this had to come to my mind by myself. Sorry to ask. The result now is in function:

def calcSatSunPhase(sat, loc, ephemeris, tEv):
  earth = ephemeris['earth']

  vecObserverSat = (sat - loc).at(tEv)
  vecSunSat = (earth + sat).at(tEv)
  phase = vecObserverSat.separation_from(vecSunSat)

  return phase

works nicely and performant (as I calculate all - up to 4000 - satellites in the database) for visualization.

[brandon-rhodes]

I'm glad it worked! I'll look at adding some documentation explaining this maneuver before the next Skyfield release. I'll be traveling next week, but I'll plan to write it up as soon as I return.