Define a ground area target

I will show you the types of satellites and ground stations I have defined in STK, as well as the report results output by STK.
At the same time, I will display the results of running the orekit code, and may I ask why there are differences. If it’s convenient for you, please help me modify the specific code so that I can have a clear understanding of the entire usage process. Thank you!

source code:
import java.io.File;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
import org.hipparchus.ode.events.Action;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.models.earth.tessellation.EllipsoidTessellator;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.events.GeographicZoneDetector;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;

/** Orekit tutorial for special event detection.

  • This tutorial shows how to easily check for visibility between a satellite and a ground station.

  • @author Pascal Parraud
    */
    public class RectangularRegionOnEarth {

    /** Private constructor for utility class. */
    private RectangularRegionOnEarth() {
    // empty
    }

    /** Program entry point.

    • @param args program arguments (unused here)
      */
      public static void main(final String args) {
      try {
      // configure Orekit
      final File home = new File(System.getProperty(“user.home”));
      final File orekitData = new File(home, “orekit-data”);
      if (!orekitData.exists()) {
      System.err.format(Locale.US, “Failed to find %s folder%n”,
      orekitData.getAbsolutePath());
      System.err.format(Locale.US, “You need to download %s from %s, unzip it in %s and rename it ‘orekit-data’ for this tutorial to work%n”,
      “orekit-data-master.zip”, “https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip”,
      home.getAbsolutePath());
      System.exit(1);
      }
      final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
      manager.addProvider(new DirectoryCrawler(orekitData));

       //  Initial state definition : date, orbit
       final AbsoluteDate initialDate = new AbsoluteDate(2023, 10, 26, 8, 00, 00.000, TimeScalesFactory.getUTC());
       final double mu =  3.986004415e+14; // gravitation coefficient
       final Frame inertialFrame = FramesFactory.getEME2000(); // inertial frame for orbit definition
       final Vector3D position  = new Vector3D(-7518565.475, 4279265.668, 4816446.308);
       final Vector3D velocity  = new Vector3D(-5463.904, -5533.9672, 50.719);
       final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
       final Orbit initialOrbit = new KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, mu);
      
       // Propagator : consider a simple Keplerian motion 
       final Propagator kepler = new KeplerianPropagator(initialOrbit);
      
       // Earth and frame
       final Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
       final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                                                    Constants.WGS84_EARTH_FLATTENING,
                                                    earthFrame);
      
       //创建一个地面区域
       GeodeticPoint point1 = new GeodeticPoint(Math.toRadians(0), Math.toRadians(30), 0.0);
       GeodeticPoint point2 = new GeodeticPoint(Math.toRadians(30), Math.toRadians(30), 0.0);
       GeodeticPoint point3 = new GeodeticPoint(Math.toRadians(30), Math.toRadians(0), 0.0);
       GeodeticPoint point4 = new GeodeticPoint(Math.toRadians(0), Math.toRadians(0), 0.0);
      
       // 创建表示地面区域的SphericalPolygonsSet
       SphericalPolygonsSet region = EllipsoidTessellator.buildSimpleZone(1.0e-10, point1, point2, point3, point4);
      
       // Event definition
       final double maxcheck  = 60;
       final double threshold =  0.001;
       final double elevation = FastMath.toRadians(0.0);
        
       final GeographicZoneDetector zoneDetector = new GeographicZoneDetector(maxcheck, threshold, earth, region, elevation)
       		.withHandler((s, detector, increasing) -> {
                   if (increasing) {
                       System.out.println("Visibility begins at " + s.getDate());
                   } else {
                       System.out.println("Visibility ends at " + s.getDate());
                       System.out.println();
                   }
                   return Action.CONTINUE; // Continue processing the event
               });
      
       // Add event to be detected
       kepler.addEventDetector(zoneDetector);
      
       // Propagate from the initial date to the first raising or for the fixed duration
       final SpacecraftState finalState = kepler.propagate(initialDate.shiftedBy(360000.));
      
       System.out.println(" Final state : " + finalState.getDate().durationFrom(initialDate));
      

      } catch (OrekitException oe) {
      System.err.println(oe.getLocalizedMessage());
      }
      }

}

You define an area that is almost rectangular. It is in fact skewed because the 30° range in longitude at 30° North latitude represents a smaller distance than the same 30° range at equator. Lets consider a satellite crossing the equator at longitude l. If l<0°, then the closest boundary point will be the South-West corner at 0°N, 0°E. So the signed distance to the boundary will be |l| = -l. Hence the g function that triggers the event will be positive if -l-5° > 0, i.e. l < -5°, zero if l=-5° and positive if l>-5°. Therefore the trigger will occur at longitude -5°. On the other side, if l > 30°, the closest boundary point will be will be the South-East corner at 0°N, 30°E. So the signed distance to the boundary will be l-30°. Hence the g function that triggers the event will be positive if (l-30°)-5° > 0, i.e. if l > 35°, zero if l=35° and negative if l<35°. Therefore the trigger will occur at longitude 35°. Basically, a positive margin enlarges the detection zone and a negative margin reduces the detection zone. As the zone is skewed (smaller on the North boundary than at equator) and as computations are done using exact spherical angles and great circles, it is a skewed trapezoid with rounded corners.

Beware that the g function is positive when the satellite is outside the area and negative when it is inside the area. So when the event handler is called with increasing set to true, it means you are leaving the region, not entering it. You have the revert the if statement in your code.

I am not sure you are computing the same thing with STK and Orekit. In the Orekit case, there are no sensors. If you have one tiny almost point sensor in STK, you also have to make sure it is pointing along the nadir direction, taking Earth flattening into account, so it is equivalent to computing the sub-satellite point. Otherwise, the satellite may see the zone despite it is not above it, just because it looks slightly on the side. Another point that may induce differences is how the Nroth boundary of your zone is handled. In Orekit, we use always great circles, i.e. the boundary would cross latitude 30° at both endpoints, but will extend slightly above 30° latitude in the middle, it is not a parallel to equator. I don’t know if STK uses the same definition and how it handles two points at the same latitude.

As we said earlier, ground zone and satellite field of view (which is what you seem to want to do with STK) is done by FootprintOverlapDetector, but it is tricky to set up and computationally intensive in Orekit, we have to improve it.

So now I want to define a satellite and a ground area, and obtain the time period of the satellite’s sub satellite point trajectory in the ground area. How do I write the specific code?