Table of contents
CSPICE_SURFPV finds the state (position and velocity) of the surface
intercept defined by a specified ray, ray velocity, and ellipsoid.
Given:
stvrtx the state of a ray's vertex.
help, stvrtx
DOUBLE = Array[6]
The first three components of `stvrtx' are the vertex's x, y,
and z position components; the vertex's x, y, and z velocity
components follow.
The reference frame relative to which `stvrtx' is
specified has axes aligned with with those of a
triaxial ellipsoid. See the description below of
the arguments `a', `b', and `c'.
The vertex may be inside or outside of this
ellipsoid, but not on it, since the surface
intercept is a discontinuous function at
vertices on the ellipsoid's surface.
No assumption is made about the units of length
and time, but these units must be consistent with
those of the other inputs.
stdir the state of the input ray's direction vector.
help, stdir
DOUBLE = Array[6]
The first three components of `stdir' are a non-zero vector
giving the x, y, and z components of the ray's direction; the
direction vector's x, y, and z velocity components follow.
`stdir' is specified relative to the same reference
frame as is `stvrtx'.
a,
b,
c respectively, the lengths of a triaxial ellipsoid's semi-axes
lying along the x, y, and z axes of the reference frame relative
to which `stvrtx' and `stdir' are specified.
help, a
DOUBLE = Scalar
help, b
DOUBLE = Scalar
help, c
DOUBLE = Scalar
the call:
cspice_surfpv, stvrtx, stdir, a, b, c, stx, found
returns:
stx the state of the intercept of the input ray on the surface of
the input ellipsoid.
help, stx
DOUBLE = Array[6]
The first three components of `stx' are the intercept's x, y,
and z position components; the intercept's x, y, and z velocity
components follow.
`stx' is specified relative to the same reference
frame as are `stvrtx' and `stdir'.
`stx' is defined if and only if both the intercept
and its velocity are computable, as indicated by
the output argument `found'.
The position units of `stx' are the same as those of
`stvrtx', `stdir', and `a', `b', and `c'. The time units are
the same as those of `stvrtx' and `stdir'.
found a logical flag indicating whether `stx' is defined.
help, found
BOOLEAN = Scalar
`found' is True if and only if both the intercept and its
velocity are computable. Note that in some cases the intercept
may computable while the velocity is not; this can happen for
near-tangency cases.
None.
Any numerical results shown for this example may differ between
platforms as the results depend on the SPICE kernels used as input
and the machine specific arithmetic implementation.
1) Illustrate the role of the ray vertex velocity and
ray direction vector velocity via several simple cases. Also
show the results of a near-tangency computation.
Example code begins here.
PRO surfpv_ex1
F1 = '(A,3E20.12)'
a = 1.D0
b = 2.D0
c = 3.D0
print, ' '
print, 'Ellipsoid radii:'
print, ' A = ', a
print, ' B = ', b
print, ' C = ', c
print, ' '
print, 'Case 1: Vertex varies, direction is constant'
print, ' '
stvrtx = [ 2.D0, 0.D0, 0.D0, 0.D0, 0.D0, 3.D0 ]
stdir = [ -1.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0 ]
print, 'Vertex:'
print, format=F1, ' ', stvrtx[0:2]
print, 'Vertex velocity:'
print, format=F1, ' ', stvrtx[3:5]
print, 'Direction:'
print, format=F1, ' ', stdir[0:2]
print, 'Direction velocity:'
print, format=F1, ' ', stdir[3:5]
cspice_surfpv, stvrtx, stdir, a, b, c, stx, found
if ( not found ) then begin
print, ' No intercept state found.'
endif else begin
print, 'Intercept:'
print, format=F1, ' ', stx[0:2]
print, 'Intercept velocity:'
print, format=F1, ' ', stx[3:5]
print, ' '
endelse
print, ' '
print, 'Case 2: Vertex and direction both vary'
print, ' '
stdir[5] = 4.D0
print, 'Vertex:'
print, format=F1, ' ', stvrtx[0:2]
print, 'Vertex velocity:'
print, format=F1, ' ', stvrtx[3:5]
print, 'Direction:'
print, format=F1, ' ', stdir[0:2]
print, 'Direction velocity:'
print, format=F1, ' ', stdir[3:5]
cspice_surfpv, stvrtx, stdir, a, b, c, stx, found
if ( not found ) then begin
print, ' No intercept state found.'
endif else begin
print, 'Intercept:'
print, format=F1, ' ', stx[0:2]
print, 'Intercept velocity:'
print, format=F1, ' ', stx[3:5]
print, ' '
endelse
print, ' '
print, 'Case 3: Vertex and direction both vary;'
print, ' near-tangent case.'
print, ' '
stvrtx[2] = c - 1.D-15
stvrtx[5] = 1.D299
stdir[5] = 1.D299
print, 'Vertex:'
print, format=F1, ' ', stvrtx[0:2]
print, 'Vertex velocity:'
print, format=F1, ' ', stvrtx[3:5]
print, 'Direction:'
print, format=F1, ' ', stdir[0:2]
print, 'Direction velocity:'
print, format=F1, ' ', stdir[3:5]
cspice_surfpv, stvrtx, stdir, a, b, c, stx, found
if ( not found ) then begin
print, ' No intercept state found.'
endif else begin
print, 'Intercept:'
print, format=F1, ' ', stx[0:2]
print, 'Intercept velocity:'
print, format=F1, ' ', stx[3:5]
print, ' '
endelse
END
When this program was executed on a Mac/Intel/IDL8.x/64-bit
platform, the output was:
Ellipsoid radii:
A = 1.0000000
B = 2.0000000
C = 3.0000000
Case 1: Vertex varies, direction is constant
Vertex:
2.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Vertex velocity:
0.000000000000E+00 0.000000000000E+00 3.000000000000E+00
Direction:
-1.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Direction velocity:
0.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Intercept:
1.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Intercept velocity:
0.000000000000E+00 0.000000000000E+00 3.000000000000E+00
Case 2: Vertex and direction both vary
Vertex:
2.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Vertex velocity:
0.000000000000E+00 0.000000000000E+00 3.000000000000E+00
Direction:
-1.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Direction velocity:
0.000000000000E+00 0.000000000000E+00 4.000000000000E+00
Intercept:
1.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Intercept velocity:
0.000000000000E+00 0.000000000000E+00 7.000000000000E+00
Case 3: Vertex and direction both vary;
near-tangent case.
Vertex:
2.000000000000E+00 0.000000000000E+00 3.000000000000E+00
Vertex velocity:
0.000000000000E+00 0.000000000000E+00 1.000000000000E+299
Direction:
-1.000000000000E+00 0.000000000000E+00 0.000000000000E+00
Direction velocity:
0.000000000000E+00 0.000000000000E+00 1.000000000000E+299
Intercept:
2.580956827952E-08 0.000000000000E+00 3.000000000000E+00
Intercept velocity:
-3.874532036208E+306 0.000000000000E+00 2.999999974190E+299
The position and velocity of the ray's vertex as well as the
ray's direction vector and velocity vary with time. The
inputs to cspice_surfpv may be considered the values of these
vector functions at a particular time, say `t0'. Thus
State of vertex: stvrtx = ( V(t0), V'(t0) )
State of direction vector: stdir = ( D(t0), D'(t0) )
To determine the intercept point, W(t0), we simply compute the
intersection of the ray originating at V(t0) in the direction of
D(t0) with the ellipsoid
2 2 2
x y z
---- + ---- + ---- = 1
2 2 2
a b c
W(t) is the path of the intercept point along the surface of
the ellipsoid. To determine the velocity of the intercept point,
we need to take the time derivative of W(t), and evaluate it at
t0. Unfortunately W(t) is a complicated expression, and its
derivative is even more complicated.
However, we know that the derivative of W(t) at `t0', W'(t0), is
tangent to W(t) at `t0'. Thus W'(t0) lies in the plane that is
tangent to the ellipsoid at t0. Let X(t) be the curve in the
tangent plane that represents the intersection of the ray
emanating from V(t0) with direction D(t0) with that tangent
plane.
X'(t0) = W'(t0)
The expression for X'(t) is much simpler than that of W'(t);
cspice_surfpv evaluates X'(t) at t0.
Derivation of X(t) and X'(t)
----------------------------
W(t0) is the intercept point. Let `n' be a surface normal at I(t0).
Then the tangent plane at W(t0) is the set of points X(t) such
that
< X(t) - I(t0), N > = 0
X(t) can be expressed as the vector sum of the vertex
and some scalar multiple of the direction vector,
X(t) = V(t) + s(t) * D(t)
where s(t) is a scalar function of time. The derivative of
X(t) is given by
X'(t) = V'(t) + s(t) * D'(t) + s'(t) * D(t)
We have V(t0), V'(t0), D(t0), D'(t0), W(t0), and `n', but to
evaluate X'(t0), we need S(t0) and S'(t0). We derive an
expression for s(t) as follows.
Because X(t) is in the tangent plane, it must satisfy
< X(t) - W(t0), N > = 0.
Substituting the expression for X(t) into the equation above
gives
< V(t) + s(t) * D(t) - W(t0), N > = 0.
Thus
< V(t) - W(t0), N > + s(t) * < D(t), N > = 0,
and
< V(t) - W(t0), N >
s(t) = - ---------------------
< D(t), N >
The derivative of s(t) is given by
s'(t) =
< D(t),N > * < V'(t),N > - < V(t)-W(t0),N > * < D'(t),N >
- -------------------------------------------------------------
2
< D(t), N >
1) If the input ray's direction vector is the zero vector, an
error is signaled by a routine in the call tree of this
routine.
2) If any of the ellipsoid's axis lengths is nonpositive, an
error is signaled by a routine in the call tree of this
routine.
3) If the vertex of the ray is on the ellipsoid, the error
SPICE(INVALIDVERTEX) is signaled by a routine in the call tree
of this routine.
4) If any of the input arguments, `stvrtx', `stdir', `a', `b' or
`c', is undefined, an error is signaled by the IDL error
handling system.
5) If any of the input arguments, `stvrtx', `stdir', `a', `b' or
`c', is not of the expected type, or it does not have the
expected dimensions and size, an error is signaled by the Icy
interface.
6) If any of the output arguments, `stx' or `found', is not a
named variable, an error is signaled by the Icy interface.
None.
None.
ICY.REQ
None.
J. Diaz del Rio (ODC Space)
-Icy Version 1.0.0, 09-AUG-2021 (JDR)
ellipsoid surface point and velocity
|