A&A 387, 1114-1122 (2002)
DOI: 10.1051/0004-6361:20020466

Transformation of spherical harmonic coefficients to ellipsoidal harmonic coefficients

D. Dechambre - D. J. Scheeres

Department of Aerospace Engineering, The University of Michigan, Ann Arbor, MI 48109-2140

Received 13 February 2002 / Accepted 8 March 2002

Abstract
Analytical expressions linking spherical harmonics gravity field expansions with ellipsoidal harmonics gravity field expansions are developed. Certain symmetries and simplifications for the transformation between the two are noted. Using the expressions, a numerical approach is developed and applied for the computation of ellipsoidal harmonic gravity coefficients using spherical harmonics coefficients as inputs. This method can be used to transform a measured spherical harmonic gravity field into an ellipsoidal harmonic gravity field for use within the Brillouin sphere of the original body. This would allow for the computation of spacecraft and natural trajectories close to the surface of an asteroid or comet using a measured gravity field.

Key words: comets: general - gravitation - minor planets, asteroids

1 Introduction

Small solar system bodies like asteroids and comets have become the target of current and forthcoming space missions. These bodies are known to have irregular shapes and complex gravity fields which can cause large perturbations to the trajectories of natural and artificial satellites in their vicinity. Thus, the measurement and use of a body's gravity field is of fundamental importance to understanding the environment close to its surface.

The classical approach for representing the gravitational field of an arbitrary body consists of expanding its gravitational potential into spherical harmonics (Kaula 1966). The distributed mass of the body is then modeled by the spherical harmonic coefficients C_lm and S_lm. These coefficients can be determined to high accuracy by evaluating the perturbations induced by the body on a spacecraft trajectory over the orbital phase of its mission. The advantage of this method is that it involves simple mathematics and converges to the correct gravity field outside the circumscribing (Brillouin) sphere about the body. In addition, finite truncations are often sufficient to match the "true potential'' of the body with good accuracy. Their drawback is that spherical harmonics expansions can exhibit severe divergence inside the circumscribing sphere. Thus they are not adequate to model the asteroid's gravity at close range.

Another approach consists of modeling the asteroid gravity field as a constant density polyhedron (Werner & Scheeres 1997). Small details of the asteroid surface can be included by modeling these regions with high resolution. The main advantage of this method is that the polyhedral potential is valid and exact for any given shape and density up to the surface of the body. Errors are thus reduced entirely to the errors in the asteroid shape determination and discretization. For Asteroid 433 Eros, a polyhedron model with 8200 plates was used to analyze close flybys of the surface (Antreasian et al. 1999). Despite the exactitude of the polyhedral potential, the constant density assumption can lead to erroneous gravity computations at close range. Improvements to match the "true gravity field'' more closely can be made by simulating density variations within the asteroid (Scheeres et al. 2000) but do not give a unique density distribution assignment. Thus, this approach cannot be used to develop a precision model of a body's gravitational field.

A third option for modeling a gravity field is the use of ellipsoidal harmonics (Garmier & Barriot 2001). While these harmonics still diverge when close to the body, their Brillouin ellipsoid lies much closer to the surface of the body. A disadvantage of the ellipsoidal harmonics is their complicated computation and the attendent difficulty to directly estimate these coefficients. Additionally, the estimation of spherical harmonic coefficients has been perfected to a high degree of accuracy, and thus one does not like to abandon their use altogether. To combine these approaches we propose a method for the transformation of measured spherical harmonic coefficients to an equivalent set of ellipsoidal harmonic coefficients. These gravity field expansions will then agree outside the Brillouin sphere of the body, and the ellipsoidal coefficients will be useable for computation of the gravity field within this sphere (but still outside of the Brillouin ellipsoid).

We first present the theory of ellipsoidal harmonics, introducing the ellipsoidal coordinates system and the Lamé functions. This theory was first developed by Byerly (1893), MacMillan (1930) and Hobson (1955), and was more recently applied by Garmier & Barriot (2001). The body of the paper describes a method for computing ellipsoidal harmonics directly from spherical harmonics. An analytical expression relating the ellipsoidal coefficients to the spherical coefficients is established. This result leads to a number of symmetry relations between ellipsoidal and spherical harmonic coefficients, and motivates the development of a numerical technique for transformation of spherical harmonic coefficients to ellipsoidal harmonic coefficients. Finally, numerical results are used to verify the cogency of our theory.

2 Harmonic expansions of the potential

2.1 Spherical harmonics expansion of the potential

The gravitational potential V of an arbitrary body can be expanded in a spherical harmonics expansion as:

 

$$V(r,\delta,\lambda)=\mu~\sum_{l=0}^{\infty}\sum_{m=0}^{l}{1\over r^{l+1}}~ P_{lm}(\sin\delta) [C_{lm}\cos(m\lambda)+S_{lm}\sin(m\lambda)]$$ (1)

where

-

r>0, $\delta\in[-{\pi\over2},{\pi\over2}]$, $\lambda\in[0,2\pi]$ are the usual spherical coordinates,

-

P_lm are the associated Legendre's polynomials and have the general expression

$$P_{lm}(\sin\delta)=(\cos\delta)^m \sum_{i=0}^{\left\lfloor{l-m\ov... ...\right\rfloor} {(-1)^i~(2l-2i)\over 2^li!(l-i)!(l-m-2i)!}(\sin\delta)^{l-m-2i},$$ (2)

-

the C_lm and S_lm have the dimension of a distance to the power l.

2.2 Ellipsoidal harmonics expansion of the potential

Similarly, the potential V can be expressed as a series of ellipsoidal harmonics involving the set of ellipsoidal coordinates $(\lambda_1,\lambda_2,\lambda_3)$.

   
2.2.1 Definition of the ellipsoidal coordinates

Ellipsoidal coordinates can be related in various ways to the Cartesian coordinates. They are defined with respect to a fundamental ellipsoid of semi-axes (a,b,c):

$${x^2\over a^2}+{y^2\over b^2}+{z^2\over c^2}=1 \quad\hbox{with}\quad a>b>c>0.$$ (3)

Based on this reference ellipsoid, we define the set of confocal quadrics by:

$${x^2\over {a^2+\theta^2}}+{y^2\over {b^2+\theta^2}}+{z^2\over {c^2+\theta^2}}=1,$$ (4)

which can be rewritten as

$${x^2\over {\lambda^2}}+{y^2\over {\lambda^2-h^2}}+{z^2\over {\lambda^2-k^2}}=1,$$ (5)

with

$$\left\{ \begin{array}{lll} \lambda^2&=a^2+\theta^2,\\ h^2&=a^2-b^2,\\ k^2&=a^2-c^2. \end{array} \right.$$ (6)

The advantage of rewriting Eq. (4) in the form (5) is that the family of confocal quadrics to the fundamental ellipsoid is now entirely characterized by the two parameters hand k. If we now fix x, y and z, Eq. (5) has three roots $\lambda_1$, $\lambda_2$ and $\lambda_3$ such that:

$${\lambda_1}^2>k^2>{\lambda_2}^2>h^2>{\lambda_3}^2>0.$$

The set $(\lambda_1,\lambda_2,\lambda_3)$ is orthogonal and is defined as the set of ellipsoidal coordinates. The coordinate $\lambda_1$ is often called the "elliptic radius'' by analogy with the spherical radius r since $\lambda_1=\hbox{Constant}$ is the equation of the surface of an ellipsoid in ellipsoidal coordinates.

According to Eq. (5), there are eight points in space corresponding to the same $(\lambda_1^2,\lambda_2^2,\lambda_3^2)$. These points are defined by $(\pm x,\pm y,\pm z)$ in Cartesian coordinates. Cartesian coordinates can be related to ellipsoidal coordinates using the following transformation:

   

$${\lambda_1^2~\lambda_2^2~\lambda_3^2\over h^2k^2}$$ x² = (7)

$${(\lambda_1^2-h^2)(\lambda_2^2-h^2)(h^2-\lambda_3^2)\over h^2(k^2-h^2)}$$ y² = (8)

$${(\lambda_1^2-k^2)(k^2-\lambda_2^2)(k^2-\lambda_3^2)\over k^2(k^2-h^2)}\cdot$$ z² = (9)

2.2.2 Introduction to the Lamé functions

In the following, we will deal with spaces whose boundaries are ellipsoids and solutions to Laplace's equation within these spaces. It is then convenient to express Laplace's equation $\nabla^2V=0$ in terms of the ellipsoidal coordinates  $(\lambda_1,\lambda_2,\lambda_3)$.

In the othogonal set of coordinates $(\lambda_1,\lambda_2,\lambda_3)$, Laplace's equation takes the form:

 

$$\left.\begin{array}{l} (\lambda_2^2-\lambda_3^2)\sqrt{(\lambda_1^... ...isplaystyle{\partial{V}\over \partial\lambda_3} \right) \end{array}\right\} =0.$$ (10)

If we look for normal solutions of the form $V(\lambda_1,\lambda_2,\lambda_3)=R(\lambda_1)M(\lambda_2)N(\lambda_3)$, then the functions R, M and N must satisfy Lamé's equation

 

$$(\lambda_i^2-h^2)(\lambda_i^2-k^2){{\rm d}^2E\over {\rm d}\lambda... ...a_i(2\lambda_i^2-h^2-k^2){{\rm d}E\over {\rm d}\lambda_i} +(K-H\lambda_i^2)E=0,$$ (11)

where H and K are constants and are called the separation parameters.

The parameters H and K can be chosen such that solutions to (11) take the general form:

$$E_n^p(\lambda)=\psi(\lambda)~P(\lambda),$$ (12)

where $\psi$ is a leading product and P is a polynomial in $\lambda^2$. The leading product $\psi$ is of one of the following forms:

$$\lambda^u~{\sqrt{\vert\lambda^2-h^2\vert}}^v~{\sqrt{\vert\lambda^2-k^2\vert}}^w \quad\hbox{with}\quad u, v, w \in \{0,1\}.$$ (13)

E_n^p is called a Lamé function of the first kind of degree n and order p. The index n corresponds to the degree of $\psi(\lambda)~P(\lambda)$ in $\lambda$. For a given degree n, there are (2n+1) Lamé functions of the first kind where $p=1,\dots,(2n+1)$.

There are four types of Lamé functions of the first kind K, L, M and N defined according to the form of the leading product $\psi$.
Let us define $\sigma$ as:

$$\sigma=\left\{ \begin{array}{ll} {1\over 2}~n & \quad\hbox{fo... ... 2}~(n-1) & \quad\hbox{for $n$\space odd.} \end{array} \right.$$ (14)

Then, for a given n, there are:

- $(\sigma+1)$ functions E_n^p of type K, $p=1,\dots,(\sigma+1)$, defined by:

 

$$K_n^p(\lambda)=a_{0p}\lambda^n+a_{1p}\lambda^{n-2}+\dots+\left\{ ... ...pace even}\\ a_{{\sigma}p}\lambda &\hbox{for $n$\space odd}\end{array} \right.$$ (15)

- $(n-\sigma)$ functions E_n^p of type L, $p=(\sigma+2),\dots,(n+1)$, defined by:

 

$$L_n^p(\lambda)=\sqrt{\vert\lambda^2-h^2\vert} \bigg[b_{0p}\lambda... ... even}\\ b_{{n-\sigma-1},p} &\hbox{for $n$\space odd}\end{array} \right.\bigg]$$ (16)

- $(n-\sigma)$ functions E_n^p of type M, $p=(n+2),\dots,(2n-\sigma+1)$, defined by:

 

$$M_n^p(\lambda)=\sqrt{\vert\lambda^2-k^2\vert} \bigg[c_{0p}\lambda... ... even}\\ c_{{n-\sigma-1},p} &\hbox{for $n$\space odd}\end{array} \right.\bigg]$$ (17)

- $\sigma$ functions E_n^p of type N, $p=(2n-\sigma+2),\dots,(2n+1)$, defined by:

 

$$N_n^p(\lambda)=\sqrt{\vert(\lambda^2-h^2)(\lambda^2-k^2)\vert}\bi... ...}\\ d_{{\sigma-1},p}\lambda &\hbox{for $n$\space odd}\end{array}\right.\bigg].$$ (18)

As suggested by the name Lamé function of the first kind, there exists a second kind of Lamé functions. The normal solution $V(\lambda_1,\lambda_2,\lambda_3)=E_n^p(\lambda_1)E_n^p(\lambda_2)E_n^p(\lambda_3)$is indeed inappropriate to describe an external potential since $E_n^p(\lambda_1)\to\infty$ when $\lambda_1\to\infty$. The Lamé functions of the second kind are then defined from the Lamé functions of the first kind by:

 

$$F_n^p(\lambda_1)=(2n+1)~E_n^p(\lambda_1) \int_{\lambda_1}^\infty{{\rm d}u\over(E_n^p(u))^2~\sqrt{(u^2-k^2)(u^2-h^2)}},$$ (19)

such that $F_n^p(\lambda_1)\to 0$ when $\lambda_1\to\infty$.

2.2.3 Ellipsoidal harmonics expansion of the potential

The Lamé products $E_n^p(\lambda_1)E_n^p(\lambda_2)E_n^p(\lambda_3)$ and $F_n^p(\lambda_1)E_n^p(\lambda_2)E_n^p(\lambda_3)$ are both normal solutions to Lamé's Eq. (11). The former is continuous within any ellipsoid $\lambda_1=\lambda_1^{{\rm ref}}$. The latter is continuous for $\lambda_1\geq\lambda_1^{{\rm ref}}$ and it vanishes when $\lambda_1\to\infty$.

We can then define the ellipsoidal harmonics expansion of the potential anywhere in space as:

  

$$V(\lambda_1,\lambda_2,\lambda_3) \mu\sum_{n=0}^{\infty}\sum_{p=1}^{2n+1}\alpha_{np}{E_n^p(\lambda_... ...ef})} E_n^p(\lambda_2)E_n^p(\lambda_3),\quad\lambda_1\leq\lambda_1^{{\rm ref}},$$ = (20)

$$V(\lambda_1,\lambda_2,\lambda_3) \mu\sum_{n=0}^{\infty}\sum_{p=1}^{2n+1}\alpha_{np}{F_n^p(\lambda_... ...ef})} E_n^p(\lambda_2)E_n^p(\lambda_3),\quad\lambda_1\geq\lambda_1^{{\rm ref}},$$ = (21)

where the $\alpha_{np}$ are constants. On the ellipsoid $\lambda_1=\lambda_1^{{\rm ref}}$, (20) and (21) have identical expressions so that the overall potential ( $\lambda_1\leq\lambda_1^{{\rm ref}}$ and $\lambda_1\geq\lambda_1^{{\rm ref}}$) is continuous.

2.3 Analogies between the two types of expansion

The indices (l,m) and (n,p) are called degree and order of the spherical harmonics expansion, respectively ellipsoidal harmonics expansion. For a given degree l=n, there are as many spherical harmonics (2l+1) as ellipsoidal harmonics (2n+1). Furthermore, there is a strong analogy in the expressions of the harmonics themselves that we will discuss in the followings.

When $\lambda_1$ is very large,

$$F_n^p(\lambda_1)={E_n^p(\lambda_1)\over \lambda_1^{2n+1}}\cdot$$ (22)

In addition $E_n^p(\lambda_1)\sim c_0~\lambda_1^n,\;\lambda_1\to\infty$, so that

$$F_n^p(\lambda_1)={c_0\over\lambda_1^{n+1}},\quad\lambda_1\to\infty,$$ (23)

is analogous to the term $${1\over r^{l+1}}$$ in the spherical harmonic expansion.

Similarly, one can relate the terms

$$P_{lm}\sin(\delta)\left\{\begin{array}{ll}\cos(m\lambda)\\ \s... ...}\right.\quad\hbox{and}\quad E_n^p(\lambda_2)E_n^p(\lambda_3)$$ (24)

to one another since they are functions of the analogous sets of variables $(\delta,\lambda)$ and $(\lambda_2,\lambda_3)$ respectively.

The explicit parallelism between both expansions leads us to investigate any relation between the ellipsoidal harmonics coefficients, $\alpha_{np}$, and the spherical harmonics coefficients, C_lm and S_lm. This will be the object of the next section.

3 From spherical harmonics to ellipsoidal harmonics

In the following, we want to express the normalized ellipsoidal harmonics coefficients, $\alpha_{np}$, as a function of the spherical harmonics coefficients C_lm and S_lm.

3.1 Orthogonality

 

One can show, using Green's theorem, that the Lamé functions of the first kind satisfy the following relation:

 

$$\int_{E_{\lambda_1}}{E_n^p(\lambda_2)E_n^p(\lambda_3)E_{n'}^{p'}(... ...ad\hbox{if $n=n'$\space and $p=p'$ } \\ 0\quad\hbox{else}\end{array} \right. ,$$ (25)

where $E_{\lambda _1}$ is the surface of the ellipsoid defined by $\lambda_1$ and ${\rm d}S$ is a surface element on $E_{\lambda _1}$. The Lamé products $E_n^p(\lambda_2)E_n^p(\lambda_3)$ are called surface harmonics. The relation (25) then corresponds to an orthogonalization property for the surface harmonics.

Applying (25) to both sides of (21) yields:

$$\mu~\alpha_{np}{F_n^p(\lambda_1)\over F_n^p(\lambda_1^{{\rm ref}}... ...da_3)\over \sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}~ {\rm d}S$$ (26)

so that $\alpha_{np}$ can be isolated as

 

$$\alpha_{np}={1\over\gamma_n^p}~{F_n^p(\lambda_1^{{\rm ref}})\over... ...a_3)\over \sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}~ {\rm d}S.$$ (27)

Now, if we replace the potential V by its spherical harmonics expansion (1), one can express $\alpha_{np}$ as a linear transformation of the C_lm and S_lm:

$$\alpha_{np}=\sum_{l=0}^{\infty}\sum_{m=0}^l (A_{np}^{lm}\cdot C_{lm} + B_{np}^{lm}\cdot S_{lm})$$ (28)

where the coefficients A_np^lm and B_np^lm are defined by

  

$$A_{np}^{lm}={1\over\gamma_n^p}{F_n^p(\lambda_1^{{\rm ref}})\over ... ...rm d}S\over {r^{l+1}\sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}}$$ (29)

$$B_{np}^{lm}={1\over\gamma_n^p}{F_n^p(\lambda_1^{{\rm ref}})\over ... ...S\over {r^{l+1}\sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}}\cdot$$ (30)

Some precautions have to be taken for the choice of the ellipsoid $E_{\lambda _1}$ in the above development. This will be developed in Sect. 3.2.3 Convergence considerations.

3.2 Symmetry

 

Consider the term

$${\rm d}S\over {r^{l+1}\sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}$$ (31)

in Eqs. (29) and (30). We have seen in Sect. 2.2.1 that there are eight points on the surface of the ellipsoid $E_{\lambda _1}$ coresponding to the same set of values $(\lambda_1^2,\lambda_2^2,\lambda_3^2)$. Let denote these points by $(P_j)_{j=1,\dots,8}$. They also have identical values for x², y², and z². Therefore:

• they have the same spherical radius r,

• their latitude is $\pm\delta$,

• their longitude is $\lambda$, $\pi-\lambda$, $\pi+\lambda$ or $2\pi-\lambda$.

Consequently, the term ${{\rm d}S\over {r^{l+1}\sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}}$ is identical for all the $(P_j)_{j=1,\cdots,8}$ as long as we choose identical elementary surface area ${\rm d}S$.

Consider now the other term,

$$P_{lm}(\sin\delta)~ \left\{ \begin{array}{l} \cos (m\lambda)\... ...ambda)\end{array} \right\} ~ E_n^p(\lambda_2)E_n^p(\lambda_3).$$ (32)

It has the same norm for each of the $(P_j)_{j=1,\cdots,8}$. By arguing on its sign, one can then deduce some interesting properties for the $\alpha_{np}$. Especially, we will see in the following that the coefficients A_np^lm and B_np^lm may vanish depending on the parity of n, p, l and m.

We will treat as an example the case where n is even, $1\leq p\leq (\sigma+1)$. For all the other cases, we will summarize the results in a table.

 

 

Table 1: Arguments of symmetry for the particular case n even and $1\leq p\leq ({n\over 2}+1)$.

			l and m even		l and m odd		l even and m odd		l odd and m even
			sign of Eq.		sign of Eq.		sign of Eq.		sign of Eq.
Pi	$\delta_{P_i}$	$\lambda_{P_i}$	$(\ref{eq:term1})$	$(\ref{eq:term2})$	$(\ref{eq:term1})$	$(\ref{eq:term2})$	$(\ref{eq:term1})$	$(\ref{eq:term2})$	$(\ref{eq:term1})$	$(\ref{eq:term2})$
P1	$\delta$	$\lambda$	+	+	+	+	+	+	+	+
P2	$\delta$	$\pi-\lambda$	+	-	-	+	-	+	+	-
P3	$\delta$	$\pi+\lambda$	+	+	-	-	-	-	+	+
P4	$\delta$	$2\pi-\lambda$	+	-	+	-	+	-	+	-
P5	$-\delta$	$\lambda$	+	+	+	+	-	-	-	-
P6	$-\delta$	$\pi-\lambda$	+	-	-	+	+	-	-	+
P7	$-\delta$	$\pi+\lambda$	+	+	-	-	+	+	-	-
P8	$-\delta$	$2\pi-\lambda$	+	-	+	-	-	+	-	+

3.2.1 Example

We are here going to discuss the sign of

  

$$P_{lm}(\sin\delta)\cos(m\lambda)E_n^p(\lambda_2)E_n^p(\lambda_3),$$ (33)

$$\hbox{and}\quad P_{lm}(\sin\delta)\sin(m\lambda)E_n^p(\lambda_2)E_n^p(\lambda_3).$$ (34)

For n even and $1\leq p\leq (\sigma+1)=({n\over 2}+1)$, E_n^p is of type K and thus it is an even polynomial in $\lambda$. Therefore, the sign of $E_n^p(\lambda_2)E_n^p(\lambda_3)$ is identical for all the $(P_j)_{j=1,\dots,8}$. We then only need to argue on the sign of $P_{lm}(\sin\delta)\cos(m\lambda)$ and $P_{lm}(\sin\delta)\sin(m\lambda)$.

 

 

Table 2: Table of the A_np^lm.

type	n even				n odd
of	l even		l odd		l even		l odd
Enp	m even	m odd	m even	m odd	m even	m odd	m even	modd
K	x	0	0	0	0	0	0	x
L	0	0	0	0	0	0	0	0
M	0	x	0	0	0	0	x	0
N	0	0	0	0	0	0	0	0

 

 

Table 3: Table of the B_np^lm.

type	n even				n odd
of	l even		l odd		l even		l odd
Enp	m even	m odd	m even	m odd	m even	m odd	m even	modd
K	0	0	0	0	0	0	0	0
L	x	0	0	0	0	0	0	x
M	0	0	0	0	0	0	0	0
N	0	x	0	0	0	0	x	0

Let us distinguish the four cases:

• For l and m even, P_lm is an even polynomial in $\sin\delta$. The terms (33) and (34) have constant norms over all the P_j's and their signs are listed in Table 1. Therefore, in the summation (30), the P_j's contributions cancel each other over the surface of the ellipsoid $E_{\lambda _1}$ and B_np^lm=0. On the other hand, the A_np^lmare a priori non zero. They can be rewritten as:
  

  $$A_{np}^{lm}={8\over\gamma_{np}}{F_n^p(\lambda_1^{{\rm ref}})\over... ...r {r^{l+1}\sqrt{(\lambda_1^2-\lambda_2^2)(\lambda_1^2-\lambda_3^2)}}}~{\rm d}S,$$ (35)

  
  
  where $(E_{\lambda_1})^+$ denotes the octant of the ellipsoid $E_{\lambda _1}$ corresponding to (X>0,Y>0,Z>0). This integral is a priori non zero, however it may vanish for some n, p, l and m.
• For l and m odd, P_lm is an even polynomial in $\sin\delta$. It follows from Table 1 that A_np^lm=B_np^lm=0.
• For l even and m odd, P_lm is an odd polynomial in $\sin\delta$. It follows from Table 1 that A_np^lm=B_np^lm=0.
• For l odd and m even, P_lm is an odd polynomial in $\sin\delta$. It follows from Table 1 that A_np^lm=B_np^lm=0.

3.2.2 Summary of results

 

By going through all the cases individually, we obtain the following tables for the A_np^lm and B_np^lm.

Using Tables 2 and 3, the relation between ellipsoidal harmonics coefficients and spherical harmonics coefficients can be simplified. Depending on the values of the degree n and the order p, Eq. (28) reduces to:

• For n even,

  
  

  $$\noalign{\hbox{-- for $1\leq p\leq {n+2\over 2}$ ,}} \alpha_{np} \sum_l \sum_m A_{np}^{2l,2m}\cdot C_{2l,2m},$$ = (36)

  $$\noalign{\hbox{-- for ${n+4\over 2}\leq p\leq (n+1)$ ,}} \alpha_{np} \sum_l \sum_m B_{np}^{2l,2m}\cdot S_{2l,2m},$$ = (37)

  $$\noalign{\hbox{-- for $(n+2)\leq p\leq {3n+2\over 2}$ ,}} \alpha_{np} \sum_l \sum_m A_{np}^{2l,2m+1}\cdot C_{2l,2m+1},$$ = (38)

  $$\noalign{\hbox{-- for ${3n+4\over 2}\leq p\leq (2n+1)$ ,}} \alpha_{np} \sum_l \sum_m B_{np}^{2l,2m+1}\cdot S_{2l,2m+1}.$$ = (39)

  
  

• For n odd,

  
  

  $$\noalign{\hbox{-- for $1\leq p\leq {n+1\over 2}$ ,}} \alpha_{np} \sum_l \sum_m A_{np}^{2l+1,2m+1}\cdot C_{2l+1,2m+1},$$ = (40)

  $$\noalign{\hbox{-- for ${n+3\over 2}\leq p\leq (n+1)$ ,}} \alpha_{np} \sum_l \sum_m B_{np}^{2l+1,2m+1}\cdot S_{2l+1,2m+1},$$ = (41)

  $$\noalign{\hbox{-- for $(n+2)\leq p\leq {3n+3\over 2}$ ,}} \alpha_{np} \sum_l \sum_m A_{np}^{2l+1,2m}\cdot C_{2l+1,2m},$$ = (42)

  $$\noalign{\hbox{-- for ${3n+5\over 2}\leq p\leq (2n+1)$ ,}} \alpha_{np} \sum_l \sum_m B_{np}^{2l+1,2m}\cdot S_{2l+1,2m}.$$ = (43)

  
  

3.2.3 Convergence considerations

 

When considering an attracting body, the spherical harmonics expansion of the potential is uniformly convergent outside the Brillouin sphere. The Brillouin sphere is defined as the smallest sphere that encloses the attracting body. Inside the Brillouin sphere, severe divergences of the spherical harmonics expansion may occur.
Similary, the ellipsoidal harmonics expansion of the potential is uniformly convergent outside the Brillouin ellipsoid. The Brillouin ellipsoid is defined as the smallest ellipsoid that encloses the attracting body. Inside the Brillouin ellipsoid, severe divergences of the ellipsoidal harmonics expansion may occur.

In Sect. 3.1, the linear transformation (28),

$$\alpha_{np}=f(C_{lm},S_{lm}),$$ (44)

has been obtained by replacing the potential V by its spherical harmonics expansion (1) into (27). Since convergence of the spherical harmonics expansion of the potential is certain only outside the Brillouin sphere, some precautions have to be taken for the choice of the ellipsoid $E_{\lambda _1}$over which the spherical harmonics expansion of the potential has to be evaluated in Eq. (27). The ellipsoid $E_{\lambda _1}$ has to be chosen such that it encloses the Brillouin sphere.

4 An application: Computation of ellipsoidal gravity field harmonics for small solar system bodies

Ellipsoidal harmonics are more appropriate than spherical harmonics to compute the gravity field of irregularly shaped attracting bodies. Since an arbitrary shape body is better fitted by an ellipsoid than by a sphere, the region of divergence of an ellipsoidal harmonics expansion is smaller than the one of a spherical harmonics expansion (see paragraph 3.2.3).

In the following, we are going to use the method described above for the evaluation of the ellipsoidal coefficients $\alpha_{np}$. Given spherical coefficients C_lm and S_lm, we are then going to "reconstruct'' the corresponding spherical harmonics expansion of the potential with an ellipsoidal harmonics expansion.

4.1 Numerical methods for ellipsoidal harmonics

Expression of the potential in Eq. (21) requires the knowledge of the ellipsoidal coefficients $\alpha_{np}$, the Lamé functions of the first and second kind and the ellipsoidal coordinates.

Lamé functions of the first kind E_n^p have explicit expressions up to n=3. For n>3, Ritter and Dobner have transformed the problem of computing the polynomial coefficients in Eqs. (15), (16), (17) and (18) into an eigenvalue problem (see Ritter 1998 and Dobner & Ritter 1998). Regarding Lamé functions of the second kind F_n^p defined by (19), they can be evaluated numerically using a Gauss-Legendre quadrature (see Garmier & Barriot 2001).
The ellipsoidal coordinates $(\lambda_1,\lambda_2,\lambda_3)$ have been defined as the roots of Eq. (5). They can be explictly expressed in terms of the Cartesian coordinates (x,y,z). However these solutions may be singular and we rather use a numerical solver such as a false position method (see Press et al. 1993) to determine the roots of (5).

Finally, the ellipsoidal coefficients $\alpha_{np}$ have been constructed from the spherical coefficients C_lmand S_lm by the linear transformation (28). We have seen in Sect. 3.2 that the coefficients A_np^lm and B_np^lm vanish in many cases. However, when a priori non zero, no explicit expression exists. Instead, we "tessellate'' the surface of the ellipsoid $E_{\lambda _1}$ into flat triangular faces (Fig. 1). The vertices of theses faces are placed on the surface of the ellipsoid with constant latitude and longitude increments. To first approximation, the integrands in expressions (29) and (30) are treated as constant over each plate so that the integral over the surface $E_{\lambda _1}$ is transformed into a finite summation over the number of plates:

 

$$A_{np}^{lm}={1\over\gamma_n^p}{F_n^p(\lambda_1^{{\rm ref}})\over ... ...t{({{\lambda_1}_i}^2-{{\lambda_2}_i}^2)({{\lambda_1}_i}^2-{{\lambda_3}_i}^2)}}}$$ (45)

 

$$B_{np}^{lm}={1\over\gamma_n^p}{F_n^p(\lambda_1^{{\rm ref}})\over ... ...\lambda_1}_i}^2-{{\lambda_2}_i}^2)({{\lambda_1}_i}^2-{{\lambda_3}_i}^2)}}}\cdot$$ (46)

Regarding the normalization constants $\gamma_n^p$ defined by Eq. (25), they can be rewritten as:

$$\gamma_n^p= 8\int_0^h \int_h^k {(\lambda_2^2-\lambda_3^2)(E_n^p(\... ..._2^2-h^2)(h^2-\lambda_3^2)(k^2-\lambda_3^2)}}{\rm d}\lambda_2 {\rm d}\lambda_3.$$ (47)

Garmier & Barriot (2001) have developed a method to numerically compute the above integral.

	Figure 1: "Tessellated surface'' of the ellipsoid $E_{\lambda _1}$.
Open with DEXTER

4.2 Numerical computation of the external potential for asteroid 433 Eros

 

In the following, we focus on the computation of the external gravity field for asteroid 433 Eros (Yeomans et al. 2000). We assume the knowledge of the spherical coefficients C_lm and S_lm for Eros. They completely describe the spherical gravity field of the asteroid. Using Eqs. (21) and (28), we will compute the ellipsoidal harmonics expansion for Eros. We will then compare the resulting ellipsoidal gravity field with the given spherical gravity field. We expect them to match outside the Brillouin sphere where spherical and ellipsoidal harmonics expansions are both convergent.

The Brillouin ellipsoid is defined as the ellipsoid that best fits the asteroid. In Eqs. (20) and (21), we choose $\lambda_1^{{\rm ref}}$ equal to the semi-major axis of the Brillouin ellipsoid. With this choice of $\lambda_1^{{\rm ref}}$, (20) and (21) respectively define the interior and the exterior potential of the Brillouin ellipsoid. For our simulations, we choose a Brillouin ellipsoid for asteroid Eros of semi-axes $17.6 \hbox{ km} \times 8.6 \hbox{ km} \times 6.1 \hbox{ km}$ and a Brillouin sphere of radius $17.6 \hbox{ km}$. These data are computed from a discretized shape model of the surface of the asteroid.

Per Sect. 3.2.3, the "tesselated'' ellipsoid $E_{\lambda _1}$ used for the numerical evaluation of the ellipsoidal coefficients $\alpha_{np}$ must be chosen such that it encloses the Brillouin sphere. In order to minimize the error between its "tesselated'' surface and its "true'' surface, we choose for $E_{\lambda _1}$ the smallest ellipsoid confocal to the Brillouin ellipsoid that encloses the Brillouin sphere: its semi-axes are $24.1 \hbox{ km} \times 18.6 \hbox{ km} \times 17.6 \hbox{ km}$.

	Figure 2: Maximum error on the computation of the potential for Eros.
Open with DEXTER

	Figure 3: Weighted error on the computation of the potential for Eros.
Open with DEXTER

	Figure 4: Errors $\epsilon _x$, $\epsilon _y$ and $\epsilon _z$ on the computation of the acceleration for Eros.
Open with DEXTER

Given a spherical harmonics expansion of degree 6 (C_lm and S_lm up to l=6), we compute the ellipsoidal potential $V_{\rm e}$ for various degrees n of the ellipsoidal harmonics expansion up to n=12. We then evaluate the error $\epsilon$ between $V_{\rm e}$ and the spherical potential of degree 6 taken as the reference, $V_{{\rm ref}}$:

$$\epsilon={\vert V_{\rm e}-V_{{\rm ref}}\vert\over V_{{\rm ref}}}\cdot$$ (48)

The evaluation of the potentials $V_{\rm e}$ and $V_{{\rm ref}}$ is performed on the faces of the "tesselated'' ellipsoid $E_{\lambda _1}$ and the resulting error $\epsilon$ is represented as a function of the degree n on Figs. 2 and 3 by means of:

-

a maximum error defined as $\max(\epsilon)$ over all the faces,

-

a weighted error defined as

$${\sum_i\epsilon_i\;{\rm d}S_i\over\sum_i{\rm d}S_i},$$

where $\epsilon_i$ denotes the error evaluated on the face i of area ${\rm d}S_i$.

The error $\epsilon$ is plotted on Figs. 2 and 3 for three "tesselated'' surface models of $E_{\lambda _1}$:

-

1520 faces, corresponding to increments in latitude $\Delta\delta$ and longitude $\Delta\lambda$ of $9\deg$,

-

6240 faces, corresponding $\Delta\delta=\Delta\lambda=4.5\deg$,

-

14 160 faces, corresponding $\Delta\delta=\Delta\lambda=3\deg$.

As the degree of the ellipsoidal harmonics expansion increases, we expect the error $\epsilon$ to decrease to zero. However, according to the plots, the reconstruction of the spherical gravity field by ellipsoidal harmonics depends on the level of discretization of the ellipsoid $E_{\lambda _1}$ for the computation of the ellipsoidal coefficients $\alpha_{np}$. Indeed evaluating the coefficients A_np^lm and B_np^lm themselves, we observe that the theoretical results summarized in Table 1 are violated when the "tesselated'' surface does not sufficiently match the true surface of the ellipsoid $E_{\lambda _1}$.

4.3 Numerical computation of the external acceleration for asteroid 433 Eros

 

Differentiating Eqs. (1) and (21) with respect to x, y and z, one can obtain the expressions of the spherical and ellipsoidal accelerations vectors. The ellipsoidal acceleration components can be analytically written as:

$$\partial V_{\rm e}\over \partial x_i \mu\sum_{n=0}^{\infty}\sum_{p=1}^{2n+1}{\alpha_{np}\over F_n^p\le... ...\partial\lambda_1}-{2n+1\over \sqrt{(\lambda_1^2-h^2)(\lambda_1^2-k^2)}}\right]$$ =

$$\times \left({\partial\lambda_1\over\partial x_i}\right)+ F_n^p(\... ...artial\lambda_3}\right)\left({\partial\lambda_3\over\partial x_i}\right)\Bigg\}$$ (49)

The $\partial E_n^p(\lambda_j)\over\partial\lambda_j$ are easily determined by rewriting the $E_n^p(\lambda_j)$ in the form $\psi(\lambda_j)P(\lambda_j)$ where $\psi$ is a leading product and P a polynomial in $\lambda_j^2$. The $\partial\lambda_j\over \partial x_i$ are obtained by differentiating the set of Eqs. (7), (8), (9) with respect to the x_i's.

With the model described in Sect. 4.2, we evaluated the errors $\epsilon _x$, $\epsilon _y$ and $\epsilon _z$ between the components of the ellipsoidal acceleration of degree n and the components of the spherical acceleration of degree 6 taken as a reference defined as:

$$\epsilon_x={\left\vert{\partial V_{\rm e}\over\partial x}-{\p... ...\right\vert\over {\partial V_{{\rm ref}}\over\partial x}}\cdot$$ (50)

The errors $\epsilon _x$, $\epsilon _y$ and $\epsilon _z$ are computed on a "tesselated model'' of ellipsoid $E_{\lambda _1}$ (of semi-axes $24.1 \hbox{ km} \times 18.6 \hbox{ km} \times 17.6 \hbox{ km}$) with 6240 faces and are plotted as a function of the degree n in Fig. 4.

5 Conclusion

Ellipsoidal harmonics are useful for the representation of asteroid and comet gravity fields. Our approach is to use our knowledge of spherical harmonics to provide a better understanding of the ellipsoidal harmonics theory. Because of the obvious analogies between spherical and ellipsoidal harmonics, we have been interested in establishing an analytical expression that relates the ellipsoidal harmonic coefficients $\alpha_{np}$ to the spherical harmonic coefficients C_lm and S_lm. Though no explicit formulation exists for Eq. (28), we have shown that it can be further simplified based on symmetry arguments and can be used to develop a numerical transformation procedure. Plots have been provided to validate the computation of ellipsoidal harmonics using spherical harmonics.

References

 
• Antreasian, P., Helfrich, C., Miller, J., et al. 1999, Preliminary Planning for NEAR's Low-Altitude Operations at 433 Eros, AAS, 99-465 In the text  
• Bierly, W. 1893, An Elementary Treatise on Fourier's Series and Spherical, Cylindrical, and Ellipsoidal Harmonics (Ginn and Company) In the text  
• Dobner, H.-J., & Ritter, S. 1998, Verified Computation of Lamé Functions with High Accuracy, Computing 60, 81 In the text  
• Garmier, R., & Barriot, J.-P. 2001, Celestial Mechanics and Dynamical Astronomy, 79, 235 In the text  
• Hobson, E. 1955, The Therory of Spherical and Ellipsoidal Harmonics (Chelsea Publishing Company) In the text  
• Kaula, W. 2000, Theory of Satellite Geodesy, Applications of Satellites to Geodesy (Dover)  
• MacMillan, W. D. 1930, The Theory of the Potential (McGraw-Hill) In the text  
• Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 1993, Numerical Recipes in C: The Art of Scientific Computing (Cambridge) In the text  
• Ritter, S. 1998, Zamm$\cdot$ Z. Angew. Math. Mech. 78, 1, 66 In the text  
• Scheeres, D. J., Khushalani, B., & Werner, R. A. 2000, Planet. Space Sci., 48, 965 In the text  
• Werner, R., & Scheeres, D. 1997, Celestial Mechanics and Dynamical Astronomy, 65, 313 In the text NASA ADS  
• Yeomans, D., Antreasian, P., Barriot, J.-P., et al. 2000, Science, 289, 2085 In the text NASA ADS