# 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().

Returns the intersection of the LOS at the sphere surface.

    * pos
view                       \ /---          intersection *                           /                             |   *   |
\ cent  /
\     /

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


Returnt the intersection coordinates.

    * pos
view                       \ /----          intersection1*                            / \ cent                       |  \ *   |
\   \    /
\   \  /
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

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.       ]]


Return the lengths to the intersections.

    * pos
view                       \ /----          intersection1*                            / \ cent                       |  \ *   |
\   \    /
\   \  /
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

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]


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  /    |
\     /

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.

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.

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


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.

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]


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

Parameters
• sphere_location – (x, y, z) 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.