phase_angle error with EarthSatellite

[smp55]

The new method phase_angle seems to have a bug when called from an EarthSatellite object. When callling EarthSatellite.at(time).phase_angle(sun), I get an error for line 550 in phase_angle:

v = u - s.position.au + self.center_barycentric.position.au
AttributeError: 'NoneType' object has no attribute 'position'

[brandon-rhodes]

Could you give a working example, along with the angle you would like returned? A phase angle is only defined from the point of view of a particular observer, so I need to see your earlier code to see how you are defining the position of the observer. (Or — are you not providing an observer at all? In that case maybe phase_angle() needs a more informative error message for this case; let me know!)

[smp55]

Good point, I have an observer defined, but I realized that isn't actually included in the EarthSatellite object. I use it elsewhere by creating a line-of-sight object by subtracting the observer position from the EarthSatellite object. But phase_angle(sun) is a method of the EarthSatellite object itself. Should there be another argument for that function? I do in fact want the solar phase angle (sun->sat->observer), and this method seemed likely. I suppose I could get it through other vector manipulations as well.

[brandon-rhodes]

Happily, phase_angle() is not exactly a method of the EarthSatellite, or else you would say sat.phase_angle() after defining a satellite sat. Rather, it's a method of the position returned by .at(). Maybe you could try giving the position both your observer and the satellite:

(sat - observer).at(t).phase_angle(sun)

See if that works.

[smp55]

Same error, unfortunately! Here's a minimum example (requires de440.bsp ephemeris file):

import pytz
from datetime import datetime
from skyfield.api import load
from skyfield.toposlib import wgs84
from skyfield.sgp4lib import EarthSatellite

ts = load.timescale(builtin=True)
planet_ephem_loc = 'files/de440.bsp'
planets = load(planet_ephem_loc)
Sun = planets['Sun']

tle = ['1 07734U 75027A   23304.59378895 -.00000022  00000-0  73408-4 0  9999',
       '2 07734 114.9879 167.4147 0025812 288.6764  71.1513 14.18833704513119']

sat = EarthSatellite(tle[0], tle[1])
obs = wgs84.latlon(34.224915, -118.056151, elevation_m=1.733)
t = ts.utc(datetime(2023, 10, 31, 13, 44, tzinfo=pytz.utc))

(sat-obs).at(t).phase_angle(Sun)

[brandon-rhodes]

Thank you for the working example! It looks like phase_angle() was only written and tested against big things like planets and the Moon, which in Skyfield should always be submitted to the .observe() method. In that case the computation works — though it requires both the observer and the satellite to know where the Earth itself is, since .observe() needs to know where in the solar system both the observer and the target are:

import pytz
from datetime import datetime
from skyfield.api import load
from skyfield.toposlib import wgs84
from skyfield.sgp4lib import EarthSatellite

ts = load.timescale(builtin=True)
planet_ephem_loc = 'files/de440.bsp'
planets = load(planet_ephem_loc)
sun = planets['Sun']
earth = planets['Earth']

tle = ['1 07734U 75027A   23304.59378895 -.00000022  00000-0  73408-4 0  9999',
       '2 07734 114.9879 167.4147 0025812 288.6764  71.1513 14.18833704513119']

sat = earth + EarthSatellite(tle[0], tle[1])
obs = earth + wgs84.latlon(34.224915, -118.056151, elevation_m=1.733)
t = ts.utc(datetime(2023, 10, 31, 13, 44, tzinfo=pytz.utc))

obs.at(t).observe(sat).phase_angle(sun)

It's an interesting question how phase_angle() might be made to work when taking the shortcut of (sat-obs), which computes a quick relative position by simply subtracting two vectors centered at the center of the Earth. The resulting position has no idea of where it's located relative to the Solar System barycenter. I'll try to think through how phase_angle() might work in those circumstances.

Meanwhile, I'll at least try in the next week or two to sit down and improve that error message! Thanks for letting me know of this use case; I've never had someone try to compute the phase angle for a satellite before, if I recall.

[smp55]

Thanks so much, that will work! It might not be related, but I notice that the .is_sunlit(ephemeris) function works with an EarthSatellite without adding the Earth to its position first. Seems like it would be a similar calculation.