% IMPORTANT NOTES:
% - figures available on request from the authors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle[11pt]{spacekap}
\hyphenation{sur-fa-ce}
\hyphenation{atmo-sphe-re}
\begin{opening}
\title{ALBEDO PERTURBATION MODELS~: \\ GENERAL FORMALISM AND APPLICATIONS
 TO LAGEOS}
\runningtitle{ALBEDO PERTURBATION MODELS}
\author{DAVID VOKROUHLICK\'Y\thanks{On leave from
Astronomical Institute, 
Charles University, \v Sv\'edsk\'a 8, 15000 Prague 5, Czechoslovakia}}
\institute{O.C.A., dept. C.E.R.G.A., Av N. Copernic, 06130 Grasse, France}
\author{PAOLO FARINELLA}
\author{DAVID LUCCHESI}
\institute{Dip. di Matematica, Univ. di Pisa, Via Buonarroti 2, 56127 Pisa,
 Italy}
\end{opening}
\begin{document}
\begin{abstract}
The force due to radiation pressure on a satellite of
arbitrary shape is written in a general form within a formalism
similar to that used in the theory of radiative transfer in
atmospheres.  Then the corresponding integrals are evaluated for the
simple case of a spherically symmetric satellite, and applied to model
the perturbation due to the Earth--reflected radiation flux on LAGEOS.
For this purpose, the optical behaviour of the Earth's surface and
atmosphere is described as a combination of Lambertian diffusion
(continents), partial specular reflection consistent with Fresnel law
(oceans) and anisotropic diffusion according to Chandrasekhar's
radiative transfer theory (clouds).
The in--plane Gauss components
$T$ and $S$ vs. mean anomaly are computed for a simple orbital
geometry and for different models of the Earth's optical properties.
A sensitive dependence is found on the assumed cloud distribution,
with significant perturbations possibly arising from oceanic specular
reflection when the satellite is close to the Earth's shadow
boundaries. \\
\keywords
Celestial mechanics --- Artificial
satellites --- Non--gravitational perturbations --- Radiation
mechanisms.
\end{abstract}

\section{Introduction}
Non--gravitational perturbations on artificial satellites have been an
active research subject for celestial mechanicians in the 50s and the
early 60s, at the beginning of the space age; later on, the widespread
use of purely numerical models for predicting and fitting the orbital
evolution of satellites and probes decreased the interest in more
general and/or abstract approaches.  However, starting in the late 70s
the situation has changed again.  The main reason has been the
availability of very accurate tracking techniques, specifically
designed nearly spherical and ``passive'' satellites and drag--free
technology --- to be applied in particular to space geodesy and
geodynamics, but also to other purposes, such as general relativity
tests.  In this context, new physical mechanisms have been analyzed
and new formalisms and algorithms have been developed, both for
predicting the instantaneous value of very tiny non--gravitational
forces, and for modelling their long--term orbital effects (for a
general review, see Milani {\it et al.}, 1987).

One of the most complex such forces is that due to radiation pressure
resulting from light rediffused or reflected from the Earth's surface
and atmosphere --- commonly referred to as {\it albedo force} or {\it
albedo perturbation}.  Though the corresponding force is typically at
least one order of magnitude smaller than that due to solar radiation
pressure, the model uncertainties are normally higher, due to the
difficulties involved in modelling the momentum flux associated with
radiation for a source such as the Earth, geometrically extended and
having a complex and variable optical behaviour.  Often in the past
simplifying approximations have been used --- e.g., uniform albedo,
purely diffusive behaviour, etc. --- which neglect some critical
features of the process resulting into detectable orbital
perturbations.  For this reason, in the last decade a significant
number of papers have been devoted to discuss more general or refined
treatments of the albedo force, with applications to spacecraft as
diverse as LAGEOS, ERS--1, Magellan and MACEK (Anselmo {\it et al.},
1983; Barlier {\it et al.}, 1986; Rubincam and Weiss, 1986; Rubincam
{\it et al.}, 1987; Borderies, 1990; Borderies and Longaretti, 1990;
Lucchesi and Farinella, 1992; Vokrouhlick\'y and Sehnal, 1992a,
1992b).  Still, however, these models are typically based on
``paradigms'' on the relevant physical mechanisms which limit their
suitability in a more general framework.

The main purposes of this paper are the following~: (i) to present very
general formulae, inspired by the treatment of radiative transfer in
planetary and stellar atmospheres, that give the force due to
radiation pressure on a satellite of arbitrary shape moving within the
radiative field of an arbitrary source, and that can then be
specialized to deal with the simple case of a spherically symmetric
satellite (Section 2); (ii) to discuss a physical model of the optical
behaviour of the Earth's surface and atmosphere, which is relatively
simple from a physical point of view but at the same time realistic
enough to incorporate some subtle critical effects (Section 3); (iii)
to apply these formulae and models to LAGEOS, a satellite for which
the assumption about spherical symmetry provides a good approximation
of its real structure and whose laser--tracked orbit has shown
residual unmodelled accelerations probably caused for a significant
part by the albedo force (Section 4).

\section{Acceleration of a satellite plunged into a radiative field}
In this section we will formulate in a general manner
the mathematical tools needed for the derivation of the perturbative
acceleration of a satellite, arising from the interaction with a
radiative field.  Of course, this is not the first attempt to deal
with this problem (for another recent treatment see e.g.  Borderies
and Longaretti, 1990), but our aim is to present the theory in a
different formalism.  Our approach --- which was already partially
used by Vokrouhlick\'y and Sehnal, 1992a,b) is inspired by the
techniques widely used in the radiative transfer theory applied to
planetary and stellar atmospheres (see e.g.  Mihalas, 1978).

Our basic physical quantity is {\it the radiative intensity} $I$,
which from the macroscopic point of view fully describes the radiative
field interacting with the satellite.  More precisely, this is true
provided one neglects the effects of polarization, which are
negligible as far as the radiative acceleration of satellites is
concerned.  Moreover, we use the frequency integrated radiative
intensity, as we neglect any frequency dependence of all the relevant
processes.  This is actually a strong approximation, which {\it a
priori} rules out the possibility of taking into account energy
transfers due to frequency--dependent processes in the Earth's surface
and atmosphere and thus, for instance, of describing by a unified
theory the albedo effect (due to visible light) and the perturbations
due to the Earth's thermal radiation.  This can be justified {\it a
posteriori} by considering that the dynamical effects of the infrared
radiation field are very different, owing to its nearly--isotropic
geometry with no significant day--night asymmetry (see e.g.  Rubincam,
1987).

Many previous papers devoted to this subject use from the very
beginning the {\it radiative flux}, which in our approach is a
quantity derived from the radiative intensity (see Eq.~(1) below).
In our opinion, a clear formulation of the basic equations using as a
fundamental quantity the radiative intensity is useful, in particular
when deriving more complicated formulae (e.g. in the case of a complex
satellite geometry or of a radiative field without any symmetry).  The
formulation of Borderies and Longaretti (1990) is close to our
approach (note that they speak about {\it the radiance} rather than
about radiative intensity; we adopt the latter terminology which is
closer to that of the classical radiative transfer theory).  They
concentrate mainly on some subtle methods aiming at integrating their
formulae analytically as far as possible (for another impressive
analytical work see Rubincam and Weiss, 1986).  Our goal here is
different --- to present a formulation (i) as compact and clear as
possible; (ii) putting no restriction on the involved radiative field;
and (iii) allowing to incorporate directly the results of the radiative
transfer theory of the stellar/planetary atmospheres.  Of course, when
the equations have to be applied to specific concrete situations,
we will need to adopt suitable numerical integration methods; on the
other hand, the analytical formalisms such as those quoted above can
be applied only after strongly restricting the generality of the model
(e.g., by assuming purely diffusive light scattering by the Earth).

Let us assume that the satellite surface is a two--dimensional surface
${\cal S}$ parameterized by any two variables $\mu, \varphi$ (for
instance, for a spherical satellite we can use $\mu = \cos \vartheta$,
$\vartheta,\varphi$ being some spherical angles).  Moreover, we assume
that the satellite is plunged into a radiative field ${\cal R}$
parameterized by another pair of variables $\nu,\phi$ (for instance,
specifying a direction; in any specific situation both sets of
parameters together with the corresponding reference frame will be
explicitly given).  Under this assumptions, the perturbative
acceleration is given by the following formula 
\begin{equation}
 \vec{a} = \frac{1}{mc} \int_{{\cal S}(\mu,\varphi)} d
 A(\mu,\varphi)\! \left[ \left( {\bf E}\! +\! {\rho}(\mu,\varphi){\bf
 R}(\vec{N}(\mu,\varphi)) \right)\cdot\vec{J}_{1} + \frac{2}{3}
 {\delta}(\mu,\varphi) \vec{J}_{2}\right]\; ,
\end{equation}

\noindent where
$$ \vec{J}_{1} = - \int_{{\cal R}(\nu,\phi)} 
\Pi_{\vec{n}(\nu,\phi)} (\vec{N}(\mu,\varphi))\,
I(\nu,\phi;\mu,\varphi)\,\tilde{\Theta}(\nu,\phi;\mu,\varphi) \, d\nu d\phi\; ,$$
$$\vec{J}_{2} = - \int_{{\cal R}(\nu,\phi)}
\Pi_{\vec{N}(\mu,\varphi)} (\vec{n}(\nu,\phi))\,
I(\nu,\phi;\mu,\varphi)\,\tilde{\Theta}(\nu,\phi;\mu,\varphi) \, d\nu d\phi\; ,$$
$$ {\bf R}(\vec{a}) = 2 \Pi_{\vec{a}} - {\bf E} \quad , \quad
\tilde{\Theta}(\nu,\phi;\mu,\varphi) =\Theta
(-\vec{N}(\mu,\varphi)\cdot\vec{n}(\nu,\phi)) \; .$$

\noindent Here $d A(\mu,\varphi)$ is the satellite surface element,
{\bf E} is the unit matrix of dimension 3, $\Pi$ is a 
projection operator defined by~: $\Pi_{\vec{a}} (\vec{b}) =
(\vec{a}\cdot\vec{b})\; \vec{a}$, $\Theta(x)$ is the Heaviside step
function, $c$ is the velocity of light and $m$ is the satellite mass,
$\vec{N}(\mu,\varphi)$ is the surface normal and $\vec{n}(\nu,\phi)$
the unit vector in the the direction of a chosen ray.
The three terms in Eq.~(1) describe respectively the linear momentum
of radiation removed from the incident beam (both by
absorption and scattering),
specular reflection (phenomenologically described by
the coefficient $\rho(\mu,\varphi)$) and diffusive scattering by the
satellite's surface (phenomenologically described by the coefficient
$\delta(\mu,\varphi)$).  We stress that this decomposition of the
interaction of radiation with the satellite surface is just another
{\it paradigmatic} assumption of the radiation pressure model, whose
validity may become questionable under specific circumstances. 
Let us also note, that we take into account only processes
connected with direct momentum exchange between radiation and satellite,
disregarding the emission of {\it thermal}
radiation from it.  The latter effect is due to the true absorption
(i.e. contingent on energy exchange between the radiation and satellite).
Although the two effects are physically related, their consequences are
different~: while the reflected radiation acts {\it immediately} and
(roughly speaking) in a direction opposite to that of the source, the
thermal re-emission of the absorbed radiation is normally subjected to
some time delay (the so--called {\it thermal inertia}, which for a
spinning satellite may show up as a phase delay) and also depends on
other parameters characterizing the structure and orientation of the
satellite itself.  This effect has been studied in detail for LAGEOS
(see Rubincam, 1987; Afonso {\it et al.}, 1989; Farinella {\it et
al.}, 1990; Scharroo {\it et al.}, 1991) and we will not deal with it
further in the present paper.

It is to be stressed that formula (1) does not take into account possible
role of the photon multiple reflections on the satellite surface.
As for the specular reflection, this rather complicated task was
treated by tricky numerical models (e.g. Renard and Koeck, 1989).
However, keeping in mind the LAGEOS application, the satellite shape
is simple enough (spherical in a good approximation) so that 
formula (1) is sufficiently general. Moreover, since no shadowing
of one satellite part by another can occur, we suppress the
satellite surface element ($\mu,\, \varphi$) dependence of the 
radiative intensity $I$ (see formulas (2-5)).

We do not give here a full formal proof of Eq.~(1), for which we refer
to Vokrouhlick\'y, 1991.  The essential hints can also be found in
Vokrouhlick\'y and Sehnal, 1992a, who also give a proof of the full
equivalence with the formulae of Borderies (1990).  Instead, we will
apply Eq.~(1) to evaluate the total acceleration in a simple
particular case, that of a spherical satellite with constant optical
parameters ($\rho, \; \delta$).  Eq.~(1) involves in fact a double
integration (say, ``over the satellite surface'' and ``over the
radiation field'').  We shall show that in the case of a spherical
satellite the integral over the satellite surface can be evaluated
analytically, leaving just one integral in the final formula.  The
introduction of a spherical angle parameterization of the satellite
surface is advantageous in this case, though we shall fix later on the
reference frame (in fact we will profit from this freedom to simplify
the computations).

For practical reasons let us treat the three terms in Eq.~(1)
separately.  The first one (${\bf E}$ in the integrand brackets)
is to be treated in the following
way~:\footnote{Note, that we implicitly neglect variations of
the radiative field at different part of the satellite.
Instead, we work with radiative field $I(\nu,\phi)$, which results
from averaging over its surface. The same applies for formulas (3) and (4).}
\begin{eqnarray}
\vec{a}_{A} & = & - \frac{R^2_{\cal L}}{mc} \int_{{\cal
 S}(\mu,\varphi)} d\mu d\varphi
 \int_{{\cal R}(\nu,\phi)} d\nu d\phi \,
\left(\vec{n}(\nu,\phi)\cdot\vec{N}(\mu,
 \varphi) \right) \vec{n}(\nu,\phi)
 I(\nu,\phi)\, \nonumber \\
 &  &\qquad \qquad \qquad \qquad \qquad \qquad  
 \tilde{\Theta}(\nu,\phi;\mu,\varphi) \nonumber \\
 & = & - \frac{R^2_{\cal L}}{mc} \int_{{\cal R}(\nu,\phi)} d\nu d\phi \,
 I(\nu,\phi) \, \vec{n}(\nu,\phi)
 \int_{{\cal S}(\mu,\varphi)} d\mu d\varphi
\left(\vec{n}(\nu,\phi)\cdot
 \vec{N}(\mu,\varphi)\right) \nonumber \\
 &  &\qquad \qquad \qquad \qquad \qquad \qquad  
 \qquad \qquad \tilde{\Theta}(\nu,\phi;\mu,\varphi) \nonumber \\
 & = & \frac{R^2_{\cal L}}{mc} \int_{{\cal R}(\nu,\phi)} d\nu d\phi \,
 I(\nu,\phi) \, \vec{n}(\nu,\phi)
 \int_{\mu=0}^{1} \int_{\varphi=0}^{2\pi} d\mu d\varphi \mu
 \nonumber \\
  & =  & \frac{A_{\perp}}{mc} \int_{{\cal R}(\nu,\phi)} d\nu d\phi \,
 I(\vec{n}(\nu,\phi))\, \vec{n}(\nu,\phi) \; .
\end{eqnarray}

\noindent where $R_{\cal L}$ and $A_{\perp}$ are the satellite radius
and cross section, respectively.  Note that we have fixed the local
satellite frame (in which the spherical angles $\vartheta,\; \varphi$
are defined) in such a way to simplify the expression in the second
line of Eq.~(2), thus resulting in the third line.  The quantities
$\nu,\; \phi$ parameterizing the radiative field (up to now
unspecified, but normally also spherical angles in possibly another
local frame) specify 'one ray' in it.  The second line in Eq.~(2)
clearly shows that for every such ray we have to fulfill the
integration over the satellite surface ${\cal S}$. In each of these
cases we have rotated the satellite local frame in such a way that its
$z$--axis coincides with the direction of the chosen ray.  Thus, for
each ray ($\nu,\; \phi$) we use a different satellite frame (and
therefore different spherical angles $\vartheta,\; \varphi$).  This
trick is allowed by the isotropy of satellite surface with respect to
the parameters $\rho, \; \delta$.

In the same manner we can treat the other two terms in Eq.~(1).
The second one, corresponding to the specularly reflected radiation,
becomes
\begin{eqnarray}
\vec{a}_{R} & = & - \frac{R^2_{\cal L}}{mc} \int_{{\cal
 S}(\mu,\varphi)} d\mu d\varphi \int_{{\cal R}(\nu,\phi)} d\nu d\phi
 \left(\vec{n}(\nu,\phi)\cdot\vec{N}(\mu, \varphi) \right)\,I(\nu,\phi)
 \nonumber \\
 &  &\qquad \qquad \qquad \qquad \qquad \qquad  \left({\bf R}(\vec{n}
 (\mu,\varphi))\cdot\vec{n}(\nu,\phi) \right)
\,\tilde{\Theta} (\nu,\phi;\mu,\varphi) \nonumber \\
 & = & - \frac{R^2_{\cal L}}{mc} \int_{{\cal R}(\nu,\phi)} d\nu d\phi
 I(\nu,\phi)
 \int_{{\cal S}(\mu,\varphi)} d\mu d\varphi
\left(\vec{n}(\nu,\phi)\cdot
 \vec{N}(\mu,\varphi)\right) \nonumber \\
 &  &\qquad \qquad \qquad \qquad \left(2 \Pi_{\vec{N}(\mu,\varphi)}
 (\vec{n}(\nu,\phi)) - \vec{n}(\nu,\phi) \right) \tilde{\Theta}(\nu,\phi;
 \mu,\varphi) \nonumber \\
 & =  & \frac{R^2_{\cal L}}{mc} \int_{{\cal R}(\nu,\phi)} d\nu d\phi
 I(\nu,\phi) \vec{n}(\nu,\phi)
 \int_{\mu=0}^{1} \int_{\varphi=0}^{2\pi} d\mu d\varphi \mu (2\mu^2 -1)
 =  0 \;.
\end{eqnarray}

This result is a classical one, and in fact has its counterpart in the
same result for the elastic specular reflection of the neutral atmosphere
molecules in the theory of the atmospheric drag.  Of course any,
but very particular cases,
deviation from the spherical shape or from the isotropy of the
specular reflection coefficient $\rho$ would lead to a nonzero result.

As for the third term (diffusively scattered radiation from the
satellite surface), we have
\begin{eqnarray}
\vec{a}_{D} & = & - \frac{2\delta R^2_{\cal L}}
{3cm}\int_{{\cal S} (\mu,\varphi)} d\mu d\varphi \int_{{\cal
R}(\nu,\phi)} d\nu d\phi \left(\vec{n}(\nu,\phi)\cdot\vec{N}(\mu,
\varphi) \right) \, \vec{N}(\mu,\varphi)\, I(\nu,\phi)\, \nonumber \\
 &  &\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad  
 \tilde{\Theta}(\nu,\phi;\mu,\varphi) \nonumber \\
& = &- \frac{2\delta R^2_{\cal L}}{3cm} \int_{{\cal R}(\nu,\phi)} d\nu
d\phi \, I(\nu,\phi) \int_{{\cal S}(\mu,\varphi)} d\mu d\varphi
\left(\vec{n}(\nu,\phi)\cdot \vec{N}(\mu,\varphi)\right)
\vec{N}(\mu,\varphi) \nonumber \\
 &  &\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad  
 \tilde{\Theta}(\nu,\phi;\mu,\varphi) \nonumber \\
& = & \frac{2\delta R^2_{\cal L}}{3cm} \int_{{\cal R}(\nu,\phi)} d\nu
d\phi I(\nu,\phi) \vec{n}(\nu,\phi) \int_{\mu=0}^{1}
\int_{\varphi=0}^{2\pi} d\mu d\varphi \mu^2 \nonumber \\
& = & \frac{4\delta A_{\perp}}{9cm} \int_{{\cal R}(\nu,\phi)} d\nu d\phi
I(\vec{n}(\nu,\phi)) \vec{n}(\nu,\phi)\; .
\end{eqnarray}

Here note one subtle point -- the arguments of the vectorial quantity
have changed proceeding from the second line to the third.  This is
also luckily due to the spherical symmetry of the satellite surface
(roughly speaking, this means that the integral of the normal vectors
to the sphere over one hemisphere is parallel to its axis).  In the
case of nonspherical shapes, the result would thus be different.

The resulting perturbing acceleration is $\vec{a} = \vec{a}_{A} +
\vec{a}_{R} + \vec{a}_{D}$, hence in total
\begin{eqnarray}
\vec{a} & = & \frac{A_{\perp}}{mc} \left(1 + \frac{4}{9} \delta \right)
 \int_{{\cal R}(\nu,\phi)} d\nu d\phi \vec{n}(\nu,\phi) I(\vec{n}(\nu,\phi))
 = \frac{A_{\perp}}{mc} \left(1 + \frac{4}{9} \delta\right) \vec{F}
 \nonumber \\
& = & \left(\begin{array}{c} {\rm satellite} \\ {
 \rm cross\, section} \end{array}
 \right)
 \left(\begin{array}{c} {\rm some\, coefficient\, characterizing\,
  the\, inte-} \\
 {\rm raction\, of\, light\, with\, the\, satellite\, surface}
 \end{array}\right) \nonumber \\
& & \left(\begin{array}{c} {\rm radiative\, flux,\, i.e.\,
 energy\, transferred} \\ {\rm through\, the\, perpendicular\,
 area\, unit} \end{array} \right)
 \left(\begin{array}{c} {\rm velocity} \\ {\rm of\,
 light} \end{array}\right)^{-1}
 \left(\begin{array}{c} {\rm satellite} \\ {
 \rm mass} \end{array}\right)^{-1}
\end{eqnarray}

The second and third lines give a clear ``physical interpretation''
of the relevant terms (notice that the product of the two terms in
third line yields the momentum flux, namely the amount of linear
momentum transferred through a perpendicular unit area).  This result
of course is not new; see e.g.  Aksenov, 1977, Milani {\it et al.},
1987 (pp. 74--75), Borderies and Longaretti, 1990 (where one has to
substitute $\nu \rightarrow \frac{1}{3} \delta$); here it is just an
example of the application of the more general Eq.~(1).

The result given by Eq.~(5) in general still holds for
nearly--spherical satellites with {\it the quasi--isotropic} surface
optical properties.  This --- admittedly, somewhat vague --- notion is
inspired by LAGEOS--type satellites, whose surface is covered by a
large number of retroreflectors.  We call the surface of a
nearly--spherical satellite to be quasi--isotropic, if, roughly
speaking, ``it looks the same from all the directions to a high
accuracy''.  For instance, let us imagine a parallel beam of radiation
impinging on LAGEOS.  Obviously, if the satellite is illuminated along
the axis of any retroreflector the radiative force will be somewhat
different with respect to the case of irradiation out of the
retroreflector axis.  But the difference will be very small, provided
the retroreflector has a small size compared with the satellite.  The
quasi--isotropy assumption just neglects this small difference.  More
precisely, we have to assume that the results of the satellite surface
integration along the different rays of the radiative field (as
explained in the paragraph following Eq.~(2)) in the formulae (3), (4)
and (5) are the same, or at least that the differences are negligible.

Notice that the peculiar optical behaviour of the retroreflectors ---
perfect backscattering of the impinging radiation --- causes their
contribution to the momentum transfer equation to be just twice that
of the absorption term (2), so that for a nearly--spherical and
quasi--isotropic satellite whose surface is covered by retroreflectors
for a fraction $f$ of the total area, the factor $\left(1+\frac{4}
{9}\delta\right)$ in Eq.~(5) needs just to be replaced by
$\left[1+f+\frac{4}{9}\delta(1-f)\right]$.  For LAGEOS, the effective
value of $f$ is much lower than the ``geometrical'' one ($\approx
0.4$), because the corner--cube retroreflectors backscatter light only
within some restricted incidence cone about their axis.  More in
general, the classical factor $\left(1 + \frac{4}{9} \delta\right)$
will be replaced by a coefficient ${\cal C}_{R}(\rho,\delta)$, where
the assumption of quasi--isotropy assures that ${\cal
C}_{R}(\rho,\delta)$ is a function of the surface coefficients $\rho,
\delta$ only and not e.g. of the orientation of the satellite rotation
axis.  The acceleration is thus simply rescaled if, for instance, the
optical parameters are not constant in time.  For LAGEOS, ${\cal C}_R$
has been treated as a constant solve--for parameter in the orbit
determination process, and its value has been found to be about $1.2$.

Of course, if the quasi--isotropy assumption fails by a significant
amount, the perturbative acceleration will be given by an expression
more complicated than (5).  Interestingly, the possibility of a
sizeable difference between the optical properties of two different
hemispheres of LAGEOS has been recently proposed by Scharroo {\it et
al.} (1991) as a way of explaining some unmodelled peaks in LAGEOS'
residual along--track acceleration.  In the present context, however,
we shall stick to the quasi--isotropy assumption.

Finally, note that the character of the radiative field is not limited
in any way in the previous discussion.  It may be either the nearly
parallel direct solar radiation or the geometrically more complex
radiative field coming from the Earth's surface and atmosphere.  In
what follows, we shall apply Eq.~(5) to the latter case.  In the next
section we shall describe a model of this radiative field as well as
its $(\nu,\phi)$ parameterization in Eq.~(5).

\section{Models of the Earth's radiative field}
We shall now apply the preceding formulae to the
particular case when the radiative field is originated by the
reflection/diffusion of the sunlight from the Earth --- the so--called
{\it albedo effect}.

The Earth is assumed to be spherical.  No general {\it a priori}
assumption is made on the optical characteristics of the Earth
surface, but they are just subject to a set of different model
assumptions, with the aim of studying their influence on the resulting
orbital perturbations of LAGEOS.  In particular, we are going to
consider the specific role of the scattering of light from the cloud
coverage.  Actually, one should always speak about the
reflection/diffusion from the Earth--atmosphere system rather than
from the Earth surface only.

Also, we assume that the Earth--scattered radiation does not interact
with any exoatmospheric near--Earth medium, so that the intensity of
the radiation is constant along the photon rays out of the atmosphere
(Mihalas, 1978).  This assumption is essential to obtain a geometrical
link between the Earth's optical models (treated basically in the
Earth surface/upper atmosphere local frames) and the radiative field
which appears in the satellite acceleration formulae (1) or (5).
Recall that in the previous Section we considered processes occurring
in a satellite--fixed local frame, while throughout this section we
will formulate the Earth's optical models in the Earth
surface/atmosphere local frames.  The previous assumption bridges the
gap between them.

Similarly to Vokrouhlick\'y and Sehnal (1992b), we will consider three
different {\it modes} of re\-flec\-tion/diffu\-sion from the Earth's
sur\-face/atmo\-sphere.  They are~:
\begin{itemize}
\item diffusion of light from the continents, which will be
assumed to occur in a simple isotropic way ({\it Lambert diffusion})
with a multiplying albedo value developed in spherical harmonics of
the geographical latitude and longitude (as given by Sehnal, 1979);
\item reflection/diffusion from the oceans~: here together
with a diffusive behaviour part of the light is assumed to be reflected
in a specular way;
\item scattering of light from compact atmosphere
formations (which we will simply name {\it clouds}), treated according
to Chandrasekhar's solutions of the radiative transfer equations for
a planar atmosphere.  The detailed assumptions will be specified
in what follows.
\end{itemize}

\begin{figure}
\vspace{76mm}
\caption[]{The spherical coordinate system centered at
an Earth surface element.  The angle $\vartheta$ is measures the
distance from the zenith direction of the chosen element.  The
origin of the $\varphi$ coordinate is irrelevant in our models.}
\end{figure}
We introduce the spherical angles $\vartheta, \; \varphi$ (and also
$\mu = \cos \vartheta$) in the local Earth surface/atmosphere frame
(these angles should not be confused with those used in the
previous section).  The $z$--axis of this local frame coincides with
the outer normal of the surface element.  The orientation of the other
two axes does not appear in the formulae and consequently it is
irrelevant.  We reserve $\mu_{\odot}, \; \varphi_{\odot}$ for the
Solar direction, while $\mu,\; \varphi$ specify the satellite
direction in this system (see Fig.~1).

Using these notations, one can write the intensity of the
isotropically diffused radiation according to Lambert's law as~:

\begin{equation}
I^{({\rm cont})}_{{\rm dif}} = \frac{1}{\pi} \Phi_0 {\cal A}
\mu_{\odot}\; , 
\end{equation}

\noindent where $\Phi_0$ is the flux (${\rm W/m}^2$) of the incoming 
solar radiation
and ${\cal A}$ is {\it the albedo} coefficient.  The $\mu_{\odot}$
factor is just due to geometrical reasons, namely is proportional to
the cross--section with respect to the solar rays.  Some criticisms to
this purely geometrical approach can be found as early as in Levin,
1962, whose arguments were explicitly quoted by Vokrouhlick\'y and
Sehnal (1992b).  Later on, the deviations from linearity with respect
to $\mu_{\odot}$ were referred to as dependence of the albedo ${\cal
A}$ on $\mu_{\odot}$ (see e.g.  L\'ala et al., 1978; Rubincam et al.,
1987; Knocke, 1989; and from an observational point of view Stephens
et al., 1981).  We feel that in spite of these studies, which usually
empirically fit some sets of observed data, the theoretical
understanding of this phenomenon is not satisfactory.  Both light
scattering in the atmosphere and details of the diffusion mechanisms
on the Earth's solid surface (e.g., small--scale mutual shadowing) are
plausible causes for a $\mu_{\odot}$ dependence of the albedo, and
more in general for a dependence on all the involved angular
parameters ($\mu_{\odot},\varphi_{\odot}\ ;\mu,\varphi$).  But little
is known about the separate importance of each of these effects and
their detailed physics.  In view of this, we have decided just to use
the normal albedo model as given by Sehnal (1979).  The albedo ${\cal
A}$ is here a general function of the geographical latitude and
longitude (expanded in spherical harmonics), but is constant with
respect to the solar zenith angle $\mu_{\odot}$.  In fact, it may be
understood as an {\it averaged} value over $\mu_{\odot}$.

As for the oceans, we assume that besides Lambertian diffusion their
surface can cause specular reflection of light.  Following Barlier
{\it et al.} (1986), we then decompose the intensity of the scattered
radiation in two parts~:
\begin{eqnarray}
\vspace{2mm}
I^{({\rm oc})}_{{\rm dif}} & = & \frac{1}{\pi} \Phi_0 {\cal A} \mu_{\odot}
(1-{\cal R}(\mu_{\odot}))\; , \nonumber \\
I^{({\rm oc})}_{{\rm spec}} & = & \Phi_0 {\cal A} {\cal R}(\mu_{\odot})
\psi_{\rm sc}(\mu_{\odot}) \delta(\mu-\mu_{\odot}) \delta
(\phi-\phi_{\odot}+\pi) \; ,
\end{eqnarray}

\noindent where
\begin{eqnarray}
\vspace{2mm}
\psi_{sc}(\vartheta_{\odot},\alpha ) & = &
\frac{\alpha^2\sin\vartheta_{\odot}
\cos\vartheta_{\odot}}{\sin\beta\left(2\cos(2\vartheta_{\odot}-
\beta )-\alpha \cos\vartheta_{\odot}\right)}\, , \\
\alpha\ & = & R_{\oplus}/r\; , \nonumber 
\end{eqnarray}

\noindent $r$ being the current satellite geocentric distance and the
angle $\beta$ fulfills the equation
$$\alpha\ \sin\vartheta_{\odot}\ = \sin(2\vartheta_{\odot}-\beta )\; .
$$

\noindent Note that, in agreement with the definition of ${\cal A}$
as a fraction of incident sunlight which is not absorbed, but
either diffused or reflected, the factor ${\cal A}$ appears also in
the intensity of the specularly reflected component.
Barlier {\it et~al.} (1986; see their Eq.~(10))on the contrary
included ${\cal A}$ only in the intensity of the diffused component,
thus implicitly adopting a different definition of albedo and
amplifying the effects of radiation pressure from specularly
reflected sunlight when compared with the present work.

For the derivation of Eq.~(8) see Barlier {\it et al.},
1986, where this factor is denoted as $J(\theta )$ (it actually
simulates the divergence of the reflected radiation field due to the
fact that the assumed mirror surface is spherical).  The $\delta$
function in Eq.~(7) is an obvious consequence of the specular
reflection assumption.  Of course, this is a somewhat idealized model~:
small scale wave phenomena on the ocean surfaces will actually cause
some more divergence of the reflected light (as with the ``sword of
the sun'' observed from a beach at sunset).  A phenomenological model
of this beam divergence can be found in Rubincam et~al., 1987, who
instead of the $\delta$ function adopt an exponential behaviour, and in
Lucchesi and Farinella, 1992, who assume that the light intensity is
uniformly distributed inside a narrow cone centered at
$\mu=\mu_{\odot}$.  The latter paper did show that the resulting
perturbations depend weakly on the aperture of the reflection cone (or
the corresponding exponential decay coefficient) provided it is less
than several degrees.  Therefore, here we shall stick to the simpler
specular reflection assumption.

The ``partitioning  function'' ${\cal R}$ in Eq.~(7) is derived
from the Fresnel reflection law (see Barlier {\it et al.}, 1986).
The implicit assumptions here include a sufficiently short
wavelength of the radiation and a negligible polarization (for a detailed
treatment consult e.g. Jackson, 1975, and Ditchburn, 1976).  Note that
the Fresnel function ${\cal R}$ fulfills the following relation

$$ \left(1-{\cal R}(\mu_{\odot}) \right)_{\mu_{\odot} \simeq 0} \simeq
\mu_{\odot} + {\cal O}(\mu^2_{\odot}) \; , $$

\noindent which results into a rapid increasing importance of the
specularly reflected part of radiation with respect to the diffusively
reflected part when the satellite is close to the Earth shadow's
boundary, namely it ``sees'' a narrow Earth crescent.  This is due to
the fact that when $\mu_{\odot} \simeq 0$, $I^{\rm (oc)}_{\rm dif}
\propto \mu^2_{\odot}$ but $I^{\rm (oc)}_{\rm spec} \propto
\mu_{\odot}$ only.  The corresponding sudden ``kick'' due to the
specularly reflected radiation is in a way the result of extreme
assumption, as the decomposition of the scattered radiation in the two
parts is not unique (see Vokrouhlick\'y and Sehnal, 1992b) and not all
such decompositions give this feature.  It should be noted that in
accepting this model we just follow another common argument made in
recent literature on LAGEOS perturbations~: that of testing the effects
of extreme assumptions whenever a reliable physical understanding is
missing.  In Sec. 4 we shall comment on some consequences of this
assumption when describing the results of our model.

As for the scattering of light from the clouds, we shall use the model
developed by Chandrasekhar (see e.g.  Chandrasekhar, 1950; van de
Hulst, 1980).  The intensity of the scattered radiation is given by
\begin{eqnarray}
I^{({\rm cl})}_{\rm scat} & = & \frac{1}{\pi} \frac{\omega_0\mu_{\odot} \Phi_0}
{\mu+\mu_{\odot}} (H(\mu;\omega_0,x)H(\mu_{\odot};\omega_0,x)
\Psi(\mu,\mu_{\odot};\omega_0,x)- \nonumber \\
 &  & x\sqrt{(1-{\mu}^2)(1-{\mu_{\odot}}^2)}
H^{(1)}(\mu;\omega_0,x) H^{(1)}(\mu_{\odot};\omega_0,x) \cos
(\varphi-{\varphi}_{\odot}))\; , \nonumber \\
 & & 
\end{eqnarray}

\noindent where
$$\Psi(\mu,\mu_{\odot};\omega_0,x) =
1-c({\omega}_0)(\mu+\mu_{\odot})-x (1-\omega_0)\mu\mu_{\odot} \; .$$

These formulae solve the radiative transfer problem in a planar
atmosphere with infinite optical thickness and scattering centers
characterized by the phase function
\begin{eqnarray}
p(\Theta) = {\omega}_0 (1+x \cos \Theta),\quad 0\leq \omega_0 \leq 1
,\quad -1\leq x \leq 1 \; ,
\end{eqnarray}

\noindent where $H(\mu)$ and $H^{(1)}(\mu)$ are the Chandrasekhar
H-functions and the associated $\mu$--moments.  They are defined as
solutions of integral equations, but several numerical techniques for
their evaluation are known (van de Hulst, 1980). $\Theta$ is the angle
between the incident and the scattered ray at one scattering center in
the atmosphere, and the phase function $p(\Theta)$ weighs the
scattering probability at this angle.

The solution (9) is two--parametric, with the parameters $\omega_0$
and $x$ controlling physically different effects~:

\begin{itemize}
\item $\omega_0$ can just be interpreted as the
albedo of the scattering centers in the atmosphere (not to be
confused with the global atmospheric albedo ${\cal A}_{eff}$;
their relation is discussed in Vokrouhlick\'y and Sehnal, 1992b, e.g.
in the case $x=0$ one gets ${\cal A}_{eff} = \frac{4\omega_0}
{(1+\omega_0)^2}$, while in the case $x \neq 0$ no such simple relation
exists),
\item the parameter $x$ is related to the asymmetry of the
light scattering from separate centers in the atmosphere~: $x=0$
implies symmetry, while $x=1$ and $x=-1$ correspond to the extreme
forward-- and back--scattering cases (see Eq.~(10)).  Since
we are interested in testing models with local asymmetry of the cloud
scattering behaviour (as it will be explained later), we will often put
$x=1$ (note however, that due to simplicity of the phase function (10), 
we cannot reach ``the extreme asymmetry'' as described by the asymmetry
parameter $g$ defined e.g in van de Hulst, 1980; for a review of more
realistic phase functions see the brilliant book van de Hulst, 1980 and
references in).
\end{itemize}

\noindent As with many theoretical models, that given by Eq.~(9)
suffers from several oversimplifications when compared to reality,
e.g.~: 
\begin{itemize}
\item infinite optical thickness of the atmosphere (i.e., clouds); 
\item simplified phase function;
\item both vertical and horizontal homogeneity of the atmosphere.
\end{itemize}

These simplifying assumptions are shortly commented in Vokrouhlick\'y
and Sehnal, 1992b, while more detailed discussions can be found in the
specialized literature (Chandrasekhar, 1950; Irvine, 1975; van de
Hulst, 1980).  Here we are not going to include possible
generalizations according any of the items listed above.  But we
stress that, as far as the first two items are concerned, these
generalizations can be introduced if the comparison with the data
shows that they are necessary to obtain a better fit.

We also note that all the modes of the light
scat\-te\-ring/diffu\-sion/re\-flec\-tion lead to formulae for the
resulting radiative intensity depending only on the angle $\varphi -
\varphi_{\odot}$ and not on $\varphi,\; \varphi_{\odot}$ independently
(this property justifies the choice not to define the orientation of
the local $x$--, $y$--axes, as noted earlier).  This is well
understandable, since the $(x,y)$ local plane was chosen as tangent to
the Earth's surface/atmosphere.  We do not have any reason to
distinguish symmetric directions with respect to the incident solar
radiation.  As far as we know, also all the models inferred from
empirical data (e.g.  Taylor and Stowe, 1984; Rubincam {\it et al.},
1987) have the same property.

We now turn to description of the geographical features of the Earth's
surfa\-ce/atmo\-sphere that we have used in the computation of the integrals in
Eq.(5).  In all our models, we used a realistic distribution of the
continents, modelled by a $2^o \times 2^o$ grid in latitude and
longitude.  Then the final Earth {\it mask} giving the morphology of
the different optical modes is created by covering the continental
grid with an atmospheric grid, which gives the assumed cloud
distribution.  Of course, the latter is rapidly variable, and a
realistic description could only be derived {\it a posteriori} from
the available knowledge of global meteorological processes.  This is
well beyond the scopes of this work.  Here we shall just distinguish
among several possible idealized cases~: 
\begin{itemize}
\item the atmospheric grid is empty (we shall refer to this model as the {\it no
cloud model}); 
\item the atmospheric grid is full ({\it full cloud model}); 
\item the atmospheric grid is empty on
the Southern hemisphere and full on the Northern one ({\it
hemisphere model}); 
\item the atmospheric grid is full
in a belt surrounding the equator up to $\pm 10^\circ$ in geographical
latitude and empty elsewhere ({\it cloud belt model});
\item $60\%$ of the atmospheric grid has been filled at random ({\it
stochastic cloud model}).
\end{itemize}

\section{Albedo perturbations on LAGEOS}
In this section we will discuss the perturbative
accelerations resulting from the different models introduced earlier,
and how they vary along one satellite revolution around the Earth.
For the sake of simplicity we shall use a simplified orbit geometry,
as already done by Rubincam {\it et al.} (1987) and Lucchesi and
Farinella (1992)~: LAGEOS' orbit is assumed to be polar ($I =
90^\circ$; the real orbit has an inclination of about $110^\circ$) and
the Sun lies in the LAGEOS orbit ascending node.
We keep the true LAGEOS eccentricity ($e=0.00444$) and 
longitude of perigee at epoch May 15, 1976 ($\omega=235.35^o$, see
Barlier {\it et~al.}, 1986), which determines the origin of
anomalies. This configuration will
allow us to estimate the maximum amplitude of the components of the
albedo effect perturbative accelerations in the radial ($S$) and
in--plane transverse ($T$) directions.  These two components are
important because they can be used to derive both the short-- and the
long--term perturbations on semimajor axis and mean anomaly; notice,
however, that the assumed geometry minimizes the out--of--plane
component $W$, which is important if one aims at estimating the
inclination and node perturbations.

The most time--consuming part of the numerical treatment of the
problem is the integration over the radiative field contributing to
the radiative flux at the satellite.  In order not to lose the
detailed information on the Earth's sur\-face/atmo\-sphere optical
behaviour, for each satellite position along the orbit we divided the
visible Earth cap into a grid of $3025$ surface elements, chosen in
such a way as to be of equal area.  This latter requirement implies that
the grid is somewhat denser near the edge of the visible cap, and this
may be important to properly model the acceleration when the satellite
is close to the Earth's shadow boundary and just ``sees'' a thin
illuminated Earth crescent.  A similar method was used by Lucchesi and
Farinella (1992), who showed that at least $\approx 100$ surface
elements are required for an accurate calculation; in our case the
presence of specular reflection makes this requirement even more
stringent.  Similar accuracy problems have been taken into account in
choosing the step size along the orbit; we have used $360$ steps per
revolution, i.e. $1^{\circ}$ mean anomaly steps.

\begin{figure}
\vspace{76mm}
\caption[]{The transverse component $T$ of the albedo
acceleration as a function of the mean anomaly $M$ over one LAGEOS
revolution, for the orbit geometry described in the text.  Different
{\it no cloud} models have been used to derive curves 1~: the (1a)
model includes the real distribution of the continents on the Earth
surface, the (1b) model considers an idealized Earth totally covered
by the oceans, and the (1c) model assumes only diffusive (continental)
reflection from the whole Earth surface.  Note the essential role of
the specular reflection effects when the satellite is close to the
boundaries of the Earth's shadow.  The dotted curves 2 bound the
results of $50$ different runs of the {\it stochastic} cloud model.
Curves 3 correspond to the full cloud model. the adopted cloud
parameters in both cases are~: $\omega_0 = 1, \; x=1$.}
\end{figure}
The results of several models are presented in Fig.~2, where the
transverse component $T$ is given as function of the mean anomaly $M$
(whose origin is taken in such a way that at the subsolar point
$M=100^\circ$).  The curve(s) labeled 1 correspond to the no cloud
model.  More precisely, the solid curve 1a corresponds to the no cloud
model as described in the Sec.~3, which takes into account the real
distribution of continents, while the dashed curves 1b and 1c
``ignore'' the distribution of the continents~: curve 1b corresponds to
the ideal case when the whole Earth is covered by an ocean with the
optical behaviour described in Sec.~3, and curve 1a corresponds to the
ideal case when the whole Earth surface is covered by a continent with
diffusive reflection.  As expected, the largest differences occur near
the satellite entry to/exit from the Earth shadow.  Curve 1c shows
there a smooth decrease, while curve 1b exhibits clearly the specular
reflection ``flash'', as discussed earlier.  The no cloud model 1a,
combining 1b and 1c according to the actual distribution of the
continents, shows sudden jumps due to interchanging of the two ``pure''
modes of reflection (due to the chosen orbit geometry, in the enlarged
window in Fig.~2 the zone of the Earth possibly relevant for
reflection is the Arctis).  Such structures have a significant
influence on the averaged value $<T>$ over an entire orbit, which is
important to yield long--term perturbations.  Such a potentially
important role of specular reflection from relatively small zones of
the Earth's surface was already pointed out by Barlier {\it et al.}
(1986). It may be noted that our assumption that the Arctic polar cap 
behaves as an ocean is somewhat artificial; however, according to
the data reported by Taylor and Stowe (1984), snow--covered
surfaces do generate some specular reflection features, albeit less
marked than liquid water, and thus may be seen as ``intermediate''
between the oceans and the continents.

Curve 3 corresponds to the full cloud model with parameters $(\omega_0
= 1,\; x=1)$.  As no specular reflection occurs in this case, the
curve is smooth.  Actually, as noted already by Vokrouhlick\'y and
Sehnal (1992b), this case closely resembles the simple one with just a
diffusion with a constant effective albedo ${\cal A}_{eff} =
1$ (see Lucchesi and Farinella, 1992, Fig. 5a, where the value ${\cal
A}_{eff} = 0.2$ was used); only the anisotropy of the local scattering
behaviour slightly disturbs the diffusion pattern.  The amplitude of
curve 1 is thus consistent with an average value of the albedo for the
Earth's surface of approximately $0.3$.

The region bounded by the dotted curves 2 corresponds to the
stochastic cloud model; in this case, there is no unique curve,
because the filling of the atmosphere grid is made at random.
Therefore, we performed $50$ runs of the program computing $50$
different cloud distributions (always with a filling fraction of
$60\%$), and then plotted the maximum and minimum values of the
along-track accelerations.  The typical difference of $10^{-11}\ m/s$
gives an idea of the consequences of the variability in the real cloud
distribution (note in particular that the occurrence of the ``Arctic
reflection flash'' is subject to the local cloud coverage).  A similar
method was used the case of the lower--altitude satellite ERS--1 by
Vokrouhlick\'y and Sehnal (1992b).  An interesting observation can be
done by comparing the ERS--1 results with those in Fig.~2 --- the
three curves 1--3 in Fig.~2 are well separated, as the region between
the dotted curves 2 does not span the whole ``gap'' between curves 1
and 3, while the opposite was true for ERS--1.  The reason has to do
both with the simpler satellite geometry and with the higher orbit of
LAGEOS.

\begin{figure}
\vspace{76mm}
\caption[]{The same as in Fig.~2, but comparing
results of the {\it no cloud} model (curve 1), the {\it
hemisphere} model (curve 2), and the {\it full cloud} model
(curve 3).  Cloud model parameters are kept as follows~: $\omega_0 =
1, \; x=1$.}
\end{figure}
In Fig.~3 we have replotted the mean anomaly dependence of the $T$
acceleration component in the case of the no cloud model (curve 1),
and the full cloud model (dashed curve 2), but we have added the curve
(3) corresponding to the hemisphere model.  As expected, the
hemisphere model signal coincides with either one of the two previous
cases when the satellite is over the poles (as its ``visible cap'' of
the Earth behaves according to only one mode of reflection), while
differs from them in the transient region, when the satellite is over
the equator.

Curves 3 in Fig.~4 compare with the no cloud (curve 1) and full cloud
(curve 2) models the results of the cloud belt model.  The dashed and
dotted curves allow an estimate of the effects of changing the $x$
isotropy parameter~: dashed curves correspond to the $x=1$ choice,
dotted curves to the $x=0$ one.  No qualitative new feature appears,
though there is a significant quantitative difference.  However, this
is not a general result, due to the assumed symmetrical orientation of
the orbital plane with respect to the Sun. Curves 3a--c also show the
dependence of the results on the effective cloud albedo~: 3a
corresponds to ${\cal A}_{eff} = 1$, 3b to ${\cal A}_{eff} = 0.8$, and
3c to ${\cal A}_{eff} = 0.6$.  A decreasing effective albedo of the
clouds clearly diminishes their role.
\begin{figure}
\vspace{76mm}
\caption[]{Here the results of the {\it no cloud}
model (solid curve 1) and of the {\it full cloud} models (curves 2)
are compared with those of the {\it cloud belt} model (curves 3).
Dashed curves correspond to strongly anisotropic light scattering from
the atmosphere centers ($x=1$), dotted curves to symmetric scattering
($x=0$).  The dependence on the effective cloud albedo is tested as
well~: curves (3a) correspond to ${\cal A}_{eff} = 1$, curves (3b) to
${\cal A}_{eff} = 0.8$, and curves (3c) to ${\cal A}_{eff} = 0.6$.}
\end{figure}

Other remarks follow from the inspection of the radial component of
the perturbative acceleration.  Its dependence on the mean anomaly is
shown in Fig.~5 in the case for the some models discussed in the
previous paragraphs.  Curve 1 corresponds to the no cloud model,
curves 2 to the full cloud model (with different values of $x$) and
curves 3 to the cloud belt model (with different values of both $x$
and ${\cal A}_{eff}$).  Similar graphs can be plotted for the other
models.  The observed differences are clearly of the order of several
times $10^{-10}\ m/s^{2}$.  As the real cloud distribution is highly
variable and indeed unknown when analysing the real tracking data
(unless a special effort is made to recover global meteorological
data), one can expect an intrinsic uncertainty on the magnitude of the
instantaneous albedo acceleration for LAGEOS of the order of
$10^{-10}\ m/s^2$.  The corresponding error in LAGEOS' position
predictions (or determinations) over a fraction of one orbital period
($\approx 4\ hr$) is of the order of some $mm$.  This is not much less
than the current satellite laser--tracking capabilities, thus
providing some motivation for modelling the albedo effect even when
short orbital arcs are analysed (see e.g.  Milani {\it et al.}, 1992).
\begin{figure}
\vspace{76mm}
\caption[]{The same as in Fig.~4, but for the radial
component $S$ of the albedo acceleration.}
\end{figure}

We do not discuss in detail here the $<T>$ values arising from
averaging the resulting values of $T$ over one orbital revolution,
contrary to Rubincam {\it et~al.}, 1987 and Farinella and Lucchesi,
1992.  Although these $<T>$ values may provide some rough estimates
of the long--term effects associated with the different models (see
some remarks on this in Sec. 5), in most cases they sensitively depend
on the assumed orientation of the satellite orbit, the Sun and
Greenwich meridian.  In view of this, we believe that using just {\it
one} particular choice of the above parameters may mislead the
interpretation.  Consider for instance the dominant effect of specular
reflection from the oceans when the satellite is close to the Earth's
shadow boundaries.  As it was pointed out earlier, any asymmetry
caused in these zones by the actual occurrence of continents and/or
clouds (not reflecting specularly) versus oceans may become a strong
contributor to a nonzero $<T>$ value.  Moreover, the albedo
distribution model shows a significant longitude dependence with the
presence of localized structures (see Sehnal, 1979, Figs. 1 and 2),
whose position could also somewhat affect the shape of the $T$ curves.
Thus the real value of $<T>$ may change in a substantial way even in
successive revolutions.  As a consequence, a reliable estimate of the
long--term albedo effects requires a long time span integration (many
thousands of revolutions) of the perturbation equations, followed by a
Fourier transform analysis to point out periodicities.

It is also important to note that the previous discussion applies only
to the intervals when LAGEOS' orbit crosses the Earth's shadow.  Out
of these periods the significant ``kicks'' caused by specularly
reflected radiation are suppressed, since the Fresnel function ${\cal
R}(\mu)$ undergoes a rapid decrease to very small values for
non--tangential illumination of the reflecting surface elements.
Actually, the LAGEOS orbital residuals show sharp peaks and dips
correlated with the Earth's shadow crossing periods.  While other
effects may be responsible for this correlation (Afonso {\it et al.},
1989; Scharroo {\it et al.}, 1991), we plan to carry out a detailed
study of the long--term albedo effects in a forthcoming paper.

\section{Conclusions}
The main results obtained in this paper can be summarized as
follows.

\begin{enumerate}
\item We have provided general integral formulae (1), similar to
the ones appearing in the context of radiative transfer theory, that
express the perturbing acceleration due to radiation pressure for a
satellite of arbitrary shape moving in an arbitrary radiation field.
We have then shown that they can applied in a straightforward way to
the simple case of a nearly spherical, quasi--isotropic satellite such
as LAGEOS.

\item In order to study the albedo orbital perturbations on
LAGEOS, we have discussed in some detail the available models for
the optical behaviour of the Earth's surface, including reflection from
the oceans and possibly anisotropic scattering from clouds.

\item The Gauss components $T$ and $S$ vs. mean anomaly have
been numerically computed for a simple orbital geometry ($90^\circ$
inclination, Sun at an equinox and in the orbital plane) and for
different models of the cloud geographical distribution and optical
properties.  A sensitive dependence has been found for both
components, with significant along--track perturbations possibly
arising from oceanic specular reflection when the satellite is close
to the Earth's shadow boundaries. This may contribute to the observed
eclipse--correlated spikes in the unmodelled $T$ acceleration residuals
found from LAGEOS orbital determination.

\item The radial $S$ component shows differences of the order of
several times $10^{-10}\ m/s^{2}$ when different cloud models are
used, implying short--term errors (when arcs of a fraction of an
orbital period are considered) in LAGEOS' position
pre\-dic\-tions/deter\-mina\-tions of the order of some $mm$.

\item In spite of the {\it caveats} mentioned in Sec.~4, it can
be noted that the difference in the orbit--averaged $<T>$ values
between the no cloud model and the hemisphere model is of the order
$10^{-12}\ m/s^{2}$.  This is consistent with the findings of Lucchesi
and Farinella (1992), who used somewhat different optical models and
another numerical algorithm.  As discussed by these authors, such a
small along--track acceleration may cause on the long term
comparatively large errors in LAGEOS' orbital predictions (or
determinations).  Moreover, this value is of just the same order of
magnitude as those corresponding to other non--gravitational
perturbations, including both drag--like forces and radiative effects
(see the papers referred in Sec.~1).  All these effects are
intrinsically complex and uncertain, and the resulting long--term
perturbations are probably mixed up in the observed acceleration
residuals.  As a consequence, we feel that it is dangerous to try to
fit these residuals by semi--empirical models including a number of
different mechanisms at the same time, as the modeling errors can
cause large biases and/or aliasing effects. Our plan is thus to
further study the albedo effect, investigating the long--term
perturbations associated with the different models introduced in this
paper.
\end{enumerate}

\acknowledgements
We are grateful to F. Mignard for his careful review and many
helpful comments. This work was partially supported by the Italian Space
Agency (ASI) and the EEC research contract no.  SC1--0011--C(GDF).
One of us (D.  V.) wishes to thank CERGA (Grasse, France) and the
University of Pisa (Italy) for kind hospitality.

\begin{thebibliography}{}
\bibitem{}Aksenov, E.P.: 1977, {\it Theory of motion of artificial satellites},
 Nauka (Moscow), in Russian 
\article{Barlier, F., Carpino, M., Farinella, P., Mignard, F., Milani, A.,
 Nobili, A.M.}{1986}{Ann.~Geophys.}{4}{193}
\article{Borderies, N.}{1990}{Celest.~Mech.}{49}{99}
\article{Borderies, N. and Longaretti, P.-Y.}{1990}{Celest.~Mech.}{49}{69}
\book{Chandrasekhar, S.}{1950}{Radiative Transfer}{Oxford Univ. Press}{Oxford}
\book{Ditchburn, R.W.}{1976}{Light}{Academic Press}{London}
\article{Farinella, P., Nobili, A.M., Barlier, F., Mignard, F.}{1990}{A\&A}{234}{546}
\book{van de Hulst, H.C.}{1980}{Multiple Light Scattering : Tables, Formulas,
 and Applications}{Academic Press}{New York}
\article{Irvine, W.\/M.}{1975}{Icarus}{25}{175}
\book{Jackson, J.D.}{1975}{Classical Electrodynamics}{J. Wiley}{New York}
\bibitem{} Knocke, P.Ch.: 1989, {\it Earth Radiation Pressure Effects on
 Satellites}, Dissertation, The University of Texas at Austin
\article{L\'ala, P., Barlier, F. and Oyharcabal, G.}{1978}{Bull. Astron.
 Inst. Czechosl.}{29}{238}
\bibitem{}Levin, E.: 1962, {\it ARS Journal}, 1328 
\article{Lucchesi, D. and Farinella, P.}{1992}{JGR}{97}{7121}
\book{Mihalas, D.}{1978}{Stellar Atmospheres}{Freeman $\&$ co.}{New York}
\book{Milani, A., Nobili, A.M. and Farinella, P.}{1987}{
 Non--Gravitational Perturbations and Satellite Geodesy}{A. Hilger}
 {Bristol}
\bibitem{}Milani, A., Carpino, M., Rossi, A. and Catastini, G.: 1992, 
 {\it JGR}, submitted
\bibitem{}Renard, P. and Koeck, C.: 1989, {\it Skin forces anlysis report},
 ESA Contract Report S374/NT/24/89
\article{Rubincam, D.P.}{1987}{JGR}{92}{1287}
\article{Rubincam, D.P. and Weiss N.R.}{1986}{Celest.~Mech.}{38}{233}
\article{Rubincam, D.P., Knocke, P.Ch., Taylor, V.R., Blackwell, S.}{1987}
 {JGR}{92}{11 662}
\article{Sehnal, L.}{1979}{Bull.~Astron.~Inst.~Czechosl.}{30}{199}
\article{Scharroo, R., Wakker, K.F., Ambrosius, B.A.C. and Noomen, R.}{1991}
 {JGR}{96}{729}
\article{Stephens, G.L., Campbell, G.G. and Vodner Haar, T.H.}{1981}{JGR}
 {86}{9737}
\article{Taylor, V.R. and  Stowe, L.L.}{1984}{JGR}{89}{4987}
\bibitem{}Vokrouhlick\'y, D.: 1991, {\it Notes on the Albedo Effect}, Dissertation,
 Charles University, Prague, in Czech 
\bibitem{} Vokrouhlick\'y, D. and Sehnal L.: 1992a, {\it Celest.~Mech.},
 in press
\bibitem{} Vokrouhlick\'y, D. and Sehnal L.: 1992b, {\it Celest.~Mech.},
 submitted
\bibitem{}
\end{thebibliography}
\end{document}
