Mars Odyssey/GRS Geometry

The purpose of this document is to summarize the geometrical considerations of the Odyssey spacecraft in orbit around Mars. All necessary quantities needed to undertake GRS data reductions are derived from first principles. Upon completion of these derivations, a minimum-list of required parameters is compiled. These parameters must be supplied, for each spectrum, by the appropriate spice calls.

Special care has been taken to cover the absolutely general case of a spacecraft in a non-circular orbit around an aspherical planet. We rely on the spacecraft state vector (vector position + velocity) for each spectrum. This eliminates the need to do an elliptical fit to the spacecraft trajectory, and preserves the full accuracy of whatever is considered in assembling the spice kernel. Any perturbations JPL considers in generating the spice kernel therefore will be propagated into our solution for the geometry at the mid-time of each spectrum. (Boris Semenov at JPL has indicated that the latency between a ``predict'' and ``after-the-fact'' will be on the order of 36 hours.)

Our approach is a vector-based approach. The equations are general. Once you choose a coordinate system in which to work, these equations are valid so long as you stay consistent by using this system throughout. In general, we will opt to use a Mars-fixed, Mars-centered (x, y, z) system, as it is most expeditious for our purposes. This coordinate system will be referred to as the ``Areocentric'' system. It is centered on the planet's barycenter (center-of-mass). The x-axis goes through the equator at the prime meridian (through longitude = 0^°). The z-axis goes through the North pole, and the y-axis completes the triad of a right-handed, orthogonal system (through the equator at longitude +90^° longitude). We adopt the standard IAU convention of longitude positive to the east, ranging from 0 to 360 degrees.

We divide the planet into 0.5^° by 0.5^° cells. These cells are bounded by meridians of longitude and parallels of latitude. To good approximation, each cell is a trapezoid. On Mars, a minute of arc is very nearly 1 km (987.275 m), so each side of a cell measures just about 30 km.

Preliminaries: Vector Review

For any vector [v\vec], we can form a unit vector by normalizing by the magnitude,

|| ® v   || = Ö   vx2 + vy2 + vz2

^ v   = ® v || ® v   ||  .

The dot product [a\vec] ·[b\vec] is a scalar given by:

® a   · ® b   = || ® a   || || ® b   || cos q ,

where q is the enclosed angle of the two vectors.

Solid angle of the detector at a cell center:

At any time t, and independent of coordinate system, the vector from the center of the planet to the center of a cell is [r\vec]_C and from the center of the planet the vector to the spacecraft is [R\vec]_S/C. (Refer to Figure 1.) The vector from a cell to the spacecraft is given by [d\vec]:

® d   = ® R   S/C  - ® r   C   .

The solid angle of the detector as viewed from the center of a cell ( º B_i,j,s) is given by

W =  A ||d||2  ,

where A is the projected cross-sectional area of the detector.

Note: We have characterized the spatial response of the detector in terms of the on-axis response R(0, 0). The resulting spatial maps therefore account for the variation in projected detector cross-section. Thus, we need not worry about how the cross-section varies with look angle. The appropriate cross-section to use is always that along the symmetry axis of the crystal, A ~ 34 cm² = 3.4 ×10^-9 km². As a sanity check, at the nadir,

Wnadir ~  3.4 ×10-9 4002 ~ 2.13 ×10-14 sr

Angle from a cell vertical to the spacecraft:

In essense, what we want here is the zenith distance to the spacecraft. (Refer to Figure 2.) For a spherical planet, the cell vertical (zenith) is just the normalized planetary radius vector to that point,

^ n   = ® r   C  || ® r   C  ||  .

This is not true for the case of an oblate spheroidal planet, nor is it true in general.

Similarly, the unit vector from the cell to the spacecraft is given by

^ d   S/C  = ® d   S/C  || ® d   S/C  ||  .

To determine the angle from cell vertical to spacecraft, merely take the dot product [^n] ·[^d]_S/C:

cosz = ® r   C  · ® d   S/C  || ® r   C  ||  || ® d   S/C  ||

An aside: For an oblate spheroidal planet, the zenith direction is given by a unit vector (l, b¢), where l is the longitude and b¢ is the areodetic latitude. In other words, the cell normal is given by the vector

® n   Cell  = é ê ê ê ê ê ë nx ny nz ù ú ú ú ú ú û = é ê ê ê ê ê ë cosb¢cosl cosb¢sin l sin b¢ ù ú ú ú ú ú û

Areodetic latitude relates to areocentric latitude b by:

b¢ = b- m sin 2b+  1 2 m2 sin2b- ¼    ,

where

m = e2/(2-e2) = f +  1 2 f2 -  1 4 f4 + ¼

and f is the flattening of the planet,

f = (a-b)/a = e2/2 + e4/8 + e6/16 + ¼    .

For Mars, f is 0.00647630 = 1/154.409, so the first two terms are all we need. Thus, m = 0.006497271 radians, or 0.37227 degree.

Also be careful here: although areodetic and areographic longitude are essentially the same, areodetic longitude usually is measured as positive west, instead of east, in which case the n_y component of [n\vec]_Cell would be -cosb¢sinl.

Calculate the atmospheric thickness along the line from the center of the cell to the spacecraft:

The airmass through which g-rays must travel to get to the spacecraft is given by sec z for a plane parallel atmosphere. Airmass is defined as the path length through the atmosphere, expressed in units of the vertical path:

X = sec z

(plane parallel)

For a spherical atmosphere, we can calculate the true airmass as a Taylor series expansion in (sec z - 1):

X = sec z - 0.018167 (sec z - 1)2 -0.002875 (sec z - 1)3-0.0008083 (sec z - 1)4  .

(spherical)

One can see that the correction at 60^° off-zenith is only on the order of 1%. Atmospheric attenuation of the gamma rays is a much stronger effect than this correction, so to good approximation we may take the airmass as just ~ sec z Barring topography, i.e. for a bowling ball planet, Ðz is the emission angle for the K array. To relate topography to atmospheric pressure, I quote Smith, et al. (2001) Mars Orbiter Laser Altimeter: experiment summary after the first year of global mapping of Mars. Jour. Geophs. Res. 106, 23689-23722. (This paper is required reading for an understanding of both topography and pressure equipotential surfaces for Mars.)

``The relationship between surface topography and atmospheric pressure provides a measure of the mass of the atmospheric column above a particular location, which is critical in the assessment of potential future landing sites... The average atmospheric pressure on Mars is ~6.1 mbars, which is close to the triple point of water. Early topographic models of Mars... were referenced to this atmospheric pressure surface. The use of a pressure surface as a reference introduced considerable error into estimates of elevation because of temporal variability in the height of the pressure surface due to seasonal variations in CO₂ content and dynamical motions in the atmosphere. Because of the importance of understanding surface pressure for models of circulation of the atmosphere and for landing-site assessment, MOLA elevations, which are referenced to a static areoid, have been related to the 6.1-mbar atmospheric pressure surface.

``To relate surface topography to atmospheric pressure, it is necessary to first compare planetary radii obtained from spacecraft occultations to those derived from MOLA. The occultations yield a measure of both planetary radius and atmospheric pressure and thereby a unique linkage between these two quantities... MOLA radii, which are considerably more accurate than radii obtained by occulatations, can then be related to occultation-derived surface pressures. By comparing MOLA radii to Viking and Mariner 9 occultations, Smith and Zuber (1988) showed that the zero point of MOLA topography corresponds to an atmospheric pressure of ~5.1 mbars at L_S=0^°. (L_S is the seasonal parameter on Mars and runs from 0^° to 360^° over the course of the Martian year; L_S=0 corresponds to the vernal equinox in the northern hemisphere.) The 6.1-mbar pressure level occurs at approximately ~1600 m relative to the zero reference of MOLA topography. However, the height of the 6.1-mbar surface needs to be adjusted, depending on the date. Seasonal variations in atmospheric pressure associated with the exchange of CO₂ between the atmosphere and polar caps is expected to produce vertical variations in the height of the 6.1-mbar surface of 1.5 to 2.5 km over the course of a Martian year.''

This paper also drives home the point that what one calls ``topography'' depends on the choice of spatial scale along the surface, something we have been conveniently sweeping under the carpet:

``... A power law of

V = 8.2815  l-1.9346 [km]2

provides a best fit to the [spatial] spectrum.... At 10-km spatial scale (degree 1000), a 10-m² RMS variance of topography is observed. The power law simply cannot be extroplated to estimate the RMS variance of surface roughness at landing-site scales. On a 100-m spatial scale only centimeters of roughness would be expected from extrapolation of the power law.''

How does this paragraph relate to the current problem? Simple. At the 10-km scale, the standard deviation of topography is about 3 m, or about 3 parts in 10^-4. Compare this to the value we calculated for m, namely ~ 0.0065.The root-mean slope at a scale » \leavevmode0 1/0 3 of our equatorial cell size is only » 0 1 /0 20 the difference between areocentric and areodetic latitude. Thus I have disproved Irina's claim that the difference between verticals on a spheroid and a sphere are negligible. Since both flavors of latitude are calculable in advance, and the difference does indeed matter, this effect should indeed be included in the cell_tab.

Angle from a cell normal to the spacecraft (the emission angle):

(Refer to Figure 2.) If the perpendicular to the topography in a cell is given by [^s]_C, the emission angle e is given by

cos e = ^ s   C  · ^ d   S/C  = ^ s   C  · é ë ® R   S/C  - ® r   C  || ® R   S/C  - ® r   C  || ù û

Irina has made a first attempt at calculating [^s]_C, although her results still are tainted (dominated) by cross-track errors. Presumably this problem will be worked out soon with the help of the María Zuber and her disciples.

Note: Although one can decompose [^s]_C into a 2-parameter (gradient + position) formalism, it is computationally more expedient to retain it as an (x, y, z) unit vector. The rationale for this is you eliminate the need for trig operations to decompose [^s]_C into its 3 components. Additionally, if later we decide to reference a cell to a spheroidal planet or to some other areoid, only recalculation of the normal to a cell's topography will be required.

We have already derived the distance from the cell to the spacecraft as ||d||, the magnitude of the vector [d\vec]. Now we proceed to calculate the nadir angle, y, and the position angle, ÐP.A., (defined in the usual sense as clockwise from North through East) for a given cell as seen from the spacecraft.

From Figure 3, it immediately follows that

cos y = ^ d   · ^ R   S/C   .

For the position angle, we employ a less obvious, but favorite trick of mine (Figure 4). Define [P\vec] to be the vector through the planet's North pole. E.g., in Cartesian coordinates it would be (0, 0, 1). The vector product

^ Q   1  = ^ P   × ^ R   S/C

is a vector 90^° around the sphere from both [P\vec] and [R\vec]. Now the vector product [^R]_S/C ×[^Q]₁ is a unit vector tangent to the surface at the nadir point, and pointing in the direction parallel to [P\vec].

Likewise, we form

^ Q   2  = ^ R   S/C  × ^ r   C

and cross it with [^R]_S/C to get a unit vector tangent to the surface at the nadir point and pointing toward [^r]_C.

Finally, we now take the dot product of the two vectors we have generated and get the result:

cos ÐP.A. = ( ^ R   S/C  × ^ Q   1  ) ·( ^ Q   2  × ^ R   S/C  )   .

This looks complicated, with four cross products and one dot product. However, recall that two of the three elements of [P\vec] are zero.

Þ We should be able to query spice directly to get the position angle, obviating the need for us to calculate it implicitly.

We now return now to the equation for nadir angle, and define the parameter z.

cos y = ^ d   · ^ R   S/C  º z .

When z < 0, an element is over the limb and can be ignored. The geometrical footprint is defined as those elements for which z = 0, and we can pick values for z such that a given percentage of the flux is contributed from within the cone specified. All we need do is to test if cos y is less than z, and we have determined which elements are in the footprint.

This leaves derivations d). and e). yet to be addressed. Derivation e). is a simple application of the barometric law, a.k.a. the hydrostatic equation, where the scale height H is defined as

H =  kT mg   .

Here, k is Boltzmann's constant, T is the temperature in degrees Kelvin, m is the mean molecular mass of the atmosphere, and g = GM/r² is the local acceleration due to gravity.

In general, T will be a function of altitude, season, position, and time of day on the planet. The biggest variation in T will be in the troposphere, which lies at the bottom of the first scale height from the ground. We should check, but my impression is that the stratospheric temperature is quite stable both diurnally and seasonally (barring a dust storm).

See the quote from Smith et al. on page 4.

d). Calculate the coordinates of the center of a cell from the spacecraft, i.e, a cell's GSH coordinates:

Matrix Approach (Not recommended)

This is rather an easy calculation, given that we know the vector [R\vec] to the spacecraft at a given time, and the vector [r\vec]_C to the center of the cell of interest. We recall our definition of [d\vec],

® d   = ® R   S/C  - ® r   C    .

All we have to do is calculate -[d\vec] in GSH coordinates. We achieve this via the following recipe.

·

First, a coordinate translation, i.e., bring the planetocentric frame origin coincident with the spacecraft frame origin.

·

Second, rotate about the +z-axis by the longitude of the ascending node of the spacecraft, W.

·

Third, rotate about the new +x-axis by an angle of +i, the inclination of the spacecraft orbit.

·

Fourth, rotate about the new +z-axis by an angle +n, the angle in the orbit plane since equator-crossing. inclination of the spacecraft orbit.

The same results can be achieved via use of quaternions. Quaternions have the added benefit of no singularities.

Review of Quaterion Algebra

Quaterions were first conceived by William Rowan Hamilton in 1843. They are hyper-complex numbers, or more precisely, complex numbers of rank 4. The set of quaternions, along with the two operations of addition and multiplication form a non-commutative division ring. This means that the quaternion product in general is not commutative, and also that the multiplicative inverse exists for every non-zero element in the set. Under the operations of multiplication and division, quaternions satisfy all the axioms of a field except for the commutative law of multiplication.

As the name suggests, a quaternion can be represented as a 4-tuple of real numbers, that is, for an element of R⁴:

q=(q0,q1,q2,q3)   ,

where q₀, q₁, q₂, q₃ are simply real numbers or scalars. The real part of a quaternion is defined as q₀. We then define a vector part, q, which is an ordinary vector in R³

q = iq1 + jq2 + kq3   .

Here i, j, and k are the standard orthonormal basis in R³. We can define a quaternion as the sum

q = q0 + q = iq1 + jq2 + kq3   .

A quaternion is basically the sum of a real number and a vector of three orthogonal imaginary parts. A strange beast-the sum of a scalar and a vector.

To multiply a (position) vector in R³ by a quaternion, we make a quaternion of the vector by entering its components as the hyperimaginary part, and setting the scalar part to 0. (A quaternion with zero scalar is entirely hyperimaginary and called a pure quaternion.) That is,

v Î R3 « v = 0 + v Î Q0 Ì Q

The conjugate of a quaternion q is q^*:

If    q = q0 + q    then    q* = q0 - q   ,

and the definition of quaternion multiplication is:

pq = p0q0 - p ·q + p0q +q0p + p×q

where the result is also a quaternion. ``·'' and ``×'' indicate vector dot and cross products.

The overall rotation of any vector v Î R³ is achieved by:

w = qvq*   ,

and multiple rotations by:

w = pqrvr*q*p* = (pqr)v(pqr)*   .

Quaternion multiplication goes from left to right, unlike matrix multiplication. We can write for any unit quaternion q,

q = q0 + q = cosq+ usinq

where q is the angle of rotation and u is the axis of rotation.

Finally, we define quaternion rotation operators as

Lq(v) = qvq*

which may be interpreted as a point or vector rotation, with respect to a fixed coordinate frame, and

Lq*(v) = q*vq

which is a coordinate frame rotation with respect to a fixed space of points or vectors.

Quaternion Approach

For the translation, we use the translation matrix T:

T= é ê ê ê ê ê ê ê ë 1 0 0 -x0 0 1 0 -y0 0 0 1 -z0 0 0 0 1 ù ú ú ú ú ú ú ú û   .

where x₀, y₀, and z₀ are the components of [R\vec]_S/C. Translation brings the origins of the two frames coincident.

Quaternion for rotation about +z-axis:        qz,W = cos[(W)/2] + ksin[(W)/2]        Quaternion for rotation about +x-axis:        qx,i = cosi/2 + isini/2        Quaternion for rotation about +z-axis:        qz,n = cos[(n)/2] + ksin[(n)/2]        The overall quaternion is then:        q º qz,W qx,i qz,n

While this sequence of rotations does solve the problem of transforming from Mars-fixed, Mars-centered coordinates to GSH-centered, GSH-fixed coordinates, there are two problems with it.

· First, the required angles, W and n are not really ``standard.'' That is, we have pretty much defined what they are, and it would take extra queries to the spice kernel to get them.

· Second, in the time elapsed between last south-to-north equator crossing and the time at which the last spectrum was acquired, several things transpire that have an effect on the relative position of the spacecraft and a cell. These must be taken into account, since neither the Mars-centered, Mars-fixed coordinate frame, nor the GSH-centered, GSH-fixed coordinate frame is an inertial frames. To list these effects: Mars has rotated beneath the spacecraft orbit due to its diurnal spin; Mars has rotated beneath the spacecraft by a smaller amount due to progress in its orbit around the Sun; the spacecraft orbit has regressed (by an even smaller amount) about the line of its nodes; if the orbit has any eccentricity whatsoever (it turns out it will) then the orbit will have precessed along the line of apsides. These last two effects are predominantly due to the oblateness (equatorial bulge) of Mars. However, there is a second set of regressions and precessions (even smaller in magnitude, but present nonetheless) due to solar torque on the spacecraft orbit. Other perturbations also exist (primarily due to Jupiter).

After a bit of head scratching, I have come up with an alternate scheme to treat the problem more accurately. It uses as inputs more ``classical'' (and accurately determined) parameters, and will account for all these additional rotations and perturbations at the same time.

With this alternative rotation scheme, we require as inputs only the sub-spacecraft (nadir) longitude and latitude at a given time, and the instantaneous spacecraft orbital inclination. This inclination is calculable from the Mars-centered, Mars-fixed spacecraft position and velocity vectors (i.e., the S/C state vector).

We proceed as follows.

·

First (as before), a coordinate translation, i.e., bring the planetocentric frame origin coincident with the spacecraft frame origin.

·

Rotate about the +z-axis by an angle l, the longitude of the subspacecraft point.

·

Rotate about the -y-axis by an angle b, the latitude of the S/C, as before.

·

Finally, rotate about the zenith-pointing +x-axis by an angle a, the instantaneous angle of the velocity vector to a parallel of latitude. Note: this angle is not the inclination of the orbit, except for the special case of at the ascending node.

Recall [P\vec], which we earlier defined to be the vector through the planet's North pole. In Mars-centered, Mars-fixed coordinates it is always (0, 0, 1). Now, the vector [P\vec] ×[R\vec]_S/C will point due eastward at any given time.

If the orbit were perfectly circular, then the angle a would directly follow from the dot product of this east-pointing vector and the spacecraft velocity vector, having normalized both vectors first. However, should there be a radial component to the spacecraft velocity, then the angle measured will be out of the plane tangent to the surface at the S/C's nadir. We must first determine the radial component of the spacecraft velocity and subtract it from the overall 3-dimensional velocity vector to get just the components of velocity perpendicular to [R\vec]_S/C.

It turns out this is rather easy. The magnitude of the radial component of velocity is given simply by

||vr|| = ^ R   S/C  · ® v

and its direction is along [R\vec]_S/C. Thus the nonradial S/C velocity is

® v   ^  = ® v   - ® v   r  = ||vr|| ^ R   S/C      ,

and the angle a follows as

a = cos-1 é ë ^ v   ^  · ® P   × ® R   S/C  || ® P   × ® R   S/C  || ù û     .

Thus we have the angle a for all time, regardless of orbital eccentricity or instantaneous perturbations due to asphericity and/or mascons in the solid body of Mars.

Summarizing the quaternions for this rotation sequence:

Quaternion for rotation about +z-axis:        q+z,l = cos[(l )/2] + ksin[(l )/2]        Quaternion for rotation about - y-axis:        q*- y,b = cos[(b)/2] -jsin[(b )/2]        Quaternion for rotation about +x-axis:        q*+x,a = cos[(a )/2] + isin[(a )/2]        The overall quaternion is then:        q º qz,l q*- y,(90° +b ) q*+x,a

Computer Implementation

Recall the definition of quaternion multiplication,