irfpy.util.intersection

Calculate intersection between LoS and sphere.

Code author: Yoshifumi Futaana

Line of sight and sphere

Given a observation point, P, and the looking vector, V, whether it crosses a sphere with radius of R at the point C.

Use los_sphere() function. To get limb position, use los_limb_sphere().

irfpy.util.intersection.los_sphere(pos, view, cent, radius, insphere='exception')

Returns the intersection of the LOS at the sphere surface.

          * pos
      view \
            \ /---\
intersection *     \
            /       \
            |   *   |
            \ cent  /
             \     /
              \---/ sphere with radius

Parameters:

• pos (numpy.array) – Position of observer

• view (numpy.array) – Viewing direction

• cent (numpy.array) – Center of the target sphere

• radius (float) – Radius of the sphere

• insphere – If the observer position is in sphere, how to handle? “exception” will raise RuntimeError, “none” will return None.

Returns:

Position of the intersection as vector3d.Vector3d instance. None will be returned if no intersection is found.

Return type:

vector3d.Vector3d instance or None.

>>> print(los_sphere([10,0,0], [-1,0,0], [-5,0,0], 3))
Vector3d( -2, 0, 0 )

>>> print(los_sphere([10,0,0], [0,1,0], [-5,0,0], 3))
None

>>> print(los_sphere([10, 0, 0], [1, 0, 0], [-5, 0, 0], 3))
None

If the given pos is in the sphere, RuntimeError will be raised as default.

>>> try:
...    los_sphere([10,0,0], [-1,0,0], [-5,0,0], 30)
...    print('RuntimeError should have been raised.')
... except RuntimeError as e:
...    print('OK')
... except:
...    print('Not correct exception.')
OK

irfpy.util.intersection.intersections_sphere(pos, view, cent, radius)

Returnt the intersection coordinates.

          * pos
      view \
            \ /----\
intersection1*      \
            / \ cent \
            |  \ *   |
            \   \    /
             \   \  /
              \--*-/ sphere with radius
                 intersections2

Parameters:

• pos – (3,) or (N, 3) shaped array of positions

• view – (3,) or (N, 3) shaped array of viweing directions (or velocity vector)

• cent – The center of the reference sphere

• radius – The radius of the reference spehre

Returns:

A tuple of (is1, is2, k). is1 and is2 is array with shape of (3,) or (N, 3). k is array with shape of (2,) or (N,2)

Samples. Let’s consider the sphere cetnered at (1, 1, 0) with a radius of 1.

>>> cent = [1, 1, 0]
>>> r = 1

If one places at (-1, 1, 0), looking toward (1, 0, 0), then one sees the interesections at (0, 1, 0) and (2, 1, 0).

>>> pos = [-1, 1, 0]
>>> view = [1, 0, 0]
>>> is1, is2, k = intersections_sphere(pos, view, cent, r)
>>> print(is1)
[0. 1. 0.]
>>> print(is2)
[2. 1. 0.]

Let’s see the next sample. For the same shere as above, but let’s see from the different position, (-1, 2, 0). The view vector is assumed to be (2, 0, 0). In this geometry, the intersection is 1 point, (1, 2, 0). The k-values are 1 (remmeber again the k-values are normalized by the view vector, in this case (2, 0, 0).).

>>> pos = [-1, 2, 0]
>>> view = [2, 0, 0]
>>> is3, is4, k = intersections_sphere(pos, view, cent, r)
>>> print(is3)
[1. 2. 0.]
>>> print(is4)
[1. 2. 0.]

What if one cannot see the sphere. Set the position (3, -1, 0) and looking toward (0, 1, 0). nan will be returned, while some warning messages are shown. >>> pos = [3, -1, 0] >>> view = [0, 1, 0] >>> p = intersections_sphere(pos, view, cent, r) >>> print(p) (array([nan, nan, nan]), array([nan, nan, nan]), array([nan, nan]))

The last example is if the position is inside the sphere, e.g. at (1.5, 1, 0), looking toward (0, 1, 0) with a length of square root of 3. In this case, k- becomes negative. But still you can calculate the intersection points.

>>> pos = [1.5, 1, 0]
>>> view = [0, np.sqrt(3), 0]
>>> is5, is6, k = intersections_sphere(pos, view, cent, r)
>>> print(is5)
[1.5       0.1339746 0.       ]
>>> print(is6)
[1.5       1.8660254 0.       ]

Vectorized sample

The above example can be vectorized, which has much higher performance.

>>> pos = [[-1, 1, 0], [-1, 2, 0,], [3, -1, 0], [1.5, 1, 0]]      # (N=4, 3) shape
>>> view = [[1, 0, 0], [2, 0, 0], [0, 1, 0], [0, np.sqrt(3), 0]]  # (N=4, 3) shape
>>> ism, isp, k = intersections_sphere(pos, view, cent, r)
>>> print(ism.shape)    # (N=4, 3) shape
(4, 3)
>>> print(ism)
[[0.        1.        0.       ]
 [1.        2.        0.       ]
 [      nan       nan       nan]
 [1.5       0.1339746 0.       ]]
>>> print(isp)    # (4, 3) shape
[[2.        1.        0.       ]
 [1.        2.        0.       ]
 [      nan       nan       nan]
 [1.5       1.8660254 0.       ]]

irfpy.util.intersection.lengths_to_intersections(pos, view, cent, radius)

Return the lengths to the intersections.

          * pos
      view \
            \ /----\
intersection1*      \
            / \ cent \
            |  \ *   |
            \   \    /
             \   \  /
              \--*-/ sphere with radius
                 intersections2

Parameters:

• pos – (3,) or (N, 3) shaped array of positions

• view – (3,) or (N, 3) shaped array of viweing directions (or velocity vector)

• cent – The center of the reference sphere

• radius – The radius of the reference spehre

Returns:

An array of (2,) or (N, 2), where 2 refers to the length of intersection1 and intersection2 normalized by the view. Always small for the first, large (or equal) for the second. If no intersection, (np.nan, np.nan) is returned.

The returned is two floating point. To calculate the interesection1 and 2, you may use intersections_sphere() function.

Samples. Let’s consider the sphere cetnered at (1, 1, 0) with a radius of 1.

>>> cent = [1, 1, 0]
>>> r = 1

If one places at (-1, 1, 0), looking toward (1, 0, 0), then one sees the interesections at (0, 1, 0) and (2, 1, 0); then k- = 1 and k+ =3. Note that the k-values are normalized by the length of view vector.

>>> pos = [-1, 1, 0]
>>> view = [1, 0, 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
1.0 3.0

Let’s see the next sample. For the same shere as above, but let’s see from the different position, (-1, 2, 0). The view vector is assumed to be (2, 0, 0). In this geometry, the intersection is 1 point, (1, 2, 0). The k-values are 1 (remmeber again the k-values are normalized by the view vector, in this case (2, 0, 0).).

>>> pos = [-1, 2, 0]
>>> view = [2, 0, 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
1.0 1.0

What if one cannot see the sphere. Set the position (3, -1, 0) and looking toward (0, 1, 0). nan will be returned, while some warning messages are shown. >>> pos = [3, -1, 0] >>> view = [0, 1, 0] >>> km, kp = lengths_to_intersections(pos, view, cent, r)

The last example is if the position is inside the sphere, e.g. at (1.5, 1, 0), looking toward (0, 1, 0) with a length of square root of 3. In this case, k- becomes negative. Still you can calculate the intersection points as above.

>>> pos = [1.5, 1, 0]
>>> view = [0, np.sqrt(3), 0]
>>> km, kp = lengths_to_intersections(pos, view, cent, r)
>>> print(km, kp)
-0.5 0.5

Vectorized sample

The above example can be vectorized, which has much higher performance.

>>> pos = [[-1, 1, 0], [-1, 2, 0,], [3, -1, 0], [1.5, 1, 0]]
>>> view = [[1, 0, 0], [2, 0, 0], [0, 1, 0], [0, np.sqrt(3), 0]]
>>> k = lengths_to_intersections(pos, view, cent, r)
>>> print(k[:, 0])
[ 1.   1.   nan -0.5]
>>> print(k[:, 1])
[3.  1.  nan 0.5]

irfpy.util.intersection.los_limb_sphere(pos, view, cent, radius)

Return the intersection of LOS or limb.

This function returns the intersection between view and the sphere. If the configuration results in the intersection, the result is the same as los_sphere().

But if the intersection is not found (view direction is into space not sphere), this function returns the limb position of the sphere in the same plane defined by the vectors of pos and view.

             * pos
             |  view
  /---\      |
 /     \     |
/       \    |
|   *   *----|  No intersection, returning the neighboring surface
\ cent  /    |
 \     /
  \---/ sphere with radius

Parameters:

• pos (numpy.array) – Position of observer

• view (numpy.array) – Viewing direction

• cent (numpy.array) – Center of the target sphere

• radius (float) – Radius of the sphere

Returns:

Position of the intersection or limb as vector3d.Vector3d instance.

Return type:

vector3d.Vector3d instance or None.

>>> print(los_limb_sphere([10,0,0], [-1,0,0], [-5,0,0], 3))
Vector3d( -2, 0, 0 )

>>> print(los_limb_sphere([10,0,0], [0,1,0], [-5,0,0], 3))
Vector3d( -4.4, 2.93939, 0 )

>>> print(los_limb_sphere([10,0,0], [-1,1,0], [-5,0,0], 3))
Vector3d( -4.4, 2.93939, 0 )

>>> print(los_limb_sphere([15,0,0], [-1,1,0], [0,0,0], 3))
Vector3d( 0.6, 2.93939, 0 )

Note for developer

The los_sphere is developed in vector3d.Vector3d class, but los_limb_sphere is written basically using np.array class.

irfpy.util.intersection.isvisible(observer_xyz, target_xyz, sphere_center_xyz, radius)

Return if the target is visible from the observer if there is a sphere.

Parameters:

• observer_xyz – Numpy array of observer location. Must be (3,) shape.

• target_xyz – Numpy array of target location. Must be (3,) shape.

• sphere_center_xyz – Numpy array of sphere center. Must be (3,) shape.

• radius – Radius.

Returns:

True if the target is visible.

Use case

It is commonly used if a spacecraft (observer) can see specific location (target) if a planetary body exist.

Theory

The observer can only see the specific target either when

• No intersection of its view angle and the sphere

• Within a virtual sphere centered at the observer with radius of the distance between the observer and the limb of the sphere

• But not inside the sphere

Sample

Suppose a sphere at (-5, 0, 0) with radius of 3. Observer is assumed to be at (5, 0, 0) in this example.

Point (-1.5, 0, 0) should be visible from the observer at (5, 0, 0).

>>> print(isvisible([5, 0, 0], [-1.5, 0, 0,], [-5., 0, 0], 3))
True

Point (-2.5, 0, 0) should not be visible from any observer because the point is inside the sphere

>>> print(isvisible([5, 0, 0], [-2.5, 0, 0], [-5, 0, 0], 3))
False

Point (10, 0, 0) should be visible, since the sphere is opposite direction for the observer.

>>> print(isvisible([5, 0, 0], [10, 0, 0], [-5, 0, 0], 3))
True

Point (-5, 3.001, 0) should not be visible, because it is below the horizon

>>> print(isvisible([5, 0, 0], [-5, 3.001, 0], [-5, 0, 0], 3))
False

Point (-10, 0, 0) should not be visible, because it is behind the sphere.

>>> print(isvisible([5, 0, 0], [-10, 0, 0], [-5, 0, 0], 3))
False

irfpy.util.intersection.arevisible(observer_xyz, targets_xyz, sphere_center_xyz, radius)

Return if the target is visible from the observer if there is a sphere.

Parameters:

• observer_xyz – Numpy array of observer location. Must be (3,) shape.

• targets_xyz – Numpy array of target location. Must be (N, 3) shape.

• sphere_center_xyz – Numpy array of sphere center. Must be (3,) shape.

• radius – Radius.

Returns:

True if the target is visible.

This is for evaluating multiple particles of visibility, i.e., parallelized version of isvisible().

Sample

Suppose a sphere at (-5, 0, 0) with radius of 3. Observer is assumed to be at (5, 0, 0) in this example.

• Point (-1.5, 0, 0) should be visible from the observer at (5, 0, 0).

• Point (-2.5, 0, 0) should not be visible from any observer because the point is inside the sphere

• Point (10, 0, 0) should be visible, since the sphere is opposite direction for the observer.

• Point (-5, 3.001, 0) should not be visible, because it is below the horizon

• Point (-10, 0, 0) should not be visible, because it is behind the sphere.

>>> print(isvisible([5, 0, 0], [-1.5, 0, 0,], [-5., 0, 0], 3))
True
>>> print(isvisible([5, 0, 0], [-2.5, 0, 0], [-5, 0, 0], 3))
False
>>> print(isvisible([5, 0, 0], [10, 0, 0], [-5, 0, 0], 3))
True
>>> print(isvisible([5, 0, 0], [-5, 3.001, 0], [-5, 0, 0], 3))
False
>>> print(isvisible([5, 0, 0], [-10, 0, 0], [-5, 0, 0], 3))
False

>>> xyz = np.array([[-1.5, 0, 0], [-2.5, 0, 0], [10., 0, 0, ], [-5, 3.001, 0], [-10, 0, 0]])

>>> print(arevisible([5, 0, 0], xyz, [-5, 0, 0], 3))
[ True False  True False False]

irfpy.util.intersection.sphere_limb_plane_projection(sphere_location, radius, resolution=360)

Return a (theta, phi) list of a sphere projected on a plane.

Parameters:

• sphere_location – (x, y, z) of the sphere.

• radius – The radius of the sphere

Returns:

A pair (theta, phi). theta (and phi) is a numpy array with a size of (360,), or specified by the parameter resolution. These values are with the unit of radians. theta ranges [0, pi], and phi ranges [-pi, pi].

Raises:

ValueError – if the sphere radius is larger than the distance to the sphere. (i.e., the observer at origin is inside the sphere)

It is used for example, to make a all sky map for ions or ENAs, but also to plot the planet.

>>> theta, phi = sphere_limb_plane_projection([300, 0, 1000], 300)
>>> print(phi.shape)
(360,)

This gives the angles (theta and phi) of the planet at (300, 0, 1000) [km] with radius of 300 [km]. You get 360 points by default to be connected (or scatter plotted).

The returned angles are in radians. The 100’th sampling point is shown as

>>> print('{:.4f}'.format(phi[100]))
0.8962
>>> print('{:.4f}'.format(theta[100]))
0.3705

You can of course plot it

>>> plt.plot(phi, theta, '.')    

Theory

It is a combination of frame rotations.

Sphere loated at (X, Y, Z) with radius R.

Observer at the origin.

Then, the projection to all sky map (theta, phi) is represented by

theta = acos(z)

phi = atan2(y, x)

(x, y, z) = ((cos PHI, -sin PHI, 0), (sin PHI, cos PHI, 0), (0, 0, 1)) . ((cos THETA, 0, sin THETA), (0, 1, 0), (-sin THETA, 0, cos THETA)) . (sin theta0 cos alpha, sin theta0 sin alpha, cos theta0)

PHI = atan2(Y, X)

THETA = acos(Z / L)

L = sqrt(X^2 + Y^2 + Z^2)

theta0 = asin(R / L)

alpha is the parameter to range from -pi to pi.