Re: [Orekit Users] Question on using getIntersectionPoint

[Luc Maisonobe <Luc.Maisonobe@c-s.fr>]

Le 17/04/2012 10:32, Petrus Hyvönen a écrit :
> Hi Luc,
> 
> Thanks for your response.
> 
> Some clarifications,
> both my points, ref_pos and geo_pos are satellites, where I would like
> to know where I have a clear visibility between them. To do so I thought
> to check if there is an Earth that the line between them intersects.

Oh, I see. then it is really the correct way to use Line.

The null result is one case when you know you have a clear path, since
the Earth is really out of the way. However, there are also cases when
the line (which in infinite in both directions) does intersect Earth but
is farther than the LEO satellite (for example when LEO is just in the
nadir direction as seen from GEO. So if the returned point is non-null,
you have to check if the LEO is between the GEO and the intersection
point (in which case you do have a clear path) or if the LEO is past the
Earth. One way to do this is to compare the abscissae of all points (the
line is oriented and has a canonical zero at the point where is is
closest to origin). Try something along these lines (Java syntax here,
not Python):

  Vector3D intersection =
    earth.getIntersectionPoint(los, geoPos, inertialFrame, time);
  if (intersection == null) {
    return CLEAR_PATH;
  }

  double xIntersection = los.getAbscissa(intersection);
  double xLeo          = los.getAbscissa(refPos);
  double xGeo          = los.getAbscissa(geoPos);
  if ((xLeo - xIntersection) * (xLeo - xGeo) < 0) {
    // leo is in between, we have a clear path
    return CLEAR_PATH;
  } else {
    // leo is past Earth, so it is behind Earth limb
    return OBSTRUCTION;
  }

Be aware that you may have atmospheric effects on the path between the
two satellites, so you should perhaps increase radius to take at least
part of the atmosphere into account.

> 
> See below,
> 
> On Tue, Apr 17, 2012 at 10:06 AM, Luc Maisonobe <Luc.Maisonobe@c-s.fr
> <mailto:Luc.Maisonobe@c-s.fr>> wrote:
> 
>     Le 16/04/2012 22:59, Petrus Hyvönen a écrit :
>     > Hi,
> 
>     Hi Petrus,
> 
>     >
>     > I'm having a question if this is the right way to use
>     getIntersectionPoint.
>     > I'm trying to simulate the coverage of a relay satellite in
>     > geostationary orbit.
>     >
>     > The geostationary satellite I have modelled as a GroundPoint with 6
>     > earth radius altitude, and then transformed to an inertialFrame
>     EME2000
>     > frame.
>     >
>     > I have my two positions as Vector3D coordinates in inertialFrame.
>     > geo_pos
>     >
>     > ref_pos
> 
> 
>     I assume ref_pos is a point on ground, is this right ? If not, see the
>     last remark at the end of this message.
> 
>     With such a model, you can use an Earth frame for both points instead of
>     an inertial frame. It does not really matter, though as conversions will
>     be done under the hood as needed.
> 
>     >
>     >
>     > Create a Line with
>     > los = Line(ref_pos, geo_pos)
>     >
>     > ip = mysc.ore_earth.getIntersectionPoint(los, geo_pos,
>     inertialFrame, time)
>     >
>     >
>     > and if this is == None (Null) then it has a clear view I would think.
> 
>     No. If the result is null, then the line that joins your two points does
>     not intersect the ore_earth BodyShape. This should not happen if your
>     ref point is on Earth (except perhaps at limb as seen from geo_pos).
> 
>     If the result is non null, then the line enters Earth surface at the
>     returned point (and exits at some other point further away).
> 
>     >
>     >
>     > The results are not what I expect them to be, no contacts in what I
>     > would expect the coverage region and some at other locations.
>     Would this
>     > be the/a correct way to use getIntersectionPoint?
> 
>     It is a strange way, but it should work.
>     As you already have a ground point ref_pos, a simpler way could be to
>     build a TopocentricFrame from it and to use its getElevation method to
>     check is it is above observation threshold.
> 
> 
> As a temporary solution I have used TopocentricFrame, projected the LEO
> satellite footprint on earth and used a slightly negative minimum
> elevation of the GEO satellite as critera in the subpoint topocentric
> frame. Not pretty but seems to work.
>  
> 
> 
>     > Is it ok that the Line
>     > is in a different frame than the OneAxisEllipsoid, as long as it is
>     > specified in the call?
> 
>     Yes, it's OK. Both los and geo_pos must be in the same specified frame,
>     but this frame can be anything. These elements are converted to the
>     BodyShape internal frame at the start of the computation.
> 
>     >
>     >
>     > I do not understand the Line object really, los.contains(geo_pos)
>     gives
>     > False, as well as los.contains(ref_pos). The los.origin gives a
>     > different value than given in the constructor?
> 
>     This is strange. Could you send some numerical values with coordinates
>     of the various points in all the involved frames (body shape and
>     inertial) ? There was also a Line class up to version 5.x of the Orekit
>     library, but now we use directly the class from Apache Commons Math
>     (which has different and incompatible meanings in the construction
>     parameters). What version do you use, Apache Commons Math 2.x/Orekit 5.x
>     or Apache commons Math 3.0/Orekit development ? Note that you cannot mix
>     these versions.
> 
> 
> I am using common math 3.0 downloaded using mvn and latest orekit 6.0
> SNAPSHOT. I am using this under Python, but I do not think/hope this is
> influencing the results, all calculations are done in the jvm.
> 
> I am not sure if the intersection thing is related to this, but it seems
> like the distance threshold for Line is preset at 1e-10, which seems to
> be lower than my practical case:
> 
> In [84]: p1 = Vector3D(163919.77, -6224471.67, -3350987.6)
> In [85]: p2 = Vector3D(18680140.69, -40551236.16, -23192.4)
> In [86]: l = Line(p1, p2)
> 
> In [87]: l.contains(p1)
> Out[87]: False
> 
> In [88]: l.distance(p1)
> Out[88]: 5.45683966183878e-10

OK. Then it is simply a numerical accuracy problem. The Line object in
Apache Commons Math do not simply retain the two points but compute some
internal parameters and do some normalization. I think the result is OK.
You can still use l.distance(p) < myOwnThreshold instead of the buil-in
contains method.

best regards,

Luc


> 
> But:
> In [110]: p1=Vector3D(0.0,0.0,0.0)
> In [111]: p2=Vector3D(0.0, 0.0, 3.0)
> In [112]: l=Line(p1,p2)
> 
> In [113]: l.contains(p1)
> Out[113]: True
> 
> In [114]: l.distance(p1)
> Out[114]: 0.0
>  
> Best Regards,
> /Petrus
>