Hi,
I’m trying to determine satellite access times to an arbitrary point on Earth, considering the satellite’s Field of View (FOV). I want to use a square FOV of 1° (0.5° half FOV) with 30º maximum off-nadir angle. However, I’ve approximated the FOV directly as a 1°x61° rectangle without off-nadir (in the code you can see the half width and half height values)
I’m encountering issues with the access calculations. When using the 1°x61° FOV, I only detect a single access event over a 10-day period. This seems improbable, and STK confirms my suspicion that the access frequency should be higher. I suspect the problem lies in my configuration of the DoubleDihedraFieldOfView, the target zone, or the FootprintOverlapDetector. All the propagator part is working correctly.
I’ve reviewed the available examples and forum discussions, but I haven’t been able to pinpoint the error in my setup. Could you please help me identify what I might be doing incorrectly in the configuration process?
Here you can find my code (I’ll skip the imports, initial_setup, logger…):
# 0. Set up constants
initial_state = [
-4305.630472e3,
5417.737575e3,
10.331028e3,
0.801232e3,
0.605594e3,
7.524735e3
] # m/s
mass = 2.455 # kg
drag_coeff = 2.2
drag_area = 0.0025 # m²
srp_coeff = 1.8
srp_area = 0.0025 # m²
start_date = AbsoluteDate(2025, 3, 16, 13, 28, 26.724, TimeScalesFactory.getUTC())
end_date = AbsoluteDate(
2025, 3, 25, 0, 0, 00.000, TimeScalesFactory.getUTC()
) # propagation for 5 days
step_size = 10.0 # seconds
file_path = "data/output/orekit_data.txt"
# 1. Set up Earth parameters
earth = CelestialBodyFactory.getEarth()
earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earthRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
earthFlattening = Constants.WGS84_EARTH_FLATTENING
earthShape = OneAxisEllipsoid(earthRadius, earthFlattening, earthFrame)
# 2. Set up Earth gravity model
gravity_provider = GravityFieldFactory.getNormalizedProvider(70, 70)
gravity_model = HolmesFeatherstoneAttractionModel(
FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravity_provider
)
# 3. Create initial orbit - EME2000 frame
pvCoordinates = PVCoordinates(
Vector3D(initial_state[0], initial_state[1], initial_state[2]),
Vector3D(initial_state[3], initial_state[4], initial_state[5]),
)
initial_orbit = CartesianOrbit(
pvCoordinates, FramesFactory.getEME2000(), start_date, Constants.WGS84_EARTH_MU
)
# 4. Set up numerical propagator
min_step = 1e-3
max_step = 500.0
init_step = 60.0
# For extracting tolerances if is needed, for now I hard-code 1e-13 when
# construting the integrator instance
# position_tolerance = 1e-4
# velocity_tolerance = 1e-5
integrator = DormandPrince853Integrator(min_step, max_step, 1e-13, 1e-13)
integrator.setInitialStepSize(init_step)
propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)
propagator.setInitialState(SpacecraftState(initial_orbit, mass))
propagator.setMu(Constants.WGS84_EARTH_MU)
# 5. Add force models
# 5.1. Add gravity model
propagator.addForceModel(gravity_model)
# 5.2. Add atmospheric drag model
cssiSpaceWeatherData = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
atmosphere = NRLMSISE00(cssiSpaceWeatherData, CelestialBodyFactory.getSun(), earthShape)
dragSC = IsotropicDrag(drag_area, drag_coeff)
drag_force = DragForce(atmosphere, dragSC)
propagator.addForceModel(drag_force)
# 5.3. Add solar radiation pressure model
sun = CelestialBodyFactory.getSun()
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(srp_area, srp_coeff)
srpForce = SolarRadiationPressure(sun, earthShape, isotropicRadiationSingleCoeff)
propagator.addForceModel(srpForce)
# 5.4 Add third body forces models
propagator.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getSun()))
propagator.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getMoon()))
# Define target point on Earth
target_pos = Vector3D(earthRadius * 0.75, earthRadius * 0.03, earthRadius * 0.66)
meridian = Vector3D(
target_pos.getX(),
target_pos.getY() + 10000.0, # Desplazamiento de 10 km al norte
target_pos.getZ()
)
tolerance = 1.0e-6 # Tolerancia más precisa
zone_radius = radians(0.001) # Aproximadamente 0.05 radianes
zone = SphericalPolygonsSet(target_pos, meridian, zone_radius, 8, tolerance)
half_height = radians(0.5)
half_width = radians(30.0)
t = FramesFactory.getEME2000().getTransformTo(FramesFactory.getITRF(IERSConventions.IERS_2010, True), start_date)
pvCoordinates = t.transformPVCoordinates(SpacecraftState(initial_orbit, mass).pVCoordinates)
# Unit nadir vector
pos_vector = pvCoordinates.position
nadir_vector = pos_vector.negate()
unit_nadir_vector = nadir_vector.scalarMultiply(1.0/nadir_vector.getNorm())
# Unit velocity vector
vel_vector = pvCoordinates.velocity
unit_vel_vector = vel_vector.scalarMultiply(1.0/vel_vector.getNorm())
# Unit across vector
unit_across_vector = unit_vel_vector.crossProduct(unit_nadir_vector)
# Create imager FOV
imagerFOV = DoubleDihedraFieldOfView(
unit_nadir_vector,
unit_across_vector,
half_height,
unit_vel_vector,
half_width,
0.0
)
sampling_step = 100.
targetDetector = FootprintOverlapDetector(imagerFOV, earthShape, zone, sampling_step)
logger = EventsLogger()
propagator.addEventDetector(logger.monitorDetector(targetDetector))
propagator.propagate(end_date)
mylog = logger.getLoggedEvents()
for event in mylog:
print(event.date.date.toString())