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')[source]¶
Returns the intersection of the LOS at the sphere surface.
* pos view \ \ /---\ intersection * \ / \ | * | \ cent / \ / \---/ sphere with radius
- Parameters:
pos (
numpy.array
) – Position of observerview (
numpy.array
) – Viewing directioncent (
numpy.array
) – Center of the target sphereradius (
float
) – Radius of the sphereinsphere – 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 orNone
.
>>> 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)[source]¶
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)[source]¶
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)
; thenk- = 1
andk+ =3
. Note that thek
-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)[source]¶
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 aslos_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
andview
.* pos | view /---\ | / \ | / \ | | * *----| No intersection, returning the neighboring surface \ cent / | \ / \---/ sphere with radius
- Parameters:
pos (
numpy.array
) – Position of observerview (
numpy.array
) – Viewing directioncent (
numpy.array
) – Center of the target sphereradius (
float
) – Radius of the sphere
- Returns:
Position of the intersection or limb as
vector3d.Vector3d
instance.- Return type:
vector3d.Vector3d
instance orNone
.
>>> 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, butlos_limb_sphere
is written basically using np.array class.
- irfpy.util.intersection.isvisible(observer_xyz, target_xyz, sphere_center_xyz, radius)[source]¶
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)[source]¶
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)[source]¶
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.