Failed to find root

Orekit Python Wrapper usage
1

Hello,
I’m new in this forum, arrived last month. Since then, I used a lot orekit and the function I found the most interesting is generator.generate(AbsoluteDate, AbsoluteDate). Unfortunately, it happens that when I put a lot of satellites and a big amount of time between the 2 dates, a crash stops the code with the following error:
org.orekit.errors.OrekitException: org.orekit.estimation.measurements.generation.EventBasedScheduler$FeasibilityAdapter@4044f07b failed to find root between 2023-10-10T20:00:08.95072920450986Z (g=-1.238782448165088E-2) and 2023-10-10T20:00:09.95059701190634Z (g=2.10063578957876E-2)
[…]
Caused by: org.hipparchus.exception.MathIllegalArgumentException: interval does not bracket a root: f(0E0) = -12.387824481650878E-3, f(999.8678073964804E-3) = -13.602691220661267E-3
[…]

This does not seem to be an error on my side, since most of the time I just need to relaunch the code to get a result.
Is this error common to generator? Is this only the case through the python wrapper?
Is there a way to avoid this, rather than generating for small duration?
thanks in advance for your help
LM

2

Hi @lmirallie welcome

This may be related to issue 1164. In fact I discovered the issue exactly in the same case: generating measurements with several satellites.

Could you try using the version from the develop branch?

3

Hi,
Thanks for your answer!
I was using python, so I translated my code to java to test with the version from the develop branch.
At the point where I’m calling generate, I have an issue with the new function generate, which is void… Shall I use the old one or is it the one that is interesting to use for my problem? How do I get the results from public void generate(final AbsoluteDate start, final AbsoluteDate end) ?

By the way, if I call generate( start, end), I still get the error mentioned above.
Thanks
Louis

4

Hi @luc and @lmirallie ,

I do not believe it comes from issue 1164 because @lmirallie is using the python wrapper which is most probably the 11.3 version of orekit and not the develop branch where issue 1164 occured.

Could you detail the detectors you are using @lmirallie ?

Cheers,
Vincent

5

Hello, thanks for your question.
I’m creating a boolean detector composed of 3 detectors: one fov detector, one eclipse detector and one groundatnight detector.
Is that enough?

6

Hello @lmirallie,

Could you please send us a small runnable example of code that exhibits the buggy behavior?

Thanks,
Maxime

7

It follows for the example in “simulations of Observations” (Simulations of observations)

earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, earthFrame)
utc = TimeScalesFactory.getUTC()
satellites = []
satellites.append({"satname" : "ADEOS 2",
                   "tle_line1" : "1 27597U 02056A   23292.82900545  .00000751  00000-0  32003-3 0  9998",
                   "tle_line2" : "2 27597  98.6603 235.9627 0001869  56.8550 303.2810 14.28048440 85804"})
satellites.append({"satname"   : "Starlink-4522",
                    "tle_line1" : "1 53388U 22097A   23271.93975095  .00016889  00000-0  10762-2 0  9999",
                   "tle_line2" : "2 53388  53.2169 171.8279 0001571 122.4076 237.7068 15.08857317 63249"})
satellites.append({"satname"   : "Starlink-2733",
                   "tle_line1" : "1 48645U 21044H   23271.94144011  .00013105  00000-0  89686-3 0  9997",
                   "tle_line2" : "2 48645  53.0563 169.3229 0001526 104.7477 255.3681 15.06408749130050"})
satellites.append({"satname"   : "Starlink-1584",
                   "tle_line1" : "1 46044U 20055T   23272.14091619  .00007340  00000-0  51088-3 0  9996",
                   "tle_line2" : "2 46044  53.0554 173.4112 0001705  89.1092 271.0092 15.06419315173853"})
satellites.append({"satname"   : "Starlink-1302",
                   "tle_line1" : "1 45375U 20019R   23272.14555380  .00018655  00000-0  12690-2 0  9998",
                   "tle_line2" : "2 45375  53.0541 213.4058 0001370  95.7398 264.3748 15.06369552195202"})
satellites.append({"satname"   : "Starlink-3299",
                   "tle_line1" : "1 50831U 22001AE  23272.27601560  .00004205  00000-0  28215-3 0  9997",
                   "tle_line2" : "2 50831  53.2172 215.3839 0001265  92.3413 267.7724 15.08840448 96952"})

station = GeodeticPoint(latitude, longitude, altitude)
stationFrame = TopocentricFrame(earth, station, "Lausanne")
sun = CelestialBodyFactory.getSun()
groundStation = GroundStation(stationFrame)
t_beginning = AbsoluteDate(2023, 10, 18, 0, 22, 0.0, utc)
t_end       = AbsoluteDate(2023, 10, 20, 8, 22, 0.0, utc)

fov = DoubleDihedraFieldOfView(center, axis1, ha1, axis2, ha2, 0.0)
fovGroundDetector = GroundFieldOfViewDetector(stationFrame, fov).withHandler(ContinueOnEvent())
fovGroundDetector = fovGroundDetector.withMaxCheck(1.)
nightDetector = GroundAtNightDetector(stationFrame, sun, GroundAtNightDetector.CIVIL_DAWN_DUSK_ELEVATION,
                                      None).withHandler(ContinueOnEvent())
eclipseDetector = EclipseDetector(sun, Constants.SUN_RADIUS, earth).withPenumbra().withHandler(ContinueOnEvent())
fovAtNightDetector = BooleanDetector.andCombine([BooleanDetector.notCombine(fovGroundDetector),
                                                 nightDetector,
                                                 eclipseDetector]).withHandler(ContinueOnEvent()).withMaxCheck(1.)

generator = Generator()
step = 60.0
for sate in satellites:
    myTLE = TLE(sate["tle_line1"], sate["tle_line2"])
    propagator = TLEPropagator.selectExtrapolator(myTLE)
    sat = generator.addPropagator(propagator)                      
    azElBuilder = AngularAzElBuilder(None, groundStation, [1., 1.], [1., 1.], sat)

    scheduler = EventBasedScheduler(azElBuilder, FixedStepSelector(step, utc), generator.getPropagator(sat), fovAtNightDetector, SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE)
    generator.addScheduler(scheduler)
measurements = generator.generate(t_beginning, t_end)

I can also give you the code in java if you prefer

8

Thanks @lmirallie,

I’m not against the code in Java because Python can be a bit hard to debug.
I couldn’t use the FOV detector because I’m missing the parameters for:

Also, I used (lat, lon, alt) of Lausanne from wikipédia…

I did reproduce your bug but I don’t know where it comes from.
What I’ve noticed is that:

  • If the eclipse detector is removed, the problem disappears
  • If there’s only one satellite, the problem disappears. As soon as there are at least two satellites, then the bug comes back.

I also tried to move the creation of the detectors in the loop since I thought there would be problems using the same detector for all propagators. But it didn’t change the outcome.

Maxime

9

Ok thanks, thats very interesting. I indeed see that if I just consider one detector, the bug does not appear. In the following, I’m including the satellites TLEs from a file like the download from celestrak.

        final double latitude  = FastMath.toRadians(46.52183d);
        final double longitude = FastMath.toRadians(6.63270d);
        final double altitude  = 300.0d;
        final double azimuth  = 90.0d - FastMath.toRadians(0.0d);
        final double elevation = FastMath.toRadians(90.0d);
        final double height    = FastMath.toRadians(160.0d);
        final double width     = FastMath.toRadians(160.0d);
        
        
        final TimeScale utc = TimeScalesFactory.getUTC();
        Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, earthFrame);
        
        GeodeticPoint station = new GeodeticPoint(latitude, longitude, altitude);
        TopocentricFrame stationFrame = new TopocentricFrame(earth, station, "Lausanne");
        CelestialBody sun = CelestialBodyFactory.getSun();
        GroundStation groundStation = new GroundStation(stationFrame);
        
        Vector3D center = new Vector3D(FastMath.cos(azimuth) * FastMath.cos(elevation), FastMath.sin(azimuth) * FastMath.cos(elevation), FastMath.sin(elevation));
        Vector3D axis2 = Vector3D.crossProduct(Vector3D.PLUS_K, center).normalize();
        Vector3D axis1 = Vector3D.crossProduct(axis2, center).normalize();
        double ha1 = width/2.;
        double ha2 = height/2.;
        
        // FOV ground detector
        DoubleDihedraFieldOfView fov = new DoubleDihedraFieldOfView(center, axis1, ha1, axis2, ha2, 0.0);
        GroundFieldOfViewDetector fovGroundDetector = new GroundFieldOfViewDetector(stationFrame, fov);
        fovGroundDetector = fovGroundDetector.withHandler(new ContinueOnEvent());
        fovGroundDetector = fovGroundDetector.withMaxCheck(1.);
        
        // night detector (his g function return true when the ground (ie the station) is at night)
        GroundAtNightDetector nightDetector = new GroundAtNightDetector(stationFrame, sun, GroundAtNightDetector.CIVIL_DAWN_DUSK_ELEVATION, null);
        nightDetector = nightDetector.withHandler(new ContinueOnEvent());

        // eclipse detector ( to check if the satellite is in the earth penumbra, thus invisible by station)
        //EclipseDetector eclipseDetector = new EclipseDetector(sun, Constants.SUN_RADIUS, earth).withPenumbra().withHandler(new ContinueOnEvent());

        // combination of the two detectors for a boolean detector that return true when both are true
        BooleanDetector fovAtNightDetector = BooleanDetector.andCombine(BooleanDetector.notCombine(fovGroundDetector), nightDetector);
        fovAtNightDetector = fovAtNightDetector.withHandler(new ContinueOnEvent()).withMaxCheck(1.);
        
        final AbsoluteDate t_beginning = new AbsoluteDate(2023, 10, 10, 18, 22, 0.0, utc);
        final AbsoluteDate t_end       = new AbsoluteDate(2023, 10, 10, 23, 22, 0.0, utc);
        
        
        double step = 60.0d;
        double[] sigma_ = {1., 1.};
        double[] baseWeight_ = {1., 1.};

        Generator generator = new Generator();
        for (String[] str : satellites) 
        {
        TLE myTLE = new TLE(str[0], str[1]);
        TLEPropagator propagator = TLEPropagator.selectExtrapolator(myTLE);
        ObservableSatellite sat = generator.addPropagator(propagator);

        AngularAzElBuilder azElBuilder = new AngularAzElBuilder(null, groundStation, sigma_, baseWeight_, sat);
        EventBasedScheduler scheduler = new EventBasedScheduler(azElBuilder, new FixedStepSelector(step, utc),
                                                     generator.getPropagator(sat), nightDetector,
                                                     SignSemantic.FEASIBLE_MEASUREMENT_WHEN_POSITIVE);
        generator.addScheduler(scheduler);
        }
        generator.generate(t_beginning, t_end);
10

Hello @lmirallie,

It seems to work with your code on develop branch (i.e. future 12.0 version).

Could you try on your side and tell me if anything goes wrong?
Telescope_Detector_Lmirallie.java (8.2 KB)

I couldn’t tell you exactly why it’s working with 12 and not with 11.3… Maybe @luc guess on issue 1164 was right after all.

Cheers,
Maxime

11

Thank a lot! It does actually seems to work on my side as well, thanks for your help!
But only for less than 10 satellites. With more, the error comes back!
Best,
Louis