I decided to make my own detector since the events are connected to the orbiting satellite.
Thanks for the guidance.
I attach the code here. I got within 2-3 seconds of STK for a polar orbiting satellite at approx 630km.
It works, but don’t take as super scientific answer.
-Tommy
class InsideConeDetector(PythonAbstractDetector):
def __init__(self, time_step: float, half_cone_angle: float,
frame: TopocentricFrame, handler: EventHandler = None):
self._frame = frame
self._half_cone_angle = half_cone_angle
if handler is None:
handler = StopOnEvent()
super().__init__(float(time_step),
PythonAbstractDetector.DEFAULT_THRESHOLD,
PythonAbstractDetector.DEFAULT_MAX_ITER,
handler)
def init(self, state: SpacecraftState, date: AbsoluteDate):
pass
def g(self, state: SpacecraftState):
# https://en.wikipedia.org/wiki/Geostationary_orbit
geo_altitude_above_equator_m = 35786000 # meters
sat_elevation = self._frame.getElevation(state.getPosition(), state.getFrame(),
state.getDate())
sat_altitude = EARTH.transform(state.getPosition(), state.getFrame(), state.getDate()).getAltitude()
sat_latitude = EARTH.transform(state.getPosition(), state.getFrame(), state.getDate()).getLatitude()
#https://gis.stackexchange.com/questions/20200/how-do-you-compute-the-earths-radius-at-a-given-geodetic-latitude
earth_a = 6378137
earth_b = 6356752.3142
earth_radius_at_latitude_m = sqrt((pow(earth_a*earth_a*cos(abs(sat_latitude)),2) +
pow(earth_b*earth_b*sin(abs(sat_latitude)),2))/
(pow(earth_a*cos(abs(sat_latitude)),2) +
pow(earth_b*sin(abs(sat_latitude)),2)))
# Half cone angle at 300km and 0km altitude
half_cone_angle_at_300km = asin((earth_radius_at_latitude_m + 300e3)/(geo_altitude_above_equator_m+earth_radius_at_latitude_m))
half_cone_angle_at_0km = asin((earth_radius_at_latitude_m)/(geo_altitude_above_equator_m+earth_radius_at_latitude_m))
max_elevation = atan(earth_radius_at_latitude_m/(earth_radius_at_latitude_m+sat_altitude))
hypotenus_sat = self._frame.getRange(state.getPosition(), state.getFrame(),state.getDate()) # hypotenus, slantrange
# Distance from line, between geo relay and nadir point, to satellite
x_sat = hypotenus_sat * cos(abs(sat_elevation))
# Height of satellite, above or below nadir point of geo relay
y_sat = hypotenus_sat * sin(abs(sat_elevation))
# Calculate radius of cone at satellite altitudes
if (sat_elevation > 0):
cone_radius_at_sat_altitude = tan(radians(self._half_cone_angle))*(geo_altitude_above_equator_m-y_sat)
cone_radius_at_300km = tan(half_cone_angle_at_300km)*(geo_altitude_above_equator_m-y_sat)
cone_radius_at_0_altitude = tan(half_cone_angle_at_0km)*(geo_altitude_above_equator_m-y_sat)
else:
cone_radius_at_sat_altitude = tan(radians(self._half_cone_angle))*(geo_altitude_above_equator_m+y_sat)
cone_radius_at_300km = tan(half_cone_angle_at_300km)*(geo_altitude_above_equator_m+y_sat)
cone_radius_at_0_altitude = tan(half_cone_angle_at_0km)*(geo_altitude_above_equator_m+y_sat)
# Check that satellite is not in earth shadow marked by cone_radius at 0km altitude
# and that satellite is not inside 300 km grazing altitude
if ((sat_elevation < -max_elevation and x_sat < cone_radius_at_0_altitude) or
(sat_elevation < -max_elevation and x_sat < cone_radius_at_300km)):
_log.debug("Dato: %s, elevation: %.2f (max: %.2f), y_sat: %.3f km, x_sat %.3f =? cone_radius_at_0km: %.3f km, cone_radius_at_300km: %.3f km" %
(state.getDate().toString(), degrees(sat_elevation), degrees(max_elevation), y_sat/1e3, x_sat/1e3,
cone_radius_at_0_altitude/1e3, cone_radius_at_300km/1e3))
return -1.0
res = cone_radius_at_sat_altitude - x_sat
#if (res > 0):
# _log.debug("Dato: %s, elevation: %.2f, y_sat: %.3f km, x_sat %.3f =? cone_radius: %.3f km, res = %.3f" %
# (state.getDate().toString(), degrees(sat_elevation), y_sat/1e3, x_sat/1e3, cone_radius_at_sat_altitude/1e3, res))
#else:
# _log.debug("Elevation: %.2f, y_sat: %.3f km, x_sat %.3f =? cone_radius: %.3f km, res = %.3f" %
# (degrees(sat_elevation), y_sat/1e3, x_sat/1e3, cone_radius_at_sat_altitude/1e3, res))
return res