% IMPORTANT NOTES:
% - Figures avilable on request from the authors.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Slow-rotation LAGEOS thermal effects, I
% Draft revised by PF, 31 December 1995
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle{elsart}
\begin{document}
\begin{frontmatter}
\title{Thermal force effects on
 slowly rotating artificial satellites. I. Solar heating}
\author[Pisa]{Paolo Farinella} and
\author[Prague]{David Vokrouhlick\'y}
\address[Pisa]{Gruppo di Meccanica Spaziale, Dipartimento di
     Matematica, Universit\`a di Pisa, Via Buonarroti 2,
     56127 Pisa, Italy}
\address[Prague]{Institute of Astronomy, Charles University,
     \v Sv\'edsk\'a 8, 15000 Prague 5,\\ Czech Republic}

\begin{abstract}
The heating of a spinning artificial satellite by natural radiation
sources such as the Sun and the Earth results into temperature
gradients arising across the satellite's surface.  The corresponding
anisotropic emission of thermal radiation leads to a recoil force,
commonly referred to as ``thermal force''.  We have developed a
quantitative theory of this effect, based on more general assumptions
than used so far to model such radiation forces on LAGEOS--like
satellites.  In particular, our theory holds for any ratio of the
three basic time scales of the problem: the rotation period of the
satellite, the orbital period around the Earth, and the relaxation
time for the thermal processes.  Thus, we avoid the simplifying
assumption of a comparatively fast rotational motion, which will fail
for LAGEOS within the next decade, owing to magnetic dissipation
effects.  A number of predictions about the future behaviour of
non--gravitational long--term orbital perturbations of LAGEOS become
possible with the new theory.  Here we have studied in particular the
Yarkovsky--Schach thermal force effects arising as a consequence of
the solar radiation flux onto the satellite, periodically interrupted
by eclipses.  Starting on about year 2005, the orbital perturbation
effects predicted by the new theory are substantially different from
those inferred in the fast--rotation case.  This holds not only for
the long--term semimajor axis effects, but also for eccentricity and
inclination perturbations.
\end{abstract}

\end{frontmatter}

\section{Introduction}
The passive artificial satellite LAGEOS is the celestial body with the
most precisely known orbit.  Current laser tracking technology allows
an accuracy of the order of $1$~cm in determining its position, with
the available data spanning about two decades.  As a result, LAGEOS
orbit analysis has led to derive valuable information on a number of
geophysical processes, such as Earth rotation, geopotential harmonics
and their time variations, tectonic plate motions, and tides.  At the
same time, the tracking data have shown that dynamical models aimed at
reconstructing/predicting the long--term evolution of LAGEOS' orbital
elements (or better, of their time derivatives, or {\it excitations})
should include a number of tiny non--gravitational forces, related to
the interaction of the satellite with the radiation and particle flux
present in space.  The study of these dynamical effects has presented
a number of complex problems, at the borderline between celestial
mechanics and environmental physics.

The most precisely determined orbital element is the semimajor axis
$a$.  Its changes ($\approx 1$~mm/day for LAGEOS, over a value
of about $12,270$~km) can be translated into a time series for the
along--track component of the perturbing force acting on the
satellite.  The data clearly show a secular decrease of $a$,
corresponding to an average drag--like perturbing acceleration of
$-3.4\times 10^{-12}$~m/s$^2$, which must have a non-gravitational
origin.  This secular effect is modulated by long--periodic
oscillations, for which a number of both gravitational and
non--gravitational mechanisms have been suggested (Rubincam, 1982;
Barlier {\it et~al.}, 1986; Milani {\it et~al.}, 1987; Scharroo {\it
et~al.}, 1991).

Rubincam (1987, 1988, 1990) proposed that the dominant role in
generating the secular part of this perturbation is played by
``thermal thrust'' (or ``drag'') effects --- a recoil force due to
thermal radiation anisotropically emitted from the satellite's
surface, which is differentially heated by external radiation sources.
Due to the relatively fast rotation of LAGEOS, the temperature was
considered as constant along any ``parallel'' on the satellite's
surface, and only latitude--dependent temperature changes were seen as
relevant --- giving rise to a ``seasonal'' variant of the classical
Yarkovsky thermal effect known in solar system dynamics since a long
time (for a review, see Burns {\it et~al.}, 1979).  According to
Rubincam's analysis, up to $70\%$ of the secular semimajor axis decay
and a significant portion ($\approx 50\%$) of the long--term
oscillations are caused by this thermal effect, with the infrared
radiation flux coming from the Earth as the source of differential
heating of the satellite surface (hereinafter this will be referred to
as {\it Yarkovsky--Rubincam effect}).  Soon afterwards, it was
realized that the eclipse--related spikes apparent in the semimajor
axis residuals may also be explained by a similar effect, triggered by
the eclipse--modulated visible sunlight flux (hereinafter this will be
called {\it Yarkovsky--Schach effect}; see Slabinski, 1988, Afonso
{\it et~al.}, 1989, Scharroo {\it et~al.}, 1991).  Other likely
contributors to the observed semimajor axis residuals include neutral
and charged particle drag, Earth--albedo radiation pressure, solar
radiation pressure at the Earth's penumbra crossings, and a possible
anisotropy in the optical properties of the satellite's surface.

At the beginning of 90's, it became clear that a further
complicating factor should be taken into account if the dynamical
models are aimed at a quantitative match of the orbital data,
including not only the semimajor axis but also the
eccentricity/perigee excitations (see e.g.  Tapley {\it et~al.},
1993).  Since thermal effects depend sensitively on the orientation of
the spin axis and since this direction is not directly observed, one
can either derive it from a dynamical model, or solve for it while
fitting the along--track orbit residuals.  Rubincam applied the
latter, inverse method in analysing the LAGEOS data.  First, he
assumed a constant orientation of the LAGEOS axis, with the
best--fitting solution corresponding well to the in--orbit injection
data reported by the flight control center at launch time.  Later on,
however, it was realized that the spin axis direction could not have
remained the same since the time of launch.  This was suggested by the
unavoidable presence of disturbing torques, and also by the
deficiencies of the constant--orientation model in reproducing the
observed semimajor axis residuals.  To derive the evolution in time of
the spin axis direction, two approaches were again tried: (i) an
empirical one, determining from the orbital data a set of
best--fitting spin axis directions over a sequence of relatively short
intervals of time, and fitting these successive positions with a
regular curve (Ries {\it et ~al.}, 1993; Robbins {\it et~al.}, 1994);
(ii) a theoretical one, based on a dynamical model for the torques
which affect the rotation of the satellite (Bertotti and Iess, 1991).
In the former approach, like in the original Rubincam work, there is
no way to distinguish between the two senses of rotation about a given
axis in space.  This implicitly led a number of authors not to
distinguish between the two corresponding sets of initial conditions,
both of which were consistent with in--orbit injection data.

However, the situation is different when the evolution of LAGEOS'
rotation is derived from a dynamical model, such as that developed by
Bertotti and Iess.  The differential equations derived by these
authors for the components of the satellite's rotation rate vector are
{\it not} invariant with respect to a change in the sense of rotation,
and therefore it is essential to choose initial conditions
corresponding to the correct sense of rotation.  This has been
recently shown by Farinella {\it et~al.} (1995), who have derived a
satisfactory, self--consistent solution for the coupled evolution of
LAGEOS' spin axis and along--track orbit residuals.  Using the same
model for the spin axis evolution, M\'etris {\it et ~al.} (1995) have
also tackled the problem of adjusting the dynamical parameters to fit
the eccentricity/perigee excitations as well, and have come to the
surprising conclusion that this is possible only provided the value
for the amplitude of the Yarkovsky--Schach effect is about $1.75$
times greater than was expected from the previous thermal models of
the satellite.  This is one of the reasons which have prompted us to
reanalyse the problem of how the thermal effects are modelled, and
thus led to this paper.

Actually, there are several issues which appear worth further study.
First, the original thermal models assumed a homogeneous satellite
body, or anyway an unrealistically simple structure for it.  On the
other hand, it has been recently shown by Slabinski (1995) that the
real complex structure of the satellite results in the fact that its
individual components (the retroreflectors, their ring mountings, the
metallic interior) have different responses and sensitivities to the
visible and infrared heating.  In other words, their individual
thermal histories are significantly different from the average ones
inferred from the ideal homogeneous case.  Slabinski's work is a clear
step forward in modelling these subtle effects, yet the theory still
needs empirical corrections when applied to the observed orbital data.

Second, all the previous analyses missed to include the
Earth--reflected sunlight as a heating source affecting the
temperature distribution across the satellite surface.  M\'etris and
Vokrouhlick\'y (1995) recently determined the main effects of this
additional radiation source.  They concluded that the only significant
perturbation, within the current accuracy of the tracking data, is a
small additional drag--like acceleration, with a constant value of
$\approx -0.3\, {\rm pm/s}^2$.

Third and most relevant here, all the previous thermal force models
adopted the simplifying approximation of a fast rotation of the
satellite, such that only the ``seasonal'' (latitude--dependent)
component of the temperature gradient is significant and therefore the
perturbing force is always directed along the spin axis.  Here ``fast
rotation'' means that the spin period is much shorter than the thermal
relaxation timescale of the satellite structure.  This was a
suitable assumption in the early years of the LAGEOS mission (as the
initial rotation period was less than one second, compared to about
one hour for the thermal relaxation time).  However, according to the
Bertotti and Iess theory and the Farinella {\it et ~al.} solution
mentioned above, LAGEOS' rotation rate is decreasing exponentially
with time under the action of magnetic torques, and in about $10$
years the spin period will reach a value of $3000$~s, comparable to
the thermal lag time.  As a consequence, a more complex thermal and
dynamical model will then be required, with the resulting force
depending also on the sense of rotation.

Therefore, we have decided to develop a more general theory for the
thermal effects acting on the LAGEOS satellite, which removes the fast
rotation assumption.  In other words, we shall not impose any {\it a
priori} constraint on the ratio between the satellite rotation period
and the thermal lag time, but allow for any value of this parameter.
In this situation, ``diurnal'' (i.e., longitudinal) temperature
gradients on the satellite surface will no longer be negligible, and
the net thermal force will tilt from the direction of the spin axis.
The corresponding ``equatorial'' components of the thermal force are
those which appear in the classical formulation of the Yarkovsky
effect (see e.g.  Peterson, 1976; Burns {\it et~al.}, 1979).  Note
that because the thermal lag time is comparable to the revolution
period of LAGEOS around the Earth, we will also have to allow for a
wide range of values of the parameter $r \equiv T_{\rm orb}/T_{\rm
rot}$, namely the ratio of the satellite's revolution period $T_{\rm
orb}$ to its rotation period $T_{\rm rot}$.  The only limitation of
our theory is that it will not be applicable for $r \leq 1$, because
the spin axis evolution of Farinella {\it et~al.} (1995), which we are
going to use, has been obtained by integration of the averaged Euler
equations, which fails at the 1:1 ($r=1$) spin--orbit resonance.

In this paper, we will develop a general theory for the satellite
thermal response to the external heating sources, and then apply it to
the specific case when the heating source is the Sun. In other words,
we are going to generalize the previous treatments of the
Yarkovsky--Schach effect.  A forthcoming paper will be devoted to a
similar general treatment of the Yarkovsky--Rubincam effect, with the
Earth IR radiation as the heating source.

Before proceeding to the mathematical formulation of the problem, it
is worth stressing a few general concepts.  There are three natural
time scales entering our analysis: (i) the satellite rotation period
$T_{\rm rot}$, (ii) the satellite revolution period $T_{\rm orb}$, and
(iii) the thermal relaxation time $\tau_R$.  An even more general
formulation would need to account also for the time scale
characterizing the evolution of the spin axis direction and of the
direction of the radiation source (the Sun in this paper).  In our
particular problem the latter time scales are much longer (one year or
more) than the previously mentioned ones, hence we can consider fixed
orientations of the satellite spin axis and of the heating source over
the longest time among $(T_{\rm rot},T_{\rm orb}, \tau_R)$.

Our goal is to derive expressions for the long--term orbital effects,
that is for the averaged excitations of the mean orbital elements.  It
is thus an obvious requirement of our theory that the revolution
period $T_{\rm orb}$ will disappear from the final solution.  But the
averaging procedure is a somewhat delicate task.  It would be a
mistake to average directly the instantaneous perturbations over both
the rotational and the orbital time scales, because by doing so we
would
obtain the fast--rotation results depending on a single parameter
$\sigma_R$, defined as $(\sigma_R/ 2\pi)\equiv \tau_R/T_{\rm orb}$.
A more careful treatment is needed so as the final results are defined
on a larger domain, parametrized by the two quantities $\sigma_R$ {\it
and} $r$.  Then, by taking the limit $r \rightarrow \infty$ we must
recover the fast--rotation approximation.

The rest of this paper is organized as follows.  In Sec.~2 we give a
general scheme for the treatment of LAGEOS' thermal response to
external heating sources (Sec.~2.1), and compute the corresponding
thermal force due to the solar radiation both with and without the
eclipse effects arising when the orbit crosses the Earth's shadow
(Secs.~2.2 and 2.3, respectively; some technical details of the
calculations will be moved to the Appendix).  In Sec.~3 we apply this
theory to the case of the LAGEOS satellite, and in particular we
discuss how and when the generalized theory will modify the previous
predictions on the long--term evolution of the orbital elements.

\section{Theory}
In the current context, we do not intend to follow Slabinski (1995) in
taking into account the detailed structure of the LAGEOS satellite
when modelling the heat transfer processes which result into the
Yarkovsky--type thermal effects.  We rather use a simple model close
to that of Afonso {\it et~al.} (1989) [see also Rubincam (1987) or
Afonso {\it et~al.} (1995)], namely a spherical satellite with
radius $R$ and homogeneous physical properties.  Like in the
fast--rotation case, this is likely to provide a fairly accurate model
for the time dependence of the thermal effects, but subsequent
empirical readjustments of the numerical parameters appearing in the
final expressions will be required to obtain a good fit of the
observed orbital residuals.  Also, in computing the long--term
dynamical effects, we are going to neglect the eccentricity of the
orbit --- a reasonable approximation, since its value for LAGEOS is $e
\simeq 0.0044$.

\subsection{ Simple thermal model of LAGEOS}
The basic assumptions of our approach are:
\begin{itemize}
\item we develop a linear theory in the
 temperature departure $\Delta T$ from its equilibrium value
 $T_0$ (i.e., terms of the order of $(\Delta T)^2$ are neglected);
\item we neglect the longitudinal and latitudinal heat flow in the
 satellite body, and take into account only the radial heat transfer
 with appropriate boundary conditions.
\end{itemize}

Thus, we assume that a ``weak'' external radiative energy flux $\Delta
\Phi(t;R,{\bf n})$, impinging onto the outer satellite surface element
of normal ${\bf n}$ at the time $t$, causes a temperature perturbation
$\Delta T(t;r,{\bf n})$ in the satellite body, i.e.
\begin{equation}
 T(t;r,{\bf n})  =  T_0 + \Delta T(t;r,{\bf n}) \; .
  \label{eq1}
\end{equation}
The general heat transfer equation is
\begin{equation}
 \kappa \nabla^2\left[\Delta T\left(t;r,{\bf n}\right)\right]
  = \rho C_p {\partial \Delta T\left(t;r,{\bf n}\right)\over
  \partial t} \label{eq2}
\end{equation}
(see e.g. Landau and Lifshitz, 1986). Here
$\kappa$ is the thermal conductivity, $\rho$ the density and
$C_p$ the specific heat. This equation links the divergence of the
heat vector ${\bf q}$ defined as
\begin{equation}
 {\bf q} = -\nabla \left[\Delta T\left(t;r,{\bf n}\right)\right]
  \; , \label{eq3}
\end{equation}
to the time changes of the temperature field $\Delta T$.

In our case the boundary condition, expressing the energy balance
between the external radiation flux impinging onto the outer surface
of the satellite and the re--emitted radiation (together with the
transferred heat), takes the following form:
\begin{equation}
 \left[\kappa {\partial \Delta T\left(t;r,{\bf n}\right)\over
  \partial r}\right]_R + 4\epsilon \sigma T_0^3 \Delta T\left(
  t;R,{\bf n}\right) = \alpha \Delta \Phi\left(t;R,{\bf n}
  \right)\; , \label{eq4}
\end{equation}
where $\epsilon$ is the surface's thermal emissivity, $\alpha$ is its
absorption coefficient and $\sigma$ is the Stefan--Boltzmann constant.
In (\ref{eq4}) we have neglected the latitudinal and longitudinal
components of the heat flux in the satellite body, taking into account
only the radial heat transfer in the satellite surface layer
(instead of projecting the heat vector ${\bf q}$ onto the normal ${\bf
n}$, in the first term we have just taken its radial component $q_R$).
Similarly, the Laplace operator in (\ref{eq2}) can be simplified to
the radial terms only, because the boundary condition (\ref{eq4})
allows one to decouple all the non--radial eigenfunctions of the
temperature diffusion problem.  The same simplifying assumption has
been directly or indirectly adopted in all previous studies on this
issue, and we shall also rely on it.

In any case, the system of equations (\ref{eq2}) and (\ref{eq4})
represents a fairly complicated mathematical problem and different
approximation methods can be envisaged for its solution.  One
possibility, not yet employed in the treatment of satellite thermal
problems, would be to rewrite the Fourier equation (\ref{eq2}) with
the boundary condition (\ref{eq4}) directly into the form of an
equation for the surface temperature $\Delta T(t;R,{\bf n})$, looking
for a proper junction of the solution at the centre of satellite.  The
latter must satisfy a system of Volterra integral equations of the
second kind, as shown by Landau and Lifschitz (1986) in case of a
one--dimensional problem.  This rigorous method, however, is quite
complex and lies beyond the scope of this paper.

In the following, we shall adopt a much simpler approximate method,
introduced by Rubincam (1987).  It provides a fairly simple,
analytical solution of the satellite thermal problem and, very
important in this context, also allows one to carry out the next step
--- the computation of orbital perturbations --- analytically.
Actually, this method leads to an exponential factor expressing the
``thermal memory'' of the body [see Eq.~(\ref{eq8})], which
corresponds to the temporal part of the eigenfunctions constituting
the Green's function of the diffusion problem (Morse and Feshbach,
1953).

A suitable tool for solving the parabolic diffusion equation
(\ref{eq2}) is the use of the Fourier transform (see Rubincam, 1987).
By labelling the Fourier--transformed quantities with a superscript
$F$, the solution for the temperature change $\Delta T$ related to the
variable heating flux $\Delta \Phi$ is approximately
\begin{equation}
 \Delta T^F(\nu) \simeq {3\alpha\over \rho C_p R}\,
  \Delta \Phi^F(\nu)\; {1 \over \nu_R + \imath \nu}\; , \label{eq5}
\end{equation}
(here $\imath$ means the complex unit) with
\begin{equation}
 \nu_R \equiv {1\over \tau_R} \simeq {12\epsilon\sigma T_0^3\over
 \rho C_p R}\; , \label{eq6}
\end{equation}
The exact numerical coefficients in these formulae are not essential
because of our simplifying homogeneity assumption and the poor
knowledge of the LAGEOS ``mean'' material constants.  The fundamental
property of solution (\ref{eq5}) is its functional dependence of
frequency.  It is also important to note that this solution suffers
from a ``high--frequency cut'', which rigorously limits its validity
to the case when the spectrum of $\Delta\Phi(t)$ has a sufficiently
fast decay above the limit
\begin{equation}
 \nu_{\rm cut} \simeq {\kappa\over \rho C_p R^2}\; .\label{eq7}
\end{equation}
In practice this means that timescales shorter than $10-100$~s should
be missing in the spectrum of the radiation flux. This condition is
satisfied in our case.

The computation of the inverse Fourier transform of the solution
(\ref{eq5}) leads to
\begin{equation}
 \Delta T(t;R,{\bf n}) = {3\alpha \over \rho C_p R}
  \int_{t' < t} dt'\,  \Delta \Phi(t';R,{\bf n}') e^{-\nu_R(t-t')}
  \; , \label{eq8}
\end{equation}
with the exponential damping kernel mentioned above.  The effect of
the satellite's rotation is hidden by the fact that in the integrand
of Eq.~(\ref{eq8}) one must consider surface elements ${\bf n}'$ at
the retarded times $t'$ which eventually become identical with the
surface element ${\bf n}$ at time $t$, appearing in the left--hand
side of this equation.  In our case, the uniform rotational motion is
to be expressed by Eq.~(\ref{eq12}) below.

Once the temperature field $\Delta T(t;R,{\bf n})$ on the satellite
surface is known, the recoil force (per unit mass) ${\bf a}$ caused by
the emitted thermal radiation can be derived.  If an isotropic
(``Lambertian'') emission law is assumed, one gets
\begin{equation}
 {\bf a}(t) = -{8\over 3}{\epsilon \sigma\over mc} T_0^3 \int
   dA(\theta,\phi)\, \Delta T[t;R,{\bf n}(\theta,\phi)]\,
  {\bf n}(\theta,\phi)\; , \label{eq9}
\end{equation}
where the surface element is $dA(\theta,\phi)=R^2\sin\theta
d\theta d\phi$. The integral in the right--hand side of
Eq.~(\ref{eq9}) has to be performed over the whole satellite
surface.

An alternative way to compute ${\bf a}$ is by interchanging the two
successive integrations in (\ref{eq8}) and (\ref{eq9}).  We thus
obtain
\begin{equation}
 {\bf a}(t) = -{8\over mc} {\alpha\epsilon\sigma\over \rho C_p R}
 T_0^3 \int_{t'<t} dt' \, e^{-\nu_R(t-t')}\, {\bf h}\left(t,t';R
  \right) \; , \label{eq10}
\end{equation}
where the ``kernel function'' ${\bf h}(t,t';R)$ is defined by
\begin{equation}
 {\bf h}\left(t,t';R \right) = \int_{\Omega_\star(t)} dA(\theta,\phi)
  \, \Delta \Phi\left(t';R,{\bf n}'\right) {\bf n}(\theta,\phi)\; .
  \label{eq11}
\end{equation}
The integration domain $\Omega_\star(t)$ in the previous integral
is implicitly defined by the illumination condition $\Delta
\Phi(t;R,{\bf n}) \geq 0$.  It can be easily expressed explicitly in
the coordinates $(\theta, \phi)$.  Note that the integration in the
right--hand side of (\ref{eq11}) is performed over the domain
$\Omega_\star(t)$ given at the time $t$, whereas the flux projection
$\Delta \Phi(t';R,{\bf n}')$ in the integrand is taken at the retarded
time $t'$.

We stress the different interpretations which can be made for the
approaches based on (\ref{eq8}) -- (\ref{eq9}) and (\ref{eq10}) --
(\ref{eq11}), respectively, to compute the thermal force.  In the
former case, which more closely corresponds to intuition, one first
establishes at a given time $t$ the temperature field on the satellite
surface [integral (\ref{eq8})] and then considers the contribution of
all surface elements to the resulting thermal force [integral
(\ref{eq9})].  In the latter case, one starts by computing the first
multipole moment of the illumination function $\Delta \Phi(t';R,{\bf
n}')$ [with the need to distinguish the time $t$ and the retarded
instant $t'$; integral (\ref{eq10})] and then takes into account the
thermal relaxation processes characterized by an exponential damping
[integral (\ref{eq11})].  We will find the second method much more
suitable for the computation of the thermal force in the general
non--stationary case (Sec.~2.3).

\subsection{ Yarkovsky--Schach effect, with no eclipses}
We now introduce the following local reference frame: (i) the unit
vector ${\bf e}_Z$, characterizing the direction of the $z$--axis,
along the instantaneous direction of the satellite's spin axis (${\bf
s}$), (ii) the unit vector ${\bf e}_X$, characterizing the direction
of the $x$--axis, chosen in such a way that the $xz$--plane contains
the instantaneous direction towards the Sun (${\bf n}_s$).  We then
define a system of spherical coordinates $(\theta,\phi)$, attached
to this reference frame in the usual manner.  We thus have for
instance ${\bf n}_s \equiv (\sin\theta_s, 0,\cos\theta_s)^{\rm T}$.
Figure~\ref{fig1} shows the geometrical quantities introduced above.

\begin{figure}
%\vspace{70mm}
\vspace{33mm}
\caption{ Reference frame and geometrical quantities introduced in the
text. ${\bf s}$ is the satellite's spin vector, ${\bf n}_s$ gives the
direction of the Sun and ${\bf n}$ is the normal to a satellite
surface element.
 \label{fig1}}
\end{figure}

To define uniquely our problem we have still to specify two
properties: (i) the rotation law of the satellite, and (ii) the
external radiation flux.  Consider a particular surface element
$dA(\theta_0, \phi_0)$, parametrized by its position
$(\theta_0,\phi_0)$ at an arbitrary initial time $t_0$.  Then the
rotation of the satellite is expressed by the position of this element
at any time $t$:
\begin{equation}
 \phi(t) = 2\pi\nu_r\, t + \phi_0 \,\, , \theta =\theta_0\; .
\label{eq12} \end{equation}
Here, $\nu_r \equiv 1/T_{\rm rot}$ is the rotational frequency of the
satellite.  As for the surface illumination, we assume a homogeneous
radiation field with intensity $I_0$, specified by the solar direction
${\bf n}_s$.  Then the scalar radiation flux through the surface
element characterized by a normal vector ${\bf n}(\theta,\phi)$ is
given by
\begin{equation}
 \Delta \Phi(R,{\bf n}) = I_0 \left[\sin\theta_s \sin\theta
 \cos\phi +\cos\theta_s \cos\theta\right]\; . \label{eq13}
\end{equation}
By definition $\Delta \Phi$ must be positive for an illuminated
surface element, otherwise its value is zero.  From the physical point
of view, $I_0$ is the ``effective value'' of the solar constant,
including its annual variations due to the Earth's orbital
eccentricity.  Note that on the time scales $(T_{\rm rot},T_{\rm
orb},\tau_R)$ relevant for our problem, the radiation flux
$\Delta\Phi$ does not depend {\it explicitly} on the time $t$, that is
we can neglect the slow, ``adiabatic'' changes of the spin axis
direction and the position of the Sun. This fact provides an important
simplification of the problem in the case when the orbit undergoes no
eclipse.

Another important implication of this assumption is that the orbital
time scale is totally absent in the physical description of the
satellite thermal processes.  It is thus natural to relate all the
time periods to the basic time scale given by the satellite's
rotation, $T_{\rm rot}$.  As a result, for the purposes of this
section we can define the auxiliary quantities $(\sigma_R'/2\pi)
\equiv (\tau_R/T_{\rm rot})$ and ${\bar e}' \equiv
\exp(-2\pi/\sigma_R')$, different from $\sigma_R$ as defined in
Sec.~1.  Changing the integration variable in Eq.~(\ref{eq8}) using
(\ref{eq12}), we obtain for the temperature field $\Delta T$
\begin{equation}
 \Delta T(R,{\bf n}) = {3\alpha\over 2\pi\nu_r\rho
  C_p R} \left(1-{\bar e}'\right)^{-1} \int_{\phi-2\pi}^{\phi}
  d\phi'\,  \Delta \Phi(R,{\bf n}')\, e^{-\sigma_R'
  (\phi-\phi')} \; . \label{eq14}
\end{equation}
This integral can be computed easily for any surface element ${\bf n}$
(we do not give here the corresponding lengthy but straightforward
formulae).  Eq.~(\ref{eq9}) then yields the components of the thermal
force per unit mass as
\begin{eqnarray}
 a_Z & = & \gamma\, \cos\theta_s \; , \label{eq15} \\
 a_X & = & \gamma\, {\sin\theta_s\over 1 + \sigma_R'^2} \; ,
  \label{eq16} \\
 a_Y & = & \gamma \,{\sin\theta_s\sigma_R'\over 1 + \sigma_R'^2}
  \; ,  \label{eq17}
\end{eqnarray}
with the normalizing factor $\gamma$ given by
\begin{equation}
 \gamma = - {16\over 3} {\pi R^2 I_0\over mc} \tau_R T_0^3 {\alpha
 \epsilon\sigma\over \rho C_p R}\; . \label{eq18}
\end{equation}
The spin--directed component $a_Z$ coincides with the usual
fast--rotation expression (Afonso {\it et~al.}, 1989; Farinella {\it
et~al.}, 1990), while the equatorial components vanish in the same
approximation ($\sigma_R' \rightarrow \infty$).  However, as soon as
the LAGEOS rotation will be slow enough, i.e. $T_{\rm rot} \simeq
\tau_R$, the latter components will no more be negligible.  Their
perturbing effects will be discussed in Sec.~3.

The perturbing acceleration ${\bf a}$ from (\ref{eq15}) --
(\ref{eq17}) does not contribute to the averaged along--track force
because it does not depend on the satellite's position along the
orbit, hence it does not affect the semimajor axis evolution (since we
are neglecting the eccentricity terms).  However, it does contribute
to the eccentricity/perigee excitation.  Denoting $\Delta {\dot h}$
the real part and $\Delta {\dot k}$ the imaginary part of the
excitation (i.e., $h=e\cos\omega$ and $k=e\sin\omega$, $\omega$ being
the argument of perigee; see Tapley {\it et~al.}, 1993; M\'etris {\it
et~al.}, 1995) we have
\begin{eqnarray}
 \Delta {\dot h} & = & {3\over 2na} \left({\bf a}\cdot{\bf \hat b}
  \right) + {\cal O}(e) \; , \label{eq19} \\
 \Delta {\dot k} & = & -{3\over 2na} \left({\bf a}\cdot{\bf \hat a}
  \right) + {\cal O}(e) \; , \label{eq20}
\end{eqnarray}
where $n \equiv (2\pi/T_{\rm orb})$ is the mean motion and $a$ the
semimajor axis of the orbit.  The two unit vectors in the previous
equations are defined as follows: ${\bf \hat a}$ is along the
direction of the ascending node, ${\bf \hat b}$ lies in the orbital
plane and $({\bf \hat a},{\bf \hat b})$ together with the orbital
angular momentum vector form a direct triad of vectors.

\subsection{ Yarkovsky--Schach effect, including eclipses}
Contrary to the previous case, the problem with the satellite's orbit
crossing the Earth's shadow is no more stationary, because the
radiation flux $\Delta\Phi$ in (\ref{eq13}) depends {\it explicitly}
on time due to the interruption of the heating flux which occurs
during eclipses.  As a consequence, the thermal force is no longer
constant along one satellite revolution around the Earth, and
semimajor axis perturbations arise.  To better explain our method of
derivation in this case, we start with the special situation of a
$N$:1 spin--orbit resonance, i.e. with the ratio $r$ between the
orbital and rotational periods coinciding with an
integer number $N$.
We stress that this assumption has no specific dynamical meaning here,
but it is just a convenient ``trick'' to simplify the mathematical
derivation.  At the end of this section and in more details in the
Appendix, we will describe how the results can be generalized to any
real value of the parameter $r$.

As mentioned in Sec.~2.1, the approach based on Eqs.~(\ref{eq10}) --
(\ref{eq11}) for calculating the thermal force is more suitable in the
non--stationary, shadow--crossing case.  We thus start with the
computation of the auxiliary vector quantity ${\bf h}(t,t')$.
Moreover, instead of time we shall use as a variable the satellite
longitude along the orbit $\lambda$, measured from the ascending node
and given by $\lambda = nt + \lambda_0$.  With the resonance trick
described above, we have
\begin{equation}
 \phi(\lambda') = \phi(\lambda) + N \left(\lambda'-\lambda\right)
 \; , \label{eq21}
\end{equation}
to be used in Eq.~(\ref{eq11}). Straightforward algebra leads to
\begin{equation}
 {\bf h}\left(t,t';R \right) = {2\pi\over 3} R^2 I_0
 \pmatrix{\sin\theta_s \cos N \left(\lambda-\lambda'\right) \cr
          \sin\theta_s \sin N \left(\lambda-\lambda'\right) \cr
          \cos\theta_s \cr}  \; . \label{eq22}
\end{equation}

Carrying out the remaining integration (\ref{eq10}) is rather lengthy
but not complicated, and leads the following result for the three
components of the thermal force per unit mass (in the same coordinate
system we introduced in Sec.~2.2)
\begin{eqnarray}
a_Z & = & \gamma \,\cos\theta_s\, \psi_Z(\lambda) \; , \label{eq23} \\
a_X & = & \gamma\,{\sin\theta_s\over 1 + N^2\sigma_R^2}\, \psi_X(
\lambda;N)\; , \label{eq24} \\
a_Y & = & \gamma\,{\sin\theta_s\over 1 + N^2\sigma_R^2}\,
\psi_Y(\lambda;N)\; . \label{eq25}
\end{eqnarray}
Note that now the thermal force depends explicitly on the satellite's
orbital longitude through the ``shadow factors'' $\psi(\lambda)$,
which can be expressed as follows
\begin{eqnarray}
\psi_Z(\lambda) &=& \psi_{\rm shad}(\lambda) + \left(1-{\bar e}
\right)^{-1} \left(e^{\Delta_1/\sigma_R} - e^{\Delta_2/\sigma_R}
\right) \; , \label{eq26} \\
\psi_X(\lambda;N) &=& \psi_{\rm shad}(\lambda) + \left(1-{\bar e}
\right)^{-1} \Bigl[e^{\Delta_1/\sigma_R}\left(\cos N\Delta_1+
N\sigma_R \sin N\Delta_1\right) \nonumber \\ & & \qquad\qquad\qquad -
\,e^{\Delta_2/\sigma_R}\left(\cos N \Delta_2+N\sigma_R\sin
N\Delta_2\right)\Bigr] \; , \label{eq27}\\
\psi_Y(\lambda;N) &=& N\sigma_R\psi_{\rm shad}(\lambda) + \left(
1-{\bar e}\right)^{-1} \Bigl[e^{\Delta_2/\sigma_R}\left(\sin N
\Delta_2-N\sigma_R \cos N\Delta_2\right) \nonumber \\ & &
\qquad\qquad\qquad - \,e^{\Delta_1/\sigma_R}\left(\sin N
\Delta_1-N\sigma_R\cos N\Delta_1\right)\Bigr] \; . \label{eq28}
\end{eqnarray}
Here $\psi_{\rm shad}(\lambda)$ denotes the usual
shadow function, equal to $0$ when the satellite is inside the
geometric shadow and $1$ out of it, ${\bar e} \equiv
\exp(-2\pi/\sigma_r)$ as before and the quantities $\Delta_1$ and
$\Delta_2$ are given by
\begin{equation}
\Delta_1 = \cases{\lambda_1 - \lambda\; , & for $\lambda_1 < \lambda$,
\cr \lambda_1 - \lambda-2\pi \; , & for $\lambda_1 > \lambda$, \cr}
\label{eq29}
\end{equation}
and similarly for the index $2$, with $\lambda_1$ and $\lambda_2$
denoting the orbital longitudes corresponding to entry into and exit
from the Earth's shadow, respectively.  We observe that the polar
force component $a_Z$ again coincides with that obtained using the
fast--rotation assumption (see for instance Afonso {\it et~al.}, 1989,
and M\'etris {\it et~al.}, 1995).  The equatorial components $a_X$ and
$a_Y$ average out in the fast--rotation case $(N \rightarrow \infty$),
as expected owing to the absence of longitudinal (``diurnal'')
temperature gradients.

Neglecting the eccentricity terms, the averaged along--track component
of the perturbing acceleration ${\bf a}(\lambda)$ is given by
\begin{equation}
 T = {1\over 2\pi} \int_0^{2\pi} d\lambda \left[-\sin\lambda\,
  {\bf \hat a} +\cos\lambda\, {\bf \hat b}\right]\cdot {\bf a}(
  \lambda) + {\cal O}(e) \label{eq30}
\end{equation}
(Milani {\it et~al.}, 1987), with the pair of unit vectors $({\bf \hat
a},{\bf \hat b})$ defined as above.

We can also compute separately the along--track acceleration terms
related to the three components of the thermal force
(\ref{eq23}) -- (\ref{eq25}). For this purpose,
we define the following auxiliary projection parameters ($i=1,2$)
\begin{eqnarray}
 A_i^X &=& s_a^X \cos\lambda_i + s_b^X \sin \lambda_i \; ,
  \label{eq31} \\
 B_i^X &=& -s_a^X \sin\lambda_i + s_b^X \cos \lambda_i \; ,
  \label{eq32}
\end{eqnarray}
where $s_a^X \equiv {\bf e}_X\cdot {\bf \hat a}$ and $s_b^X\equiv {\bf
e}_X\cdot {\bf \hat b}$ are the projections of the base vector ${\bf
e}_X$ onto the pair of vectors $({\bf \hat a},{\bf \hat b})$.
Analogous parameters can be introduced for the base vectors ${\bf
e}_Y$ and ${\bf e}_Z$.  In the last case, we have ${\bf e}_Z = {\bf
s}$, so that the $A$ and $B$ parameters correspond to the similar ones
introduced by Scharroo {\it et~al.} (1991) and M\'etris {\it et~al.}
(1995).

Using  (\ref{eq23}), (\ref{eq26}) and (\ref{eq30}) we obtain
\begin{equation}
 T_Z = {\gamma\over 2\pi} {\sigma_R\cos\theta_s\over 1 +
  \sigma_R^2} \left[A_1^Z - A_2^Z + \sigma_R\left(B_1^Z-B_2^Z
  \right) \right] \label{eq33}
\end{equation}
for the mean drag--like acceleration related to the spin--oriented
thermal force component.  Note that this result (i) does not depend on
the parameter $N$, and (ii) coincides with the fast--rotation
expression (see Scharroo {\it et~al.}, 1991).

In a similar way, the equatorial force components yield the following
averaged along--track perturbing accelerations:
\begin{eqnarray}
 T_X & = & {\gamma\over 2\pi} {\sin\theta_s\over \left[
  1+\left(N-1\right)^2\sigma_R^2\right] \left[1+\left(N+1\right)^2
  \sigma_R^2\right]} \times \nonumber \\
 & & \qquad \Bigl\{\left[1+\left(N^2+1\right)\sigma_R^2
  \right] \left(A_1^X - A_2^X\right) \nonumber \\
 & & \quad\qquad\qquad\qquad\qquad +\,\sigma_R \left[1-\left(
  N^2-1\right)\sigma_R^2\right]\left(B_1^X-B_2^X\right)\Bigr\}
  \; , \label{eq34} \\
 T_Y & = & {\gamma\over 2\pi} {N\sigma_R \,\sin\theta_s\over
  \left[1+\left(N-1\right)^2\sigma_R^2\right] \left[1+\left(N+1
  \right)^2\sigma_R^2\right]} \times \nonumber \\
 & & \qquad \left\{\left[1+\left(N^2-1\right)
  \sigma_R^2\right] \left(A_1^Y - A_2^Y\right) + 2\sigma_R\left(
  B_1^Y- B_2^Y\right)\right\} \; . \label{eq35}
\end{eqnarray}
Both the components $T_X$ and $T_Y$ vanish in the fast--rotation
approximation $N \rightarrow \infty$, as expected.  Note, however,
that the former component is characterized by a much faster decay with
rotational frequency ($T_X \propto (1/N^2)$) than the latter one ($T_Y
\propto (1/N)$).  Therefore, the semimajor axis perturbation due to
the $T_Y$ component is more important for the LAGEOS case.  We shall
comment in more detail on the orbital effects associated with $T_X$
and $T_Y$ in Sec.~3.

We do not give here the laborious results for the eccentricity
excitations in the eclipse case.  However, we have verified that in
analogy with the fast--rotation case (see M\'etris {\it et~al.},
1995), the contribution of the extra terms due to the eclipses are
negligible with respect to the contribution of the basic terms
(\ref{eq19}) and (\ref{eq20}).

In the rest of this section, we shall comment on the generalization of
the previous results for any value of the ratio $r$ between the
revolution period $T_{\rm orb}$ and the rotation period $T_{\rm rot}$.
The difficulty is only of a technical character, and is related to the
infinite domain of integration in the right--hand side of
Eq.~(\ref{eq10}).  The simplest way to deal with this problem is that
of determining a suitable time scale for the periodicity of the vector
${\bf h}$.  In the case of the exact $N$:1 resonance discussed above,
we may take one revolution period as the natural ``averaging scale'',
and summing up the geometrical series due to presence of the
exponential factor in the integrand we obtain
\begin{equation}
 {\bf a}(\lambda) = - {4\over \pi mc} {\alpha\epsilon\sigma T_0^3\over
  n\rho C_p R} \left(1-{\bar e}\right)^{-1} \int_{\lambda-
  2\pi}^{\lambda}
  d\lambda' \, e^{-\sigma_R(\lambda-\lambda')}\, {\bf h}\left(
  \lambda ,\lambda' ;R\right) \; , \label{eq36}
\end{equation}
(we have again used the satellite's orbital longitude $\lambda$
instead of the time $t$ as the independent variable).  Alternatively,
we can assume a rational representation of the parameter $r$, i.e., $r
\equiv N/M$.  Because the rational numbers are dense in the set of the
real numbers, by continuity our results will hold also for the general
situation of any real value of $r$.  Because the auxiliary vector
${\bf h}(t,t')$ from (\ref{eq11}) takes now the form
\begin{equation}
 {\bf h}\left(t,t';R \right) = {2\pi\over 3} R^2 I_0
 \pmatrix{\sin\theta_s \cos r \left(\lambda-\lambda'\right) \cr
          \sin\theta_s \sin r \left(\lambda-\lambda'\right) \cr
          \cos\theta_s \cr}  \; , \label{eq37}
\end{equation}
we realize that in this case the suitable time scale for its
periodicity is $M\times T_{\rm orb}$, rather than $T_{\rm orb}$
itself.  Thus modify the previous formula (\ref{eq36}) into
\begin{equation}
 {\bf a}(\lambda) = -{4\over \pi mc} {\alpha\epsilon\sigma T_0^3\over
   n\rho C_p R} \left(1-{\bar e}^M\right)^{-1} \int_{\lambda
   -2\pi M}^\lambda
   d\lambda' \, e^{-\sigma_R(\lambda-\lambda')}\, {\bf h}\left(
   \lambda,\lambda';R\right) \; . \label{eq38}
\end{equation}
Then we must take the time span $M\times T_{\rm orb}$ as the averaging
interval for the along--track perturbing acceleration.  Clearly, the
method is applicable if the final averaged result depends on the
parameter $r$ only, and not on the individual values of $N$ and $M$.
Lengthy calculations confirm this property and show that the final
result is identical with that given in
Eqs.~(\ref{eq33}) -- (\ref{eq35}), provided the integer number $N$ is
substituted by the ratio $r$ (see the Appendix for technical details).

Of some interest might be also the long--term excitation of the mean
orbital inclination $I$ (see Farinella {\it et~al.}, 1990; Tapley {\it
et~al.}, 1993).  We shall not repeat the details of the computation
procedure, which is the same as sketched above, and give only
the final results, separately for the inclination excitations due to
the three components of the thermal force.
The spin--oriented force component yields
\begin{equation}
 \left({dI\over dt}\right)_Z = {\gamma\over \pi na} {\cos\theta_s\over
  1+\sigma_R^2}\left({\bf s}.{\bf \hat c}\right) \sin\left({\Delta_{
  \rm sh}\over 2}\right) \left(\cos\lambda_s - \sigma_R \sin
  \lambda_s\right) \; . \label{eq39}
\end{equation}
Here, $\lambda_s \equiv (\lambda_1+\lambda_2)/2-\pi$ is the solar
longitude measured along the satellite orbit from its node and
$\Delta_{\rm sh} \equiv (\lambda_2-\lambda_1)/2$ is the angular width
of the Earth's shadow measured in the orbital plane.  This result is
identical with that given by Farinella {\it et~al.} (1990), an obvious
consequence of the fact that this force component exactly coincides
with the result of the fast--rotation approximation for the thermal
effects.  As for the contributions of the equatorial thermal force
components, they are
\begin{eqnarray}
 \left({dI\over dt}\right)_X & = & {\gamma\over \pi na} {\left({\bf
  e}_X.
  {\bf \hat c}\right)\sin\theta_s\over \left[1+\left(r-1\right)^2
  \sigma_R^2\right] \left[1+\left(r+1\right)^2\sigma_R^2\right]}
  \sin\left({\Delta_{\rm sh}\over 2}\right)\times  \nonumber \\
 & & \;\; \left\{\left[1+\left(r^2+1\right)\sigma_R^2
  \right]\cos\lambda_s-\sigma_R \left[1-\left(r^2-1\right)
  \sigma_R^2\right]\sin\lambda_s\right\} \; , \label{eq40} \\
 \left({dI\over dt}\right)_Y & = & {\gamma\over \pi na} {\left({\bf
  e}_Y.
  {\bf \hat c}\right)\, r\sigma_R\,\sin\theta_s\over \left[1+\left(
  r-1\right)^2\sigma_R^2\right] \left[1+\left(r+1\right)^2\sigma_R^2
  \right]} \sin\left({\Delta_{\rm sh}\over 2}\right)\times
  \nonumber \\
 & & \;\; \left\{\left[1+\left(r^2-1\right)\sigma_R^2
  \right]\cos\lambda_s-2\sigma_R\sin\lambda_s\right\} \; .
  \label{eq41}
\end{eqnarray}

\section{Application to LAGEOS perturbations}
In the framework of LAGEOS orbital studies, thermal effects have been
discussed for the first time by Rubincam (1982).  In that paper,
Rubincam correctly identified the thermal effects due to the Earth's
IR radiation as a potentially important perturbation mechanism,
although he did not give any detailed treatment for them.  He also
reported a suggestion of M.~Schach that solar heating combined with
eclipse passages may result into a net perturbing along--track force,
and gave a qualitative description of this effect.  Interestingly,
this phenomenon had been already discussed {\it and measured} by
Boudon {\it et~al.} (1979), in the framework of the French
accelerometeric mission CACTUS.

In this paper, we limit our discussion to the orbital effects related
to the interruption of the solar heating during the eclipse passages,
i.e. to the Yarkovsky--Schach effect.  From a physical point of view,
the effect is due to the simple fact that the satellite cools down
after it enters the Earth shadow, causing a change of the thermal
force, and heats up again after the satellite exits the shadow.  In
other words, following Rubincam (1982) we may say that the finite
thermal lag corresponds to a ``rotated position'' of the effective
Earth shadow, implying that the along--track component of the thermal
force does not average out over one revolution.  In our formulation,
the cooling and heating of the satellite due to the shadow passages
results in the presence of the ``shadow factors'' $\psi(\lambda)$
given by (\ref{eq26}) -- (\ref{eq28}) in Eqs.~(\ref{eq23}) -- (\ref{eq25})
for the instantaneous thermal force.

\begin{figure}
%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl2.pcx}
%\end{picture}
%\end{center} 
\vspace{33mm}
\caption{Shadow factors $\psi(\lambda)$ appearing in the three
components of the thermal force, plotted over one satellite revolution
around the Earth.  We have assumed a rotational period equal to the
orbital one.  The $\psi_Z(\lambda)$ factor appearing in the polar
component is shown by the solid curve, whereas the dashed and
dash--dotted curves represent the $\psi_X(\lambda)/(1+\sigma_R^2)$ and
the $\sigma_R\psi_Y(\lambda)/ (1+\sigma_R^2)$ factors appearing in the
equatorial components, respectively.  The Earth's shadow is assumed to
be centered at $\lambda = \pi$, and the thermal lag parameter
$\sigma_R$ has been set to $2.4$ radians.
\label{fig2}}
\end{figure}

In Figure~\ref{fig2} we plot these factors for $r=1$, i.e. the exact
spin--orbit resonance, and for a thermal lag parameter $\sigma_R =
2.4$ radians, a value close to that of LAGEOS according to analysis of
M\'etris {\it et~al.} (1995).  We centered the shadow position at
$\lambda = \pi$ with $\lambda_1=2.6$, which corresponds approximately
to the maximum width of the shadow for LAGEOS' orbit.  The solid curve
for $\psi_Z(\lambda)$ can be compared with Figure~5 in Afonso {\it
et~al.} (1989), while the shadow factors related to the two equatorial
force components have been never derived before.  Although
qualitatively similar, the three shadow factors are different in the
details.

\begin{figure}
%\vspace{20mm}
%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl2b.pcx}
%\end{picture}
%\end{center}
%\vspace{45mm}

%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl2c.pcx}
%\end{picture}
%\end{center}
\vspace{33mm}
\caption{The shadow factor $\psi_X(\lambda;N)/(1+N^2\sigma_R^2)$
associated with the $X$-component of the thermal force has been
plotted vs. $\lambda$ for different values of $N$, the ratio between
the orbital and the rotational period: $N=1$ (curve~1), $N=3$
(curve~3), $N=5$ (curve~5), and $N=10$ (curve~10).  The parameters in
the upper part a) of the Figure are the same as in Fig.~2, whereas in
in part b) we have taken $\sigma_R = 0.3$ radians.
\label{fig2bis}}
\end{figure}

Figure~\ref{fig2bis} shows the dependence of the $\psi_X(\lambda;N)$
function on the parameter $N$, i.e. on the number of spin cycles
completed during one satellite revolution around the Earth.  In Part
a) the thermal lag parameter $\sigma_R$ has been set again at $2.4$
radians, the value estimated for LAGEOS by M\'etris {\it et~al.}
(1995).  Besides having a decreasing amplitude, for higher values of
$N$ the curves approach a periodic function, with only a temporary
influence of the eclipse passage.  These properties lead to a
rapid decrease of the mean value for increasing spin rates, in
agreement with the expectation that the equatorial contribution to
mean along--track perturbation should average out during one
revolution when the rotation is fast enough.  The qualitative features
shown by the plots can be easily understood by considering that the
force is approximately oriented towards the region on the satellite's
surface with the lowest temperature, and that during a rotational
cycle this region migrates on the satellite surface; at the same time,
the eclipse entry/exit suddenly change the intensity of the heating
flux which causes a temperature gradient to arise on the surface.
Because the large value of the thermal inertia parameter $\sigma_R$ in
part a) of Fig.~\ref{fig2bis} partially hides this phenomenon, in part
b) we have plotted the same quantity but changing the value of
$\sigma_R$ to $0.3$ radians.  We observe that for a fast rotation
(large $N$), the shadow passage excites a high--frequency and strongly
damped wave pattern related to the satellite rotation.  However,
soon this pattern decays to the stationary value of the shadow factor,
due to the low thermal lag parameter $\sigma_R$.

The Yarkovsky--Schach effect is well known to provide a significant
contribution to the peaks/dips occurring in the sequence of
observed LAGEOS semimajor axis residuals when the orbit crosses the
Earth's shadow (Afonso {\it et~al.}, 1989; Scharroo {\it et~al.},
1991; Farinella {\it et~al.}, 1995).  Because here we are primarily
interested in the future influence of the equatorial thermal force
components, we have plotted in Figure~\ref{fig3} the averaged
along--track Yarkovsky--Schach perturbing acceleration starting from
Jan.~1, 2005.  Before this date the LAGEOS rotation is fast enough
that only the polar force component (which remains unaltered in the
fast--rotation case) plays a significant role.

\begin{figure}
%\vspace{7mm}
%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl3.pcx}
%\end{picture}
%\end{center}
\vspace{33mm}
\caption{Predictions of the averaged along--track acceleration on
LAGEOS due to the Yarkovsky--Schach thermal effect vs. time (in days),
starting on Jan.~1, 2005.  The solid curve (1) has been derived from
the generalized theory derived in this paper, including the effect of
the equatorial thermal force components, while the dashed curve (2)
corresponds to the usual fast--rotation approximation.  Very large
differences arise as the satellite's spin rate slows down due to
magnetic dissipation.  At the end of the time span considered in this
plot the 1:1 spin--orbit resonance is reached, and the Bertotti--Iess
dynamical theory used to predict the evolution of LAGEOS' rotational
state ceases to be applicable.  Following M\'etris {\it et~al.}
(1995), we have adopted the parameter values $\gamma = -175\, {\rm
pm/s}^2$ and $\sigma_R = 2.4$ radians.
\label{fig3}}
\end{figure}

The dashed curve in this Figure shows the predictions obtained from
the current fast--rotation case, which according to our results in
about a decade will start to fail in a very substantial way (peaks
instead of dips and {\it vice versa}).  After that, the generalized
theory accounting for the equatorial thermal force components $T_X$
and $T_Y$ will be an essential tool for LAGEOS orbital studies.  Our
simulation has been stopped after about 7 years, i.e. when the 1:1
spin--orbit resonance will be met according to the Farinella {\it
et~al.} (1995) rotational evolution predictions.  When this will
happen, a more complex non--averaged theory will be needed to model
the (possibly chaotic) subsequent spin evolution.  However, it is
clear that the equatorial thermal force components will continue to
play an important role in LAGEOS' long--term orbital evolution.

Figures~\ref{fig4} show our predictions for the eccentricity
excitation effects due to the Yarkovsky--Schach effect.  M\'etris {\it
et~al.} (1995) have recently shown that these effects account for a
significant portion of the observed orbit residuals.  It is thus
important to understand whether the two additional force components
will, in the future, significantly modify the results obtained with
the fast--rotation approximation.

\begin{figure}
%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl4a.pcx}
%\end{picture}
%\end{center}
%\vspace{45mm}

%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl4b.pcx}
%\end{picture}
%\end{center}
\vspace{33mm}
\caption{Predictions on the excitation of LAGEOS' orbital eccentricity
vector due to the Yarkovsky--Schach thermal effect vs. time (in days),
starting on Jan.~1, 1999.  Part a) of the Figure refers to the real
component ($\dot h$) and part b) to the imaginary component ($\dot
k$), both expressed in milliarcseconds per year.  The solid curves
have been derived from the generalized theory derived in this paper,
including the effect of the equatorial thermal force components, while
the dashed curves correspond to the usual fast--rotation
approximation.  Note that the differences between the two sets of
predictions are small for the real part but substantial for the
imaginary part of the eccentricity excitation.  The same parameter
choice as in the previous semimajor axis predictions has been used
here.
\label{fig4}}
\end{figure}

The results shown in Figure~\ref{fig4} suggest that there will be
only a minor contribution of the equatorial thermal force components
to the real eccentricity excitation component (which is affected in a
dominant way by the spin--oriented force component), but that the
opposite will hold for the imaginary component.

Finally, in Figure~\ref{fig5} we have plotted the predicted
inclination excitation due to the Yarkovsky--Schach effect, which was
the subject of Farinella's {\it et~al.} (1990) study in the
fast--rotation approximation case.  As before, that case (dashed
curve, label 2) is now compared to the more general one with the
equatorial components accounted for.  We have used here two different
parameter choices (corresponding to the solid curves 1 and 1'), to
illustrate a point which may become important.  Our results indicate
that an analysis of the inclination effects in the future may allow us
to decorrelate the two parameters appearing in the thermal force
model, i.e. the normalizing factor $\gamma$ and the phase lag
$\sigma_R$.  As we remarked in Sec.~1, the exact values of these
parameters depend in a fairly sensitive way on the detailed thermal
properties and structure of the satellite.  As pointed out by
M\'etris {\it el~al.} (1995), the same purpose may be achieved by an
analysis of the eccentricity excitation components.  But other
perturbing effects overlap the thermal ones in a different way for
different orbital elements, so the inclination data will possibly
supply important additional information.

\begin{figure}
%\begin{center}
%\begin{picture}(100,100)(145,-682)
%\special{em:graph sl5.pcx}
%\end{picture}
%\end{center}
\vspace{33mm}
\caption{Predictions of the averaged inclination excitation on LAGEOS
due to the Yarkovsky--Schach thermal effect vs. time (in days),
starting on Jan.~1, 2000.  The solid curves (1 and 1') have been
derived from the generalized theory derived in this paper, including
the effect of the equatorial thermal force components, while the
dashed curve (2) corresponds to the usual fast--rotation
approximation.  Again significant differences appear as the satellite
rotation period approaches the 1:1 resonance with the revolution
period (at the end of the plotted time span).  The following parameter
values have been used here: $\gamma = -175\, {\rm pm/s}^2$, $\sigma_R
= 2.4$ radians for curves~1 and 2 (in agreement with M\'etris {\it
et~al.}, 1995), and $\gamma = -82\, {\rm pm/s}^2$, $\sigma_R = 1.45$
radians for curve~1' (following Scharroo {\it et~al.}, 1991).  Both
these combinations match equally well the along--track residuals, but
they show slightly different patterns in the inclination series.
\label{fig5}}
\end{figure}

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

\begin{enumerate}
\item We have developed a fairly general theory for the thermal
effects on the dynamics of LAGEOS--type artificial satellites,
avoiding the usual fast--rotation approximation and applicable for any
ratio between the rotational, orbital and thermal relaxation time
scales.  In this paper we have studied in particular the
Yarkovsky--Schach perturbation effects arising from the solar heating
of the satellite and its interruptions by eclipses.  The second paper
of this series will deal with the Yarkovsky--Rubincam effect caused by
Earth IR radiation heating.

\item In the case of LAGEOS, a growing influence of the components of
the thermal force not aligned with the spin axis is predicted in the
near future, due to the rapid despin of the satellite by magnetic
dissipation effects.  Applying the theory developed here will be
essential after about year 2005 to correctly model the
Yarkovsky--Schach perturbations affecting the long--term semimajor
axis decay of lAGEOS' orbit.

\item Important future discrepancies from the fast--rotation case are
predicted also for the perturbations affecting the imaginary part of
the eccentricity excitation and the inclination excitation of LAGEOS.
These effects in some cases are sensitive to the individual values of
two thermal parameters, which cannot be decorrelated from an analysis
of the semimajor axis residuals only.
\end{enumerate}

\noindent{\it Acknowledgements}
P.F. acknowledges partial support from the Italian Space Agency (ASI)
and from the Italian Ministry for University and Scientific Research
(MURST).  D.V. is grateful to the University of Pisa (Italy) for its
kind hospitality and to the Observatoire de la C\^ote d'Azur (Dept.
CERGA, Grasse, France), where he could perform this work thanks
the ``H.~Poincar\'e'' research fellowship.

\section{Appendix}
As announced in Sec.~2.3, this appendix is devoted to some technical
details of the averaging method used to derive the long--term
along--track force.  The method itself is outlined in Sec.~2.3, and is
based on taking a rational approximation for the $r$ parameter:
$r\equiv N/M$, with $N$ and $M$ integer numbers.  Since the
component of the perturbing acceleration directed along the spin axis
does not depend on $N$ and/or $r$, we will limit our discussion here
to the other two components.

First, we define the auxiliary functions
\begin{eqnarray}
 {\hat C}(r,{\bar e}) &=& {1-{\bar e}\cos(2\pi r) \over
  1-2{\bar e}\cos(2\pi r)+{\bar e}^2} \; , \label{eq42}\\
 {\hat S}(r,{\bar e}) &=& {{\bar e}\sin(2\pi r) \over
  1-2{\bar e}\cos(2\pi r)+{\bar e}^2} \; . \label{eq43}
\end{eqnarray}
Then, we proceed to perform the double integration, over
the averaging time\-scale and the thermal relaxation timescale, to
obtain the long--term perturbing acceleration.  We treat separately
the two equatorial components of the thermal force.  After several
rearrangements of the integrals, they can be reduced to
\begin{eqnarray}
 T_X &=& {\gamma\over 2\pi} {\sin\theta_s\over 1+ r^2
  \sigma_R^2} \int_0^{2\pi} d\lambda \left[-\sin\lambda\, {\bf
  \hat a} + \cos\lambda\, {\bf \hat b}\right].{\bf e}_X \times
  \nonumber \\
 &  & \qquad\qquad\qquad\quad \left[{\hat C}\left(r,{\bar e}
  \right)\psi_C\left(\lambda\right)-{\hat S}\left(r,{\bar e}\right)
  \psi_S\left(\lambda\right)\right] \; , \label{eq44} \\
 T_Y &=& {\gamma\over 2\pi} {\sin\theta_s\over 1+ r^2
  \sigma_R^2} \int_0^{2\pi} d\lambda \left[-\sin\lambda\, {\bf
  \hat a} +\cos\lambda\, {\bf \hat b}\right].{\bf e}_Y \times
  \nonumber \\
 &  & \qquad\qquad\qquad\quad \left[{\hat C}\left(r,{\bar e}\right)
  \psi_S\left(\lambda\right)+{\hat S}\left(r,{\bar e}\right)
  \psi_C\left(\lambda\right)\right] \; . \label{eq45}
\end{eqnarray}
Here the $\gamma$ coefficient is the same as in (\ref{eq18}), and
the shadow factors $\psi_C$ and $\psi_S$ are defined as
\begin{eqnarray}
\psi_C(\lambda) &=& \left\{1-{\bar e}\left[\cos (2\pi r) - r\sigma_R
  \sin (2\pi r)\right]\right\} \psi_{\rm shad}(\lambda) \nonumber \\
 & & \quad + e^{\Delta_1/\sigma_R} \left[ \xi^-_1\left(\lambda;r,
  \sigma_R\right)\cos r\Delta_1 + \xi^+_1\left(\lambda;r,\sigma_R
  \right)\sin r\Delta_1 \right] \nonumber \\
 & & \quad - e^{\Delta_2/\sigma_R} \left[ \xi^-_2\left(\lambda;r,
  \sigma_R\right)\cos r\Delta_2 + \xi^+_2\left(\lambda;r,\sigma_R
  \right)\sin r\Delta_2 \right] \; , \label{eq46} \\
\psi_S(\lambda) &=& \left\{r\sigma_R-{\bar e}\left[\sin (2\pi r) +
  r\sigma_R \cos (2\pi r)\right]\right\} \psi_{\rm shad}(\lambda)
  \nonumber \\
 & & \quad - e^{\Delta_1/\sigma_R} \left[ \xi^-_1\left(\lambda;r,
  \sigma_R\right)\sin r\Delta_1 - \xi^+_1\left(\lambda;r,\sigma_R
  \right)\cos r\Delta_1\right] \nonumber \\
 & & \quad + e^{\Delta_2/\sigma_R} \left[ \xi^-_2\left(\lambda;r,
  \sigma_R\right)\sin r\Delta_2 - \xi^+_2\left(\lambda;r,\sigma_R
  \right)\cos r\Delta_2\right] \; . \label{eq47}
\end{eqnarray}
We also introduced
\begin{equation}
 \xi^-_i(\lambda;r,\sigma_R) = \cases{\cos (2\pi r) - r\sigma_R \sin
  (2\pi r)\; , & for  $\lambda_i > \lambda$, \cr 1 \; , & for
  $\lambda_i < \lambda$, \cr} \label{eq48}
\end{equation}
and
\begin{equation}
 \xi^+_i(\lambda;r,\sigma_R) = \cases{\sin (2\pi r) + r\sigma_R \cos
  (2\pi r)\; , & for $\lambda_i > \lambda$, \cr r\sigma_R \; , & for
  $\lambda_i < \lambda$. \cr} \label{eq49}
\end{equation}
For integer values of the parameter $r$
the previous results can be simplified into those given in
(\ref{eq33}) -- (\ref{eq35}).  A
straightforward integration of (\ref{eq44}) and (\ref{eq45}) yields
\begin{eqnarray}
 T_X & = & {\gamma\over 2\pi} {\sin\theta_s\over \left[
  1+\left(r-1\right)^2\sigma_R^2\right] \left[1+\left(r+1\right)^2
  \sigma_R^2\right]} \times \nonumber \\
 & & \quad \Bigl\{\left[1+\left(r^2+1\right)\sigma_R^2
  \right] \left(A_1^X - A_2^X\right) \nonumber \\
 & & \quad\qquad\qquad\qquad\qquad +\,\sigma_R \left[1-\left(
  r^2-1\right)\sigma_R^2\right]\left(B_1^X-B_2^X\right)\Bigr\}
  \; , \label{eq50} \\
 T_Y & = & {\gamma\over 2\pi} {r\sigma_R \,\sin\theta_s\over
  \left[1+\left(r-1\right)^2\sigma_R^2\right] \left[1+\left(r+1
  \right)^2\sigma_R^2\right]} \times \nonumber \\
 & & \quad \left\{\left[1+\left(r^2-1\right)
  \sigma_R^2\right] \left(A_1^Y - A_2^Y\right) + 2\sigma_R\left(
  B_1^Y- B_2^Y\right)\right\} \; , \label{eq51}
\end{eqnarray}
as given in Sec.~2.3 and used in the subsequent analysis of the LAGEOS
orbital effects.

\begin{thebibliography}{9}

\bibitem{A89}Afonso G., Barlier F., Carpino M., {\it et~al.}, Orbital
effects of LAGEOS seasons and eclipses, {\it Ann. Geophys.} {\bf 7},
501--514, 1989.

\bibitem{A95}Afonso G., Gomes R.S. and Florczak M.A., Asteroid
fragments in Earth--crossing orbits, {\it Planet.  Space Sci.} {\bf
43}, 787--795, 1995.

\bibitem{Bar}Barlier F., Carpino M., Farinella P., {\it et~al.},
Non--gravitational perturbations on the semimajor axis of LAGEOS, {\it
Ann. Geophys.} {\bf 4A}, 193--210, 1986.

\bibitem{BI}Bertotti B. and Iess L., The rotation of LAGEOS,
 {\it J. Geophys. Res.} {\bf 96}, 2431--2440, 1991.

\bibitem{Bou}Boudon Y., Barlier F., Bernard A., {\it et~al.},
Synthesis of flight results of the Cactus accelerometer for
accelerations below $10^{-9}$ g, {\it Acta Astronautica} {\bf 6},
1387--1398, 1979.

\bibitem{BNS}Burns J.A., Lamy P.L. and Soter S., Radiation forces on
small particles in the Solar system, {\it Icarus} {\bf 40}, 1--48,
1979.

\bibitem{F90}Farinella P., Nobili A.M., Barlier F. and Mignard F.,
Effects of thermal thrust on the node and inclination of LAGEOS, {\it
Astron. Astrophys.} {\bf 234}, 546--554, 1990.

\bibitem{FVB}Farinella P., Vokrouhlick\'y D. and Barlier F., The
rotation of LAGEOS and its long-term semimajor axis decay: A
self-consistent solution, {\it J. Geophys.  Res.}, submitted, 1995.

\bibitem{LL}Landau L.D. and Lifschitz E.M., {\it Hydrodynamics},
Nauka, Moscow, 1986 (in Russian).

\bibitem{MV2}M\'etris G. and Vokrouhlick\'y D., Thermal force
perturbations of the LAGEOS orbit: The albedo radiation part, {\it
Ann. Geophys.}, submitted, 1995.

\bibitem{MVR}M\'etris G., Vokrouhlick\'y D., Ries J.C. and R.J. Eanes,
Non--gravitational phenomena and the LAGEOS eccentricity excitations,
{\it J. Geophys.  Res.}, submitted, 1995.

\bibitem{MNF}Milani A., Nobili A.M. and Farinella P., {\it
Non--Gravitational Perturbations and Satellite Geodesy}, A. Hilger,
Bristol, 1987.

\bibitem{MF}Morse P.M. and H. Feshbach, {\it Methods of Theoretical
Physics}, McGraw--Hill, New York, 1953.

\bibitem{P76}Peterson C., A source mechanism for meteorites controlled
by the Yarkovsky effect, {\it Icarus} {\bf 29}, 91--111, 1976.

\bibitem{Rie93}Ries J.C., R.J. Eanes and M.M. Watkins, Spin vector
influence on LAGEOS ephemeris, presented at the {\it Second Meeting of
IAG Special Study Group 2.130}, Baltimore, 1993.

\bibitem{Rob94}Robbins J.W., R.G. Williamson and D.P. Rubincam, LAGEOS
I and II spin axis evolution, {\it EOS Trans.  AGU} {\bf 75}, 109,
1994.

\bibitem{R82}Rubincam D.P., On the secular decrease in the semimajor
axis of LAGEOS, {\it Celest.  Mech.} {\bf 26}, 361--382, 1982.

\bibitem{R87}Rubincam D.P., LAGEOS orbit decay due to infrared
radiation from Earth, {\it J. Geophys.  Res.} {\bf 92}, 1287--1294,
1987.

\bibitem{R88}Rubincam D.P., Yarkovsky thermal drag on LAGEOS, {\it J.
Geophys.  Res.} {\bf 93}, 13805--13810, 1988.

\bibitem{R90}Rubincam D.P., Drag on the LAGEOS satellite, {\it J.
Geophys.  Res.} {\bf 95}, 4881--4886, 1990.

\bibitem{S91}Scharroo R., Wakker K.F. and Ambrosius B.A.C., On the
along-track acceleration of the LAGEOS satellite, {\it J. Geophys.
Res.} {\bf 96}, 729--740, 1991.

\bibitem{S88}Slabinski V.J., LAGEOS acceleration due to intermittent
solar heating during eclipses priods, {\it Bull.  Am. Astron.  Soc.}
{\bf 20}, 902, 1988.

\bibitem{S95}Slabinski V.J., A numerical solution for LAGEOS thermal
thrust: The rapid--spin case, {\it Celest.  Mech.}, submitted, 1995.

\bibitem{T93}Tapley, B.D., Schutz B.E., Eanes R.J., Ries J.C and
Watkins M.M., LAGEOS laser ranging contributions to geodynamics,
geodesy, and orbital dynamics, in: {\it Contributions of Space Geodesy
to Geodynamics: Earth Dynamics}, AGU, pp. 147--173, 1993.

\end{thebibliography}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

