\documentclass [b5paper,twoside,12pt]{article}

\usepackage{url}
\usepackage{amsmath}
\usepackage{cite}
\usepackage{subfigure}
\usepackage{color}


\usepackage{graphicx,amsfonts}
\usepackage{here}

\usepackage{bm}
\usepackage{color}

\usepackage{subfigure}

\usepackage{fancyhdr}
\pagestyle{fancy} \footskip = 35pt

\setlength{\textwidth}{28cc} \setlength{\textheight}{42cc}

\begin{document}

\setcounter{page}{1}

\pagestyle{myheadings} \markboth{First Author, Second Author} {On the spacecraft attitude stabilization...}

\title{\bfseries \Large On the spacecraft attitude stabilization in the orbital frame  }

\author{\bfseries
First Author\thanks{Russia, St. Petersburg, 198504, Stary Peterhof, Universitetsky pr., 28, St.-Petersburg State University, Faculty of Mathematics and Mechanics, e-mail:  antipov\_k@rambler.ru}
\qquad Second Author\thanks{Russia, St. Petersburg, 198504, Stary Peterhof, Universitetsky pr., 28, St.-Petersburg State University, Faculty of Mathematics and Mechanics, e-mail: aatikhonov@rambler.ru} }

\date{}
\maketitle


\vspace*{-7cm}

\begin{center}
\textbf{Theoret. Appl. Mech.}, Belgrade 2012,  A manuscript \\[0pt]

\medskip

\rule[5mm]{12.5cm}{0.2mm}\\[-4.5mm]
% \rule[15mm]{12cm}{0.15mm}\\[-10mm]
\end{center}

\vspace*{5cm}

\begin{abstract}

\noindent
 The paper deals with spacecraft in the circular near-Earth orbit.
 The spacecraft interacts with geomagnetic field by the moments of Lorentz
 and magnetic forces. The octupole approxi\-mation of the Earth's magnetic field
 is accepted. The spacecraft electromagnetic parameters, namely the electrostatic
 charge moment of the first order and the eigen magnetic moment are the controlled
 quasiperiodic functions. The control algorithms for the spacecraft electro\-magnetic
 parameters, which allows to stabilize the spacecraft attitude position in the orbital
 frame are obtained. The stability of the spacecraft stabilized orientation is proved
 both analytically and by PC computations. \\[2ex]
 \noindent \textbf{Keywords}: spacecraft, attitude stabilization, geomagnetic field.  \\[1ex]
 \noindent \textbf{Mathematics Subject Classification}: At the file: \\
\noindent http://www.ams.org/mathscinet/msc/pdfs/classifications2010.pdf an author finds for a paper dealing with phase transformations following keys: 74N15 (analysis of microstructure), 74A35 (polar materials) etc. A similar is applied to other fields.

\end{abstract}




\newpage

\section{Introduction}

The forces of spacecraft electrodynamic interaction with the geomagnetic field considerably influence on the spacecraft attitude dynamics and so can be used for designing control systems of spacecraft attitude orientation.


\section{The moment of Lorentz forces}

A spacecraft, whose center of mass moves in the Newtonian central Earth's gravitational field in the Keplerian circular orbit of the radius $R$, is considered. It is assumed that the spacecraft has the electrostatic charge $Q$ distributed within some volume $V$ with the density $\sigma$: $Q=\int\limits_V\!{\sigma}\,dV$.




\subsection{About equations and tables}
According to AMSTeX equations may be split in diverse ways. One is following:
\begin{equation}{{
% \nonumber to remove numbering (before each equation)
\begin{split}
& \mathbf{U}(X_1,X_2,t=0) = \mathbf{0}, \\[1ex]
& \frac{\partial \mathbf{U}}{\partial t} (X_1= \pm 1/2, X_2, t=0) = \pm V_1^{*} \mathbf{e}_{1},  \\[1ex]
& \frac{\partial \mathbf{U}}{\partial t} (X_1, X_2= \pm 1/2, t=0) = \pm V_2^{*} \mathbf{e}_{2},  \\[1ex]
& \frac{\partial \mathbf{U}}{\partial t} ( -1/2 < X_1 < 1/2, -1/2 < X_2 < 1/2 , t=0) = \mathbf{0},
\end{split}
}}
\end{equation}

Tables should be written in the following way.
\begin{table}[!h]
  \centering
  \caption{Experimental frequencies for a woven composite plate.}\label{Table1}
  \vspace*{1mm}
  \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
   \hline
          & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\[1.25ex] \hline
Freq.[Hz] & 37.3 & 85.8 & 110.9 & 211.6 & 222.3 & 235.6 & 257.1 & 268.2 \\[1.25ex] \hline
\end{tabular}
\end{table}

\begin{table}[!h]
  \centering
  \caption{Estimates of the elastic constants for the 3.175 mm thick
aluminum plate.}\label{Table2}
  \vspace*{2mm}
\begin{tabular}{|c|c|c|c|}
\hline
Initial E & Initial & Predicted E  & Predicted G  \\
\hline
[GPa]& [GPa]& [GPa] & [GPa]  \\
\hline
120 & 50 & 71.04& 25.32\\
\hline
115 & 45 & 71.04 & 25.33
\\
\hline
110 & 50 & 71.06 & 25.34
\\
\hline
105 & 55 & 	71.02 & 	25.31
\\
\hline
100	 & 40	 & 71.02	 & 25.32
\\
\hline
95	 & 55	 & 71.04	 & 25.33
\\
\hline
90	 & 50	 & 71.06 & 	25.35
\\
\hline
85	 & 35	 & 71.03 & 	25.32
\\
\hline
80	 & 30	 & 71.06	 & 25.34
\\
\hline
75	 & 35	 & 71.04	 & 25.33
\\
\hline
70	 & 25	 & 71.03	 & 25.32
\\
\hline
65	 & 30	 & 71.05	 & 25.33
\\
\hline
65 & 35	 & 71.03	 & 25.32
\\
\hline
60	 & 30	 & 71.01 & 	25.31
\\
\hline
55	 & 35	 & 71.07 & 	25.32
\\
\hline
\end{tabular}
\end{table}

\begin{table}[!h]
  \centering
  \caption{Estimates of the elastic constants for the 1.6 mm thick
(CFRP) composite plate.}\label{Table3}
\vspace*{2mm}
\begin{tabular}{|c|c|c|}
\hline
 & Initial values  & Predicted values  \\
\hline
& [GPa]& [GPa]   \\
\hline
$E_{11}$ & 	85 & 136.7
\\
\hline
$G_{13}$ & 	5.32 & 5.08
\\
\hline
$G_{23}$ & 	2.59 & 3.16
\\
\hline		
\end{tabular}
\end{table}

\clearpage

\subsection{Figures in eps format}
\begin{figure}[!h]
\centerline{\includegraphics[width=10cm]{fig20.eps}} %{Fig15-Shivaiah.eps}}
\label{fig1}
\caption{A good looking figure acceptable for TAM}
\end{figure}
\vspace*{-10mm}
\begin{figure}[!h]
\centerline{\includegraphics[width=8.5cm]{Fig5.eps}}
\caption{A figure not acceptable for TAM}
\label{fig1}
\end{figure}

\begin{figure}
\centering \subfigure[Snapshot of the space behaviour of $U_1$ at time $\tau = 0.22$.]{\label{fig:U1-symm-022}
%% label for first subfigure
\includegraphics[width = 7.5cm]{SYMM_A2_U1_022.eps}}
\subfigure[Snapshot of the space behaviour of $U_2$ at time $\tau = 0.22$.]{\label{fig:U2-symm-022}
%% label for first subfigure
\includegraphics[width = 7.5cm]{SYMM_A2_U2_022.eps}}
\subfigure[Snapshot of the space behaviour of $U_1$ at time $\tau = 1.48$.]{
\label{fig:U1-symm-148} %% label for second subfigure
 \includegraphics[width = 7.5cm]{SYMM_A2_U1_148.eps}}
 \caption{An acceptable form of figure with subfigures.}
\end{figure}

\clearpage

\section{The third section}
Material symmetry group $\aleph$ of an elastic anisotropic
material with Hooke's tensor $\mathcal{D}$ is defined by all
orthogonal 2-tensors satisfying the relationship: $
    \mathcal{D} = \mathbf{H}\, \diamondsuit \, \mathcal{D}, \quad
    ( \mathbf{H} \in \aleph ),$
where the Rayleigh product explicitly reads: ...

\hspace*{0.5cm}
\parbox[t]{12cm}{
\noindent \textit{Overall symmetry definition}

Given elastic and thermal symmetries of the matrix and $N$
ellipsoidal inclusions (whose semiaxes are defined by the rotation
tensors $\mathbf{R}_{\Lambda} $) as follows: $\;
\mathcal{D}_{\Lambda} = \mathbf{H}_{\Lambda}^e \, \diamondsuit \,
\mathcal{D}_{\Lambda},  \;
    \mathbf{H}_{\Lambda}^e \in \aleph_{\Lambda}^e,  $
as well as $\;  \bm{\alpha}_{\Lambda} = \mathbf{H}_{\Lambda}^t \,
\diamondsuit \, \bm{\alpha}_{\Lambda}, \;
    \mathbf{H}_{\Lambda}^t \in \aleph_{\Lambda}^t,  \; (\Lambda \in
    \{ 0,1, \ldots N \})$
find 2-tensors $\mathbf{H}_{eff}^e \in \aleph_{eff}^e$ and
$\mathbf{H}_{eff}^e \in \aleph_{eff}^e$  such that
\begin{equation}\label{EffElastSym}
    \mathcal{D}_{eff} = \mathbf{H}_{eff}^e \, \diamondsuit \,
\mathcal{D}_{eff},  \quad
    \mathbf{H}_{eff}^e \in \aleph_{eff}^e
\end{equation}
%$\qquad  $ \\[-0.5ex]
as well as...


Groups $\aleph_{eff}^e$ and $\aleph_{eff}^t$ are then called
\emph{effective elastic symmetry group} and \emph{effective
thermal symmetry group} respectively. }

\noindent Obviously, the real task is to find $\aleph_{eff}^e$ and
$\aleph_{eff}^t$ when $\aleph_{\Lambda}^e$ and
$\aleph_{\Lambda}^t, \;(\Lambda \in
    \{ 0,1, \ldots N \}),$ are given.


In the special case when ODF depends only on one angle, say
$\theta,$ averaging over the other two Euler angles leads to the
representation
 \begin{equation} \label{P-2k}
\omega(\mathbf{R}) = \sum_{k=0}^\infty c_{2k} \, P_{2k}(cos
\theta),
\end{equation}
where $P_k, \, k \in\{0,2,4,\ldots\}$ with values $P_0(x) = 1$,
$P_2(x) = (3x^2-1)/2$, $P_4(x) = (35x^4-30x^2+3)/8$, $P_6(x) =
(231x^6-315x^4+105x^2-5)/16, \ldots$ are Legendre functions of
even order.


\medskip

\noindent \textit{Example 1. }

First, suppose that an elastomeric matrix is weakened by some
identical parallel spheroidal voids with symmetry axis aligned
with a Cartesian axis $z_3=x_3$...


\medskip

\section*{Concluding remarks}
Results of this paper may be shortly summarized as follows:
\begin{itemize}
    \item conclusion 1
    \item conclusion 2
    \item conclusion 3
\end{itemize}

\subsubsection*{Acknowledgements}
\noindent This research has been done within the project
MM-144038 financially supported by the Ministry of Science and
Technology of Serbia.


\begin{thebibliography}{99}

\small
 \bibitem{Vujanovic-Jones} {B.D.Vujanovic and S.E.Jones,}
         Variational methods in  nonconservative phenomena, Academic Press, Boston, 1988.

\bibitem{E.Reissner}
        {E. Reissner}, \textit{On one-dimensional finite strain beam theory: the plane   problem},  J.Appl.Math.Phys.(ZAMP), \textbf{23}, (1972), 795-804.

\bibitem{Smith-Johnson}
        {J.Smith and H.Johnson,}  \textit{Waves and turbulence}, In Proc. Euromech. Symp. On Fluid Mechanics, ed.H.Schmidt, Gauthier-Villars, Paris (1993), 34-37.
\end{thebibliography}

\clearpage

\listoffigures
\listoftables
\end{document}

DOVDE

Also, a system $Cxyz$ (basis vectors ${\vec i},{\vec j},{\vec k}$) of spacecraft principal central axes of inertia rigidly bound with the spacecraft is used. The orientation of the orbital coordinate system with respect to the system $O_{\text E}X_*Y_*Z_*$ is defined on the basis of equalities
\begin{equation*}
\begin{array}{c}
\vec{i_*}=-\sin u\, \vec{\xi_0}+\cos u \,\vec{\zeta_0},\ \vec{j_*}=\cos i \cos u\,
\vec{\xi_0}-\sin i \,\vec{\eta_0}+\cos i \sin u\, \vec{\zeta_0},\ \\  \vec{k_*}=\sin i
\cos u\, \vec{\xi_0}+\cos i \,\vec{\eta_0}+\sin i \sin u\, \vec{\zeta_0},
\end{array}
\end{equation*}
where $i=\widehat {({\vec{k_*}},{\vec {\eta_0}})}$ is an orbit inclination; $u=\widehat
{({\vec{i_*}},{\vec {\zeta_0}})}$ is an argument of a latitude and $u=\omega_0 t$ where
$\omega_0$ is an orbital angular velocity of the spacecraft mass center.

The orientation of the axes $xyz$ with respect to the axes $\xi\eta\zeta$ is defined by a matrix $\mathbf A$ of direction cosines $\alpha_i$, $\beta_i$, $\gamma_i \ (i=1,2,3)$ so that there exist the equalities
$${\vec\xi}_0= \alpha_1 {\vec i} + \alpha_2 {\vec j} + \alpha_3
{\vec k}, \quad {\vec\eta}_0= \beta_1 {\vec i} + \beta_2 {\vec j} + \beta_3 {\vec k},
\quad {\vec\zeta}_0= \gamma_1 {\vec i}+ \gamma_2 {\vec j}+ \gamma_3 {\vec k}.
$$
 If we determine the spacecraft orientation in the orbital coordinate system by "airborne" angles $\varphi ,\theta ,\psi$ (Fig.~\ref{Fig2}), then the elements of the matrix $\mathbf A$ will have the form
$$\aligned
\alpha_1 &= \cos\psi\cos\theta, \quad
\alpha_2= -\cos\varphi\sin\psi+ \sin\varphi\cos\psi\sin\theta,  \\
\alpha_3 &= \sin\varphi\sin\psi+ \cos\varphi\cos\psi\sin\theta, \quad
\beta_1=\sin\psi\cos\theta,                     \\
\beta_2 &=\cos\varphi\cos\psi +\sin\varphi\sin\psi\sin\theta,   \quad
\beta_3 =-\sin\varphi\cos\psi+ \cos\varphi\sin\psi\sin\theta,   \\
\gamma_1 &= -\sin\theta, \quad \gamma_2= \sin\varphi\cos\theta, \quad \gamma_3
=\cos\varphi\cos\theta.
\endaligned
$$
%\begin{figure}[h]
%\centering \includegraphics[width=5.5cm,%
%keepaspectratio]{f2} \vspace{-0.5cm}\caption{} \label{Fig2}
%\end{figure}

When the spacecraft is moving with respect to the Earth's magnetic field, the interaction of the shield charge with the geomagnetic field results in Lorentz forces excitation. The principal moment of these forces with respect to the spacecraft mass center is defined by the formula
$$\vec M_{L} =
\int\limits_V\!\sigma\vec\rho\times(\vec v\times\vec B)\,dV,
$$
where ${\vec \rho}$ is the radius-vector of the element $dV$ of the shield with respect to the spacecraft mass center, and $\vec v$ is the velocity of the element $dV$ with respect to the geomagnetic field.
It was demonstrated in \cite {Vestnik1999} that for any real spacecraft and especially for spacecrafts operating in modes close to the oriented attitude motion, the moment $\vec M_{L}$ can be approximated by the expression
\begin{equation}\label{eq2}
\vec M_{L}= \vec{P}\times
  \mathbf A^{\top}(\vec v_C\times\vec B),
\end{equation}
where $\vec{P}=Q\vec\rho_0, \quad {\vec \rho_0}=x_0{\vec i}+ y_0{\vec j}+ z_0{\vec k} =
Q^{-1}\int\limits_V\!\sigma{\vec\rho}\,dV$ is the radius-vector of the center of charge of the spacecraft with respect to its mass center,
\begin{equation}\label{eq33}
\vec v_C=\dot{\vec{R}}-{\vec \omega_{E}}\times{\vec R}=
  R(\omega_0-\omega_{E}\cos i) \, {\vec {\xi_0}}+
  R\omega_{E}\sin i\cos u \, {\vec {\eta_0}}
 \end{equation}
is the velocity of the spacecraft mass center with respect to the geomagnetic filed,
 ${\vec
\omega_{E}}= \omega_{E} {\vec {k_*}}$ is the angular velocity of the Earth's daily rotation. The Earth's magnetic field can be considered in this case to be uniform in the spacecraft volume. Therefore, the value of $\vec B$ in (\ref{eq2}) coincides with the value of $\vec B$ at the spacecraft mass center. Considering the vector
 $\vec T=\mathbf A^{\top}(\vec v_C\times\vec B)$, let us represent the moment $\vec M_{L}$ in the form
\begin{equation*}
\vec M_{L}= \vec P\times\vec T\ .
\end{equation*}

Such approximation of the moment ${\vec M}_{L}$ is admissible in the presence of not very small quantities of
$|{\vec \rho_0}|$ due to the negligible influence of the charge-distribution character within the spacecraft volume on its dynamics. In conditions of the octupole approximation of the geomagnetic field, the vector $\vec B$ is defined by the formula obtained in \cite{KRes2002}:
\begin{equation}
\vec{B}=\sum_{n=1}^{3}\vec{B}^{(n)}=-  {\textrm{grad}}\sum_{n=1}^{3}
\frac{R_{E}^{n+2}}{r^{2n+1}}\text{{\Large M}}^{(n)}\underset{n}{\underbrace{\cdots }}%
\,(\otimes ^{n}\vec{r}),  \label{eq32}
\end{equation}
where $\vec r$ is the radius-vector of a point in the near-Earth space with respect to the Earth mass center,
$\text{{\Large M}}^{(n)}$ are multipole tensors of the first, second and third ranks that are dipole, quadrupole and octupole magnetic moments respectively, earlier obtained in \cite{KRes2002}.

In projections on the orbital coordinate system we have:
\begin{equation}
\begin{pmatrix}
B_{\xi }^{(n)} \\
B_{\eta }^{(n)} \\
B_{\zeta }^{(n)}%
\end{pmatrix}%
=-\left( \frac{R_{E}}{R}\right) ^{n+2}%
\begin{pmatrix}
n\,\text{{\Large M}}_{0}^{(n)}\underset{n}{\underbrace{\cdots }}\text{T}%
_{\lambda }^{(n)} \\
n\,\text{{\Large M}}_{0}^{(n)}\underset{n}{\underbrace{\cdots }}\text{T}%
_{\theta }^{(n)} \\
-(n+1)\,\text{{\Large M}}_{0}^{(n)}\underset{n}{\underbrace{\cdots }}\text{T}%
_{r}^{(n)}%
\end{pmatrix}%
,  \label{eq41}
\end{equation}
where
\begin{align}
\hbox{T}_r & = (\cos u,\, \sin u,\, 0)^\top,
& \hbox{T}_r^{(n)} &= \otimes^{n}\hbox{T}_{r}, \label{eq44} \\
 \hbox{T}_\lambda & =(-\sin u,\, \cos u,\, 0)^\top,
 & \hbox{T}_\lambda^{(n)} & = \hbox{T}_\lambda\otimes^{n-1}\hbox{T}_{r},\label{eq45}\\
  \hbox{T}_\theta & =(0,\, 0,\, 1)^\top,
 & \hbox{T}_\theta^{(n)} & = \hbox{T}_\theta\otimes^{n-1}\hbox{T}_{r},\label{eq46}
\end{align}

$\text{{\Large M}}^{(n)}_0$ - the multipole tensor of $n$-th rank, reformed with matrix $\mathbf{\Gamma}=(\gamma_{ij})$  $(i=\overline{1,3}, j=\overline{1,3})$ according to the formula
$$
\left(\text{{\Large M}}^{(n)}_0\right)_{ i_1,i_2,\ldots,i_n}=\sum^3_{j_1,j_2,\ldots,j_n=1}\gamma_{ i_1,j_1}
\gamma_{i_2,j_2}\dots\gamma_{ i_n,j_n}
{\text{\Large M}}^{(n)}_{ j_1,j_2,\ldots,j_n}.
$$
The components of matrix $(\gamma_{ij})$ are expressed in terms of the orbit inclination $i$ and the hour angle of the ascending node $\phi=\omega_{E}t$ in the following way:
\begin{align}
\gamma_{11} & =\cos \phi, & \gamma_{12}  &= -\sin\phi,& \gamma_{13} &
= 0,\nonumber
\\
\gamma_{21}  &=\cos i \sin\phi, & \gamma_{22} & =\cos i \cos\phi,& \gamma_{23} & = \sin i, \nonumber
\\
\gamma_{31} & =-\sin i \sin \phi, & \gamma_{32}  &= -\sin i \cos \phi,& \gamma_{33} & = \cos i.\nonumber
\end{align}

For completeness of contents we shall cite also the specified expressions for multipole tensors \cite{KRes2002}:
%\vspace{-1cm}
%\begin{align*}
%M_1 &= g_1^1,&
%M_2 &= h_1^1,&
%M_3 &= g_1^0, \\
%\end{align*}
\[ M_1 = g_1^1, \qquad \quad
M_2    = h_1^1, \qquad \quad
M_3    = g_1^0,
 \]
\begin{center}
\rule[5mm]{8cm}{0.2mm}\\[-4.5mm]
\end{center}
%\vspace{-1cm}
%\begin{center}\rule[1cm]{8cm}{0.5pt}\end{center}
\vspace*{-2cm}

\begin{align*}
M_{{11}}&=\frac{1}{2}\left(\sqrt{3}g_2^2-g_2^0\right),& M_{{12}}&=\frac{\sqrt{3}}{2}\,h_2^2,&
M_{{13}}&=\frac{\sqrt{3}}{2}\,g_2^1,%\\
\\[1ex]
M_{{22}}&=-\frac{1}{2}\left(g_2^0+\sqrt{3}\,g_2^2\right),&
M_{{23}}&=\frac{\sqrt{3}}{2}\,h_2^1,&
M_{{33}}&=g_2^0,& %\\
\\[1ex]
M_{{111}}&=\frac{\sqrt {6}}{12}\left(\sqrt{15}\,g_3^3-3\,g_3^1\right),&
M_{{112}}&=\frac{\sqrt {6}}{12}\left(\sqrt{15}\,h_3^3-h_3^1\right),%\\
\\[1ex]
M_{{113}}&=\frac{1}{6}\left(\sqrt{15}\,g_3^2-3\,g_3^0\right),&
M_{{122}}&=-\frac{\sqrt {6}}{12}\left(\sqrt{15}\,g_3^3+g_3^1\right),%\\
\\[1ex]
M_{{123}}&=\frac{\sqrt {15}}{6}\,h_3^2,&
M_{{133}}&=\frac{\sqrt {6}}{3}\,g_3^1%\\
\\[1ex]
M_{{222}}&=-\frac{\sqrt {6}}{12}\left(\sqrt{15}\,h_3^3+3\,h_3^1\right),&
M_{{223}}&=-\frac{1}{6}\left(\sqrt{15}\,g_3^2+3\,g_3^0\right),%\\
\\[1ex]
M_{{233}}&=\frac{\sqrt {6}}{3}\,h_3^1,&
M_{{333}}&=g_3^0,
\end{align*}

In detailed form the components of vector $\vec B$ in octupole approximation are the following:
\begin{equation*}
\begin{split}
B_{\xi } =&-\left( \dfrac{R_{E}}{R}\right) ^{3}\left( \left( f_{{0,1,0,1,1,0%
}}-f_{{1,0,0,0,0,1}}\right) M_{1}+\left( f_{{1,0,0,0,1,0}}+f_{{0,1,0,1,0,1}%
}\right) M_{2}\right.
\\
&\left. +M_{3}f_{{0,1,1,0,0,0}}\right) -2\left( \dfrac{R_{E}}{R}\right)
^{4}\left( 2\,M_{13}f_{{0,2,1,0,0,1}}+\left( M_{22}-M_{11}\right) f_{{%
1,1,0,0,0,2}}\right.
\\[1ex]
&\left. +\left( M_{22}-M_{11}\right) f_{{0,0,0,1,1,1}}+\left(
-M_{22}+M_{33}\right) f_{{1,1,0,0,0,0}}+M_{23}f_{{0,0,1,0,1,0}}\right.
\\[1ex]
&\left. +M_{12}f_{{0,0,0,1,0,0}}-M_{13}f_{{0,0,1,0,0,1}}+\left(
2\,M_{11}-2\,M_{22}\right) f_{{0,2,0,1,1,1}}\right.
\\[1ex]
&\left. +2\,M_{{12}}f_{{1,1,0,0,1,1}}+2\,M_{12}f_{{1,1,0,2,1,1}%
}+2\,M_{13}f_{{1,1,1,1,1,0}}\right.
\\[1ex]
\end{split}
\end{equation*}%
\begin{equation}  %\begin{equation*}
\begin{split}
&\left. +\left( M_{22}-M_{11}\right) f_{{1,1,0,2,0,2}}+2\,M_{{23}}f_{{%
1,1,1,1,0,1}}-2\,M_{12}f_{{0,2,0,1,0,0}}\right.
\\[1ex]
&\left. -2\,M_{12}f_{{0,0,0,1,0,2}}+\left( M_{11}-M_{33}\right) f_{{%
1,1,0,2,0,0}}+4\,M_{12}f_{{0,2,0,1,0,2}}\right.
\\[1ex]
&\left. -2\,M_{23}f_{{0,2,1,0,1,0}}\right) -3\left( \dfrac{R_{E}}{R}\right)
^{5}\left( \left( 3\,M_{112}-M_{222}\right) f_{{1,2,0,0,1,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{133}-3\,M_{122}\right) f_{{1,2,0,0,0,1}}-6\,M_{123}f_{%
{0,3,1,2,1,1}}+6\,M_{123}f_{{0,1,1,2,1,1}}\right.
\\[1ex]
&\left. +\left( 2\,M_{{223}}-2\,M_{113}\right) f_{{1,0,1,1,1,1}}+\left(
-3\,M_{122}+M_{{111}}\right) f_{{0,3,0,3,1,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{233}-3\,M_{{112}}\right) f_{{1,2,0,2,1,0}}+\left(
M_{222}-3\,M_{112}\right) f_{{1,0,0,2,1,2}}\right.
\\[1ex]
&\left. -4\,M_{123}f_{{1,0,1,1,0,2}}-6\,M_{{123}}f_{{1,2,1,1,0,0}}\right.
\\[1ex]
&\left. +\left( -6\,M_{122}-3\,M_{133}+3\,M_{{111}}\right) f_{{1,2,0,2,0,1}%
}\right.
\\[1ex]
&\left. +\left( -2\,M_{111}+6\,M_{122}\right) f_{{0,1,0,1,1,2}}+\left(
-3\,M_{113}+3\,M_{223}\right) f_{{0,1,1,2,0,2}}\right.
\\[1ex]
&\left. +\left( -M_{111}+3\,M_{122}\right) f_{{1,2,0,0,0,3}}+\left(
-M_{111}+3\,M_{122}\right) f_{{0,1,0,3,1,2}}\right.
\\[1ex]
&\left. +\left( -3\,M_{111}+9\,M_{122}\right) f_{{1,2,0,2,0,3}}+\left(
-3\,M_{233}+M_{222}\right) f_{{1,2,0,0,1,0}}\right.
\\[1ex]
&\left. +\left( 3\,M_{223}-M_{333}\right) f_{{0,3,1,0,0,0}}+\left(
-2\,M_{223}+M_{333}\right) f_{{0,1,1,0,0,0}}\right.
\\[1ex]
&\left. -M_{133}f_{{1,0,0,0,0,1}}+M_{233}f_{{1,0,0,0,1,0}}+\left(
3\,M_{233}-3\,M_{112}\right) f_{{0,3,0,3,0,1}}\right.
\\[1ex]
&\left. +\left( -3\,M_{233}+3\,M_{222}-6\,M_{112}\right) f_{{0,3,0,1,0,1}%
}+\left( -3\,M_{233}+3\,M_{112}\right) f_{{0,1,0,3,0,1}}\right.
\\[1ex]
%\end{split}
%\end{equation*}%
%\begin{equation}
%\begin{split}
%\quad \quad
&\left. +\left( 3\,M_{112}-M_{222}\right) f_{{0,3,0,3,0,3}}+\left(
-3\,M_{122}+M_{111}\right) f_{{1,0,0,2,0,3}}\right.
\\[1ex]
&\left. +\left( -3\,M_{113}+M_{333}\right) f_{{0,3,1,2,0,0}}+\left( M_{{222}%
}-3\,M_{112}\right) f_{{0,1,0,3,0,3}}\right.
\\[1ex]
&\left. +\left( 9\,M_{112}-3\,M_{222}\right) f_{{0,3,0,1,0,3}}+\left(
3\,M_{113}-3\,M_{{223}}\right) f_{{0,3,1,0,0,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{233}+4\,M_{{112}}-2\,M_{222}\right) f_{{0,1,0,1,0,1}%
}+\left( -6\,M_{112}+2\,M_{222}\right) f_{{0,1,0,1,0,3}}\right.
\\[1ex]
&\left. +\left( 2\,M_{223}-2\,M_{{113}}\right) f_{{0,1,1,0,0,2}}+\left(
M_{133}+2\,M_{122}-M_{{111}}\right) f_{{1,0,0,2,0,1}}\right.
\\[1ex]
&\left. +\left( M_{112}-M_{233}\right) f_{{1,0,0,2,1,0}}+\left(
3\,M_{122}-3\,M_{133}\right) f_{{0,3,0,1,1,0}}\right.
\\[1ex]
&\left. +\left( -2\,M_{122}+3\,M_{133}\right) f_{{0,1,0,1,1,0}}+\left(
-3\,M_{133}+M_{111}\right) f_{{0,1,0,3,1,0}}\right.
\\[1ex]
&\left. +\left( 3\,M_{133}-M_{111}\right) f_{{0,3,0,3,1,0}}+\left(
-M_{333}+3\,M_{113}\right) f_{{0,1,1,2,0,0}}\right.
\\[1ex]
&\left. +2\,M_{123}f_{{1,0,1,1,0,0}}+\left( 9\,M_{112}-3\,M_{{222}}\right)
f_{{1,2,0,2,1,2}}\right.
\\[1ex]
&\left. +\left( 6\,M_{113}-6\,M_{{223}}\right) f_{{1,2,1,1,1,1}%
}+12\,M_{123}f_{{1,2,1,1,0,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{111}-9\,M_{122}\right) f_{{0,3,0,1,1,2}}-6\,M_{{123}%
}f_{{0,3,1,0,1,1}}\right.
\\[1ex]
&\left. +4\,M_{123}f_{{0,1,1,0,1,1}}+\left( 3\,M_{{113}}-3\,M_{223}\right)
f_{{0,3,1,2,0,2}}\right)   \label{eq151}
\end{split}
\end{equation}%
\begin{equation} % \begin{equation*}
\begin{split}
B_{\eta }=& \left( \dfrac{R_{E}}{R}\right) ^{3}\left( M_{1}f_{{0,0,1,0,1,0}%
}+M_{2}f_{{0,0,1,0,0,1}}-M_{3}f_{{0,0,0,1,0,0}}\right)
\\[1ex]
&-2\left( \dfrac{R_{E}}{R}\right) ^{4}\left( \left( M_{22}-M_{11}\right) f_{{%
0,1,1,0,1,1}}-2\,M_{12}f_{{0,1,1,0,0,2}}+M_{13}f_{{0,1,0,1,0,1}}\right.
\\[1ex]
&\left. +M_{12}f_{{0,1,1,0,0,0}}-M_{23}f_{{0,1,0,1,1,0}}+\left(
M_{33}-M_{11}\right) f_{{1,0,1,1,0,0}}\right.
\\[1ex]
%\end{split}
%\end{equation*}%
%\begin{equation*}
%\begin{split}
&\left. +\left( M_{11}-M_{22}\right) f_{{1,0,1,1,0,2}}-2\,M_{12}f_{{%
1,0,1,1,1,1}}+2\,M_{13}f_{{1,0,0,2,1,0}}\right.
\\[1ex]
&\left. +2\,M_{23}f_{{1,0,0,2,0,1}}-M_{13}f_{{1,0,0,0,1,0}}-M_{23}f_{{%
1,0,0,0,0,1}}\right)
\\[1ex]
&-3\left( \dfrac{R_{E}}{R}\right) ^{5}\left( \left(
-3\,M_{112}+3\,M_{233}\right) f_{{0,0,1,2,0,1}}-4\,M_{123}f_{{1,1,0,2,0,0}%
}\right.
\\[1ex]
&\left. +\left( 3\,M_{112}-M_{222}\right) f_{{0,0,1,2,0,3}}+\left(
M_{222}-3\,M_{112}\right) f_{{0,2,1,0,0,3}}\right.
\\[1ex]
&\left. +\left( -M_{222}+2\,M_{112}+M_{233}\right) f_{{0,2,1,0,0,1}}+\left(
M_{133}-M_{122}\right) f_{{0,2,1,0,1,0}}\right.
\\[1ex]
&\left. +\left( 3\,M_{113}-3\,M_{223}\right) f_{{0,2,0,3,0,2}}+\left(
M_{223}-M_{113}\right) f_{{0,2,0,1,0,2}}\right.
\\[1ex]
&\left. +\left( M_{111}-3\,M_{122}\right) f_{{0,0,1,2,1,2}}+\left(
M_{111}-3\,M_{133}\right) f_{{0,2,1,2,1,0}}\right.
\\[1ex]
&\left. +\left( M_{222}-3\,M_{112}\right) f_{{0,2,1,2,0,3}}+\left(
-M_{111}+3\,M_{122}\right) f_{{0,2,1,2,1,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{112}-3\,M_{233}\right) f_{{0,2,1,2,0,1}}-4\,M_{123}f_{{%
0,0,0,1,1,1}}\right.
\\[1ex]
&\left. +\left( 2\,M_{112}-2\,M_{233}\right) f_{{1,1,1,1,1,0}}+\left(
-2\,M_{111}+2\,M_{133}+4\,M_{122}\right) f_{{1,1,1,1,0,1}}\right.
\\[1ex]
&\left. +\left( 2\,M_{111}-6\,M_{122}\right) f_{{1,1,1,1,0,3}}+\left(
-6\,M_{112}+2\,M_{222}\right) f_{{1,1,1,1,1,2}}\right.
\\[1ex]
&\left. +\left( 4\,M_{113}-4\,M_{223}\right) f_{{1,1,0,2,1,1}}-4\,M_{123}f_{{%
1,1,0,0,0,2}}+8\,M_{123}f_{{1,1,0,2,0,2}}\right.
\\[1ex]
&\left. +\left( -M_{111}+3\,M_{122}\right) f_{{0,2,1,0,1,2}}+2\,M_{123}f_{{%
0,2,0,1,1,1}}\right.
\\[1ex]
&\left. +\left( -3\,M_{113}+M_{333}\right) f_{{0,2,0,3,0,0}}+\left(
-M_{333}+M_{223}+2\,M_{113}\right) f_{{0,2,0,1,0,0}}\right.
\\[1ex]
&\left. -M_{233}f_{{0,0,1,0,0,1}}-M_{133}f_{{0,0,1,0,1,0}}+\left(
-2\,M_{113}+2\,M_{223}\right) f_{{1,1,0,0,1,1}}\right.
\\[1ex]
&\left. +\left( 2\,M_{113}-2\,M_{223}\right) f_{{0,0,0,1,0,2}}+\left(
-3\,M_{113}+3\,M_{223}\right) f_{{0,0,0,3,0,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{113}-M_{333}\right) f_{{0,0,0,3,0,0}}+\left(
M_{333}-2\,M_{113}\right) f_{{0,0,0,1,0,0}}\right.
\\[1ex]
%\end{split}
%\end{equation*}%
%\begin{equation}
%\begin{split}
&\left. +6\,M_{123}f_{{0,0,0,3,1,1}}+\left( 3\,M_{133}-M_{111}\right) f_{{%
0,0,1,2,1,0}}\right. % \label{eq152}
\\[1ex]
&\left. -6\,M_{123}f_{{0,2,0,3,1,1}}+2\,M_{123}f_{{1,1,0,0,0,0}}\right)
\end{split}
\end{equation}%
\begin{equation*}
\begin{split}
B_{\zeta }=& 2\left( \dfrac{R_{E}}{R}\right) ^{3}\left( M_{1}f_{{0,1,0,0,0,1}%
}-M_{2}f_{{0,1,0,0,1,0}}+M_{1}f_{{1,0,0,1,1,0}}+M_{2}f_{{1,0,0,1,0,1}%
}\right.
\\[1ex]
&\left. +M_{3}f_{{1,0,1,0,0,0}}\right) +3\left( \dfrac{R_{E}}{R}\right)
^{4}\left( M_{33}f_{{0,0,0,0,0,0}}+2\,M_{13}f_{{0,0,1,1,1,0}}\right.
\\[1ex]
&\left. +\left( -M_{33}+M_{11}\right) f_{{0,0,0,2,0,0}}+2\,M_{23}f_{{%
0,0,1,1,0,1}}-2\,M_{12}f_{{0,2,0,0,1,1}}\right.
\\[1ex]
&\left. +\left( M_{22}-M_{11}\right) f_{{0,0,0,2,0,2}}+2\,M_{12}f_{{%
0,0,0,2,1,1}}+\left( -M_{11}+M_{33}\right) f_{{0,2,0,2,0,0}}\right.
\\[1ex]
&\left. +\left( M_{11}-M_{22}\right) f_{{0,2,0,0,0,2}}-2\,M_{12}f_{{%
1,1,0,1,0,0}}+\left( M_{11}-M_{22}\right) f_{{0,2,0,2,0,2}}\right.
\\[1ex]
&\left. +\left( 2\,M_{11}-2\,M_{22}\right) f_{{1,1,0,1,1,1}}+4\,M_{12}f_{{%
1,1,0,1,0,2}}+2\,M_{13}f_{{1,1,1,0,0,1}}\right.
\\[1ex]
&\left. -2\,M_{23}f_{{1,1,1,0,1,0}}-2\,M_{12}f_{{0,2,0,2,1,1}}-2\,M_{13}f_{{%
0,2,1,1,1,0}}-2\,M_{23}f_{{0,2,1,1,0,1}}\right.
\\[1ex]
&\left. +\left( M_{22}-M_{33}\right) f_{{0,2,0,0,0,0}}\right) +4\left( \dfrac{%
R_{E}}{R}\right) ^{5}\left( \left( 3\,M_{122}-3\,M_{133}\right) f_{{%
0,3,0,0,0,1}}\right.
\\[1ex]
&\left. +\left( 3\,M_{222}-9\,M_{112}\right) f_{{0,3,0,2,1,2}}+\left(
-M_{111}+3\,M_{133}\right) f_{{1,2,0,3,1,0}}\right.
\\[1ex]
&\left. +\left( -M_{111}+3\,M_{122}\right) f_{{1,0,0,3,1,2}}-3\,M_{233}f_{{%
0,1,0,0,1,0}}-6\,M_{123}f_{{1,2,1,0,1,1}}\right.
\\[1ex]
&\left. +6\,M_{123}f_{{0,3,1,1,0,0}}-6\,M_{123}f_{{1,2,1,2,1,1}}+\left(
-3\,M_{222}+9\,M_{112}\right) f_{{1,2,0,1,0,3}}\right.
\\[1ex]
&\left. +\left( M_{111}-3\,M_{122}\right) f_{{1,2,0,3,1,2}}+\left(
-3\,M_{133}+3\,M_{122}\right) f_{{1,2,0,1,1,0}}\right.
\\[1ex]
&\left. +\left( -3\,M_{233}+3\,M_{222}-6\,M_{112}\right) f_{{1,2,0,1,0,1}%
}-12\,M_{123}f_{{0,3,1,1,0,2}}\right.
\\[1ex]
&\left. +\left( 6\,M_{223}-6\,M_{113}\right) f_{{0,3,1,1,1,1}}+3\,M_{233}f_{{%
1,0,0,1,0,1}}\right.
\\[1ex]
&\left. +\left( 9\,M_{122}-3\,M_{111}\right) f_{{0,1,0,2,0,3}}+\left(
-3\,M_{112}+M_{222}\right) f_{{0,3,0,0,1,2}}\right.
\\[1ex]
&\left. +3\,M_{133}f_{{1,0,0,1,1,0}}+6\,M_{123}f_{{1,0,1,2,1,1}}+\left(
M_{111}-3\,M_{122}\right) f_{{0,3,0,0,0,3}}\right.
\\[1ex]
&\left. +\left( -6\,M_{223}+6\,M_{113}\right) f_{{0,1,1,1,1,1}}+12\,M_{123}f_{%
{0,1,1,1,0,2}}\right.
\\[1ex]
&\left. +\left( -M_{333}+3\,M_{113}\right) f_{{1,0,1,2,0,0}}+\left(
3\,M_{233}-3\,M_{112}\right) f_{{1,2,0,3,0,1}}\right.
\\[1ex]
&\left. +\left( 3\,M_{113}-3\,M_{223}\right) f_{{1,2,1,2,0,2}}+\left(
-3\,M_{233}+3\,M_{112}\right) f_{{0,3,0,2,1,0}}\right.
\\[1ex]
&\left. +\left( -3\,M_{133}+M_{111}\right) f_{{1,0,0,3,1,0}}+\left(
M_{333}-3\,M_{113}\right) f_{{1,2,1,2,0,0}}\right.
\\[1ex]
&\left. +\left( 3\,M_{233}-3\,M_{112}\right) f_{{0,1,0,2,1,0}}+\left(
3\,M_{223}-M_{333}\right) f_{{1,2,1,0,0,0}}\right.
\\[1ex]
&\left. +3\,M_{133}f_{{0,1,0,0,0,1}}+\left( 3\,M_{113}-3\,M_{223}\right) f_{{%
1,2,1,0,0,2}}\right.  \label{eq153}
\\[1ex]
&\left. +\left( 3\,M_{233}-M_{222}\right) f_{{0,3,0,0,1,0}}+\left(
-3\,M_{233}+3\,M_{112}\right) f_{{1,0,0,3,0,1}}\right.
\\[1ex]
&\left. +\left( 3\,M_{111}-9\,M_{122}\right) f_{{1,2,0,1,1,2}}+\left(
-3\,M_{112}+M_{222}\right) f_{{1,0,0,3,0,3}}\right.
\\[1ex]
\end{split}
\end{equation*}%
\begin{equation}
\begin{split}
&\left. +\left( -3\,M_{222}+9\,M_{112}\right) f_{{0,1,0,2,1,2}}+\left(
-3\,M_{133}+3\,M_{111}-6\,M_{122}\right) f_{{0,1,0,2,0,1}}\right.
\\[1ex]
&\left. -6\,M_{123}f_{{0,1,1,1,0,0}}+\left( 3\,M_{223}-3\,M_{113}\right) f_{{%
1,0,1,2,0,2}}\right.
\\[1ex]
&\left. +\left( 3\,M_{111}-9\,M_{122}\right) f_{{0,3,0,2,0,3}}+\left(
3\,M_{112}-M_{222}\right) f_{{1,2,0,3,0,3}}\right.
\\[1ex]
&\left. +\left( 6\,M_{122}-3\,M_{111}+3\,M_{133}\right) f_{{0,3,0,2,0,1}%
}+M_{333}f_{{1,0,1,0,0,0}}\right) ,
\end{split}
\end{equation}%
where $f_{i_1,i_2,i_3,i_4,i_5,i_6}=\sin^{i_1}(u)\cos^{i_2}(u)\sin^{i_3}(i)\cos^{i_4}(i)\sin^{i_5}(\phi)\cos^{i_6}(\phi).$

%��� ���������� �������� ������� ${\vec M}_{\text �}$ �� ��� $x, y, z$ ����������
%��������� (\ref{eq2}) � ����
%$$\vec M_{\text �}= Q(x_0\vec i+y_0\vec j+z_0\vec k)\times
%\left[\begin{pmatrix} \alpha_1 & \beta_1 & \gamma_1  \\
%               \alpha_2 & \beta_2 & \gamma_2  \\
%               \alpha_3 & \beta_3 & \gamma_3
%               \end{pmatrix}
%\begin{pmatrix} v_{C\eta}B_\zeta\\
%         -v_{C\xi}B_\zeta \\
%         v_{C\xi}B_\eta  -v_{C\eta}B_\xi\end{pmatrix}\right].
%$$
It is shown in \cite{KRes2003} that the moment ${\vec M}_{L}$ of Lorentz forces can considerably exceed gravitational and other perturbing moments and can be used as a restoring moment in the system of spacecraft orientation.

\section{Control synthesis}
Let the program orientation of the spacecraft in the orbital coordinate system is prescribed by some value of the matrix of direction cosines $\mathbf A=\mathbf A_0=$const. In particular, the matrix $\mathbf A_0$ can be unitary. Substituting $\mathbf A=\mathbf A_0$ into expression (\ref{eq2}), we obtain the values of the vector $\vec{P}$ such that ${\vec M}_{L}$ becomes zero in the program orientation of the spacecraft, i.e., is a restoring moment in the neighborhood of the orientation $\mathbf A_0$. It is obvious that the vector $\vec{P}$ must comply with the condition $\vec{P}=k_L\vec{T}_0\, ,$ where $k_L=k_L(t)$ is an arbitrary scalar function, $\vec T_0=\mathbf A_{0}^{\top} (\vec v_C\times\vec
B)$. Components of the vector $\vec P$ are specified by the formulas
\begin{equation}\label{eq3}
\begin{aligned}
 P_x &=k_L(t)[\alpha_{10}v_{C\eta}B_{\zeta}-
 \beta_{10}v_{C\xi}B_{\zeta}+
 \gamma_{10}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi})], \\
P_y &=k_L(t)[\alpha_{20}v_{C\eta}B_{\zeta}-
 \beta_{20}v_{C\xi}B_{\zeta}+
 \gamma_{20}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi})], \\ %%%%%%%%%%%%%%%%%%%%%%%%%%% (5)
P_z &=k_L(t)[\alpha_{30}v_{C\eta}B_{\zeta}-
 \beta_{30}v_{C\xi}B_{\zeta}+
 \gamma_{30}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi})]. \\
\end{aligned}
\end{equation}
Hence, if the coordinates of the center of charge of the spacecraft will change according to law (\ref{eq3}),
then the moment ${\vec M}_L$ will be restoring in the neighborhood of the prescribed position and can be used for maintaining the prescribed orientation of the spacecraft. Equalities (\ref{eq3}) can be considered as a control law for the position of the center of charge for performing the prescribed orientation of the spacecraft.
Further in the paper we shall assume for distinctness that $k_L(t)=k_L=\,$ const.

The dependencies of the maximum values of $|\vec \rho_{0}|$ (and, hence, of $|\vec P|$) for different parameters of orbits were computed in \cite{Vest2006}; also, there was demonstrated that $|\vec \rho_{0}|$ assumes values that are small comparing with the sizes of the spacecraft or comparable with them. Therefore, the practical realization of the spacecraft attitude orientation using the moment of Lorentz forces does not present serious difficulties.

We shall note the following functional particular features of the systems using ${\vec M}_{L}$ for the spacecraft attitude stabilization. The first of them is that ${\vec M}_{L}$ is orthogonal to the vector $\vec T = \mathbf A^\top(\vec v_C\times\vec B)$ and, therefore, it is impossible to create a control moment of Lorentz forces directed along the vector $\vec T$.
The second one is that $|{\vec M}_{L}|$ is proportional to $|\vec T|$. So the use of ${\vec M}_{L}$ for the spacecraft attitude stabilization is most effective for spacecrafts on low inclined orbits because the averaged by orbital period values of $|\vec T|$ for these orbits are larger than the same values for the orbits with large inclinations.

It should be emphasized that the program control of vector $\vec P$ by ensuring the validity of equalities (\ref{eq3}) is based on the use of predetermined laws of variation for vectors $\vec v_C$ and $\vec B$ as the functions of time and so it does not require the measuring of any physical properties on the board of spacecraft during the stabilization process.

Let us consider the moment of magnetic interaction $\vec M_M$. Neglecting the spacecraft magnetization in geomagnetic field we define $\vec M_M$ according to the formula \cite{Kov.1975, Alp.1978}:
\begin{equation} \label{eq69}
\vec M_M= \vec I\times\mathbf
A^{\top}\vec B,
\end{equation}
where vector $\vec B$ is defined by its projections on the axes of the orbital coordinate system $C\xi\eta\zeta$ (Fig.\ref{Fig1}), the eigen magnetic moment vector $\vec I$ is defined by its projections on the spacecraft principal central axes of inertia $Cxyz$ (Fig.\ref{Fig2}). As is known, the moment of magnetic interaction can be used for the stabilization of the spacecraft angular position. Indeed, the moment of magnetic interaction is restoring in vicinity of equilibrium position $\mathbf A_{0}$ by analogy with the Lorentz moment if vector $\vec I$ changes in program way according to the law $ \vec I = k_{M}\vec B_{0}\, ,$ where $\vec B_0=\mathbf A_{0}^{\top}\vec B$ and $k_M$ is some scalar coefficient of proportionality.
In scalar form vector $\vec I$ becomes
\begin{eqnarray}\label{eq128}
I_{x}=k_{M}(t)\bigl[\alpha_{10}B_{\xi}+\beta_{10}B_{\eta}+\gamma_{10}B_{\zeta}\bigr]\nonumber\\
I_{y}=k_{M}(t)\bigl[\alpha_{20}B_{\xi}+\beta_{20}B_{\eta}+\gamma_{20}B_{\zeta}\bigr]\\
I_{z}=k_{M}(t)\bigl[\alpha_{30}B_{\xi}+\beta_{30}B_{\eta}+\gamma_{30}B_{\zeta}\bigr]\nonumber
\end{eqnarray}

Equalities (\ref{eq128}) may be considered as the control law for the spacecraft eigen magnetic moment vector in execution of the orientation defined by matrix $\mathbf A_\mathrm{0}$. Further in the paper we shall assume for distinctness that  $k_M(t)=k_M=\,$const.

The magnetic control systems can be successfully used on long time functioning spacecrafts as they are sufficiently simple, possesses high reliability and don't need fuel consumption \cite{Kov.1975, Alp.1978}.

At the same time the magnetic control systems possesses some specific features which limits their possibilities.
Firstly, it follows from (\ref{eq69}) that the moment of magnetic interaction is orthogonal to vector $\vec B$ and therefore it is impossible to create the control moment directed along vector $\vec B$.
Secondly, $|\vec M_M|$ is proportional to $|\vec B|$. Therefore the magnetic control systems are used mainly on the orbits with large inclinations since the averaged by orbital period values of $|\vec B|$ for these orbits are larger than the same values for the orbits with small inclinations.

It should be emphasized that the program control of vector $\vec I$ by ensuring the validity of equalities (\ref{eq128}) is based on the use of predetermined law of variation for vector $\vec B$ as the functions of time and so it does not require the measuring of vector $\vec B$ on the board of spacecraft during the stabilization process. Let us remind that the octupole approximation is used for vector $\vec B$, describing the geomagnetic field.

After comparing the mentioned above functional features of magnetic control systems and systems using the moment ${\vec M}_{L}$ and shortcomings that follow from them, it is easy to see that these shortcomings disappear when a unified electrodynamic control system for spacecraft attitude orientation is created using simultaneously both restoring moments:
\begin{equation} \label{eq12}
{\vec M}_{L}=k_{L}\vec T_{0}\times\vec T,
\end{equation}
\begin{equation} \label{eq13}
{\vec M}_{M}=k_{M}\vec B_{0}\times\vec B.
\end{equation}

Indeed, the first shortcomings consisting in the presence of the directions such that the uncontrolled rotation of the spacecraft is possible along them (the directions $\vec T$ and $\vec B$) disappear, when ${\vec M}_{L}$ and ${\vec M}_{M}$ are united in one electrodynamic control system. This is explained by the fact that the vectors $\vec T$ and $\vec B$ are always orthogonal. The second of the mentioned above shortcomings are naturally compensate each other and therefore the orbit inclination is not the limiting factor for the control possibilities if the satellite is supplied with electrodynamic control system.


\section{The differential equations of motion}

Let us consider the question of spacecraft attitude stabilization in the direct equilibrium position in the orbital coordinate system, i.e., in such position when the axes $x,\,y,\,z$ coincide with the axes $\xi,\,\eta,\,\zeta$ and hence $\varphi =\theta =\psi =0$. At the same time
\begin{equation} \label{eq15}
\alpha_1=\beta_2=\gamma_3=1,\; \omega_x=\omega_z=0,\omega_y=\omega_0 .
\end{equation}

Substituting (\ref{eq15}) into (\ref{eq3}) we obtain the law of controlled changing for the center of charge coordinates for the case of the direct equilibrium position in the form:
\begin{equation}\label{eq14}
x_0=k_{L}\,v_{C\eta}B_{\zeta},\ \ \ \ y_0=-k_{L}\,v_{C\xi}B_{\zeta},\ \
\ \ \ z_0=k_{L}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi}) .
\end{equation}

On the base of expression (\ref{eq12}) the projections $M_{L x},\ M_{L y}, \ M_{L z}$ of the restoring moment $\vec{M}_{L}$ under conditions of control (\ref{eq14}) have the following form:
\begin{equation}
\begin{split}
 \label{eq131}
    M_{L x}=& Qk_{L}\Bigl[B_{\zeta}((\beta_2-\gamma_3)v_{C\xi}-\alpha_2v_{C\eta})(v_{C\xi}B_{\eta}
    -v_{C\eta}B_{\xi})\Bigr.\\
    & \Bigl. -\gamma_2(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi})^2+
    v_{C\xi}B_{\zeta}^2(\beta_3v_{C\xi}-\alpha_3v_{C\eta})\Bigr],
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\label{eq132}
M_{L y} = & Qk_{L}\Bigl[B_{\zeta}((\alpha_1-\gamma_3)v_{C\eta}-\beta_1v_{C\xi})(v_{C\xi}B_{\eta}
-v_{C\eta}B_{\xi})\Bigr.\\
&\Bigl. +\gamma_1(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi})^2+
    v_{C\eta}B_{\zeta}^2(\beta_3v_{C\xi}-\alpha_3v_{C\eta})\Bigr],
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\label{eq133}
M_{L z}= & Qk_{L}B_{\zeta}\Bigl[
    (\gamma_1v_{C\xi}+\gamma_2v_{C\eta})(v_{C\xi}B_{\eta}
    -v_{C\eta}B_{\xi})\Bigr. \\
    & \Bigl.+ B_{\zeta}(v_{C\xi}v_{C\eta}(\alpha_1-\beta_2)+
    (v_{C\eta}^2\alpha_2-\beta_1v_{C\xi}^2))\Bigr].
\end{split}
\end{equation}

Similarly substituting (\ref{eq15}) into (\ref{eq128}) we obtain the law of controlled changing for the coordinates of vector $\vec I$ in the form:
\begin{equation}\label{eq129}
I_{x}=k_{M}\,B_{\xi},\ \ \ \ I_{y}=k_{M}\,B_{\eta},\ \
\ \ \ I_{z}=k_{M}\, B_{\zeta}.
\end{equation}
On the basis of expression (\ref{eq13}) the projections $M_{M x},\ M_{M y}, \ M_{M z}$ of the restoring moment $\vec{M}_{M}$
in conditions of control (\ref{eq129}) have the following form:
\begin{eqnarray*} %{c}
M_{M x}=k_{{M}} \left[ B_{{\eta }} \left( \alpha_{{3}}B_{{\xi }}+\beta_{{3}}B_{{\eta }}+
\gamma_{{3}}B_{{\zeta }} \right) -B_{{\zeta }} \left( \alpha_{{2}}
B_{{\xi }}+\beta_{{2}}B_{{\eta }}+\gamma_{{2}}B_{{\zeta }} \right)
 \right],\nonumber\\
%\end{equation*}%\\[-1.5truecm]
%\begin{equation*} %{c}
M_{M y}=k_{{M}} \left[ B_{{\zeta }} \left( \alpha_{{1}}B_{{\xi }}+\beta_{{1}}B_{{\eta }}+
\gamma_{{1}}B_{{\zeta }} \right) -B_{{\xi }} \left( \alpha_{{3}}
B_{{\xi }}+\beta_{{3}}B_{{\eta }}+\gamma_{{3}}B_{{\zeta }} \right)
 \right],\nonumber\\
%\end{equation*}%\\[-2truecm]
%\begin{equation*} %{c}
M_{M z}=k_{{M}} \left[ B_{{\xi }} \left( \alpha_{{2}}B_{{\xi }}+\beta_{{2}}B_{{\eta }}+
\gamma_{{2}}B_{{\zeta }} \right) -B_{{\eta }} \left( \alpha_{{1}}
B_{{\xi }}+\beta_{{1}}B_{{\eta }}+\gamma_{{1}}B_{{\zeta }} \right)
 \right].\nonumber
\end{eqnarray*}

Let the spacecraft attitude control system has also the damping moment $\vec{M}_D.$ Without going into details of technical realization of the oscillations damping system we shall consider the moment $\vec{M}_D$ of model type, for example proportional to the spacecraft relative angular velocity in the orbital frame: $\vec{M}_D=-\vec{H}_D\vec{\omega'},$ where $\
\vec{H}_D=\mbox{diag}(h_1,h_2,h_3), \ h_i>0 \ (i=1,2,3),\ \vec{\omega'}=\vec{\omega}-\vec{\omega_0}$, and $\vec{\omega}$ is the  absolute spacecraft angular velocity.
Then
\begin{equation}\label{eq130}
{M}_{Dx}=-h_1(\omega_x-\omega_0\beta_1),
\ {M}_{Dy}=-h_2(\omega_y-\omega_0\beta_2), \ {M}_{Dz}=-h_3(\omega_z-\omega_0\beta_3).
\end{equation}

We shall take into account the gravitational moment $\vec{M}_G$ as the largest of disturbing moments acting upon spacecraft.

Let us prove that the suggested control (\ref{eq14}), (\ref{eq129}) in presence of damping (\ref{eq130}) in the attitude control system ensures the existence and asymptotic stability of the satellite's direct equilibrium position in the case of small orbit inclinations and the total stability in other cases, i.e. solve the problem of the spacecraft attitude stabilization in the orbital frame.

Differential equations of the spacecraft attitude motion under the influence of control moments (\ref{eq12}),
(\ref{eq13}) are constructed according to the Euler-Poisson scheme:
\begin{equation}\label{eq6}
\gathered \frac d{dt}(\mathbf J{\vec\omega})+{\vec\omega}\times
(\mathbf J{\vec\omega})={\vec M}_{L}+{\vec M}_{M}+{\vec M}_{D}+{\vec M}_{G},                     \\
\frac {d{\vec\xi}_0}{dt}={\vec\xi}_0\times{\vec\omega} - \omega_0{\vec\zeta}_0, \quad
\frac {d{\vec\eta}_0}{dt}={\vec\eta}_0\times{\vec\omega}, \quad \frac
{d{\vec\zeta}_0}{dt}={\vec\zeta}_0\times{\vec\omega}+ \omega_0{\vec\xi}_0,
\endgathered
\end{equation}
where $\mathbf J=\operatorname{diag}(A,B,C)$ is the spacecraft inertia tensor in the coordinate system $Cxyz$.

While considering small spacecraft oscillations in the neighborhood of the direct equilibrium position, the assumption that the angles $\varphi,\theta,\psi$ and their time derivatives are small is true. At the same time, the moments $\vec{M}_{L}, \ \vec{M}_{M}, \ \vec{M}_{G},\ \vec{M}_{D}$ can be extended into series according to degrees of these small quantities. As a result, we obtain their projections accurate within the terms of the second order of smallness in the form
\begin{align}\label{eq134}
{M}_{L x} &
=l_{11}(t)\varphi+l_{12}(t)\theta+l_{13}(t)\psi,
&{M}_{M x} &
=b_{11}(t)\varphi+b_{12}(t)\theta+b_{13}(t)\psi,\nonumber
\\{M}_{L y} & =l_{12}(t)\varphi+l_{22}(t)\theta+l_{23}(t)\psi,
&{M}_{M y} & =b_{12}(t)\varphi+b_{22}(t)\theta+b_{23}(t)\psi,
\\{M}_{L z} & =l_{13}(t)\varphi+l_{23}(t)\theta+l_{33}(t)\psi,
&{M}_{M z} & =b_{13}(t)\varphi+b_{23}(t)\theta+b_{33}(t)\psi,\nonumber
\end{align}
\begin{align}\label{eq1341}
{M}_{G x} &= 3\omega_0^2(C-B)\varphi,&
{M}_{D x} & = -h_1\dot{\varphi},\nonumber
\\
{M}_{G y} &= 3\omega_0^2(C-A)\theta,&
{M}_{D y} & = -h_2\dot{\theta},
\\
{M}_{G z} &=0,& {M}_{D z} & = -h_3\dot{\psi},\nonumber
\end{align}
where
\[ l_{11}(t)=-Qk_{L}\left[(v_{C\eta}B_{\xi}-v_{C\xi}B_{\eta})^2+
v_{C\xi}^2B_{\zeta}^2\right],\]
\[
l_{12}(t)=l_{21}(t)=-Qk_{L}v_{C\xi}v_{C\eta}B_{\zeta}^2,\]\[
l_{13}(t)=l_{31}(t)=Qk_{L}v_{C\eta}B_{\zeta}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi}),\]
\[
l_{22}(t)=-Qk_{L}\left[(v_{C\eta}B_{\xi}-v_{C\xi}B_{\eta})^2+v_{C\eta}^2B_{\zeta}^2
\right],\]
\[
l_{23}(t)=l_{32}(t)=-Qk_{L}v_{C\xi}B_{\zeta}(v_{C\xi}B_{\eta}-v_{C\eta}B_{\xi}),\]
\[ l_{33}(t)=-Qk_{L}(v_{C\xi}^2+v_{C\eta}^2)B_{\zeta}^2 .\]

\[ b_{11}(t)=-k_{M} (B_{\eta}^2+B_\zeta^2),\ b_{12}(t)=b_{21}(t)=k_{M}B_\xi B_{\eta}\, ,\]
\[ l_{13}(t)=l_{31}(t)=k_{M}B_\xi B_{\zeta},\
b_{22}(t)=-k_{M}(B_{\xi}^2+B_\zeta^2)\, ,\]
\[
b_{23}(t)=b_{32}(t)=k_{M}B_\eta B_{\zeta},\ b_{33}(t)=-k_{M}(B_{\xi}^2+B_\eta^2)\, .\]

The Euler equations of motion (\ref{eq6}) in the matrix form will look like
\begin{equation}\label{eq7}
\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)+\vec X=0,
\end{equation}
\[
\mbox{where \ }\ \mathbf{J}=\left(
  \begin{array}{ccc}
    A&0&0 \\
    0&B&0\\
    0&0&C
  \end{array}\right),\ \  \mathbf{H}=\left(
  \begin{array}{ccc}
    h_1&0&\omega_0(A-B+C) \\
    0 & h_2& 0\\
    -\omega_0(A-B+C)&0&h_3
  \end{array}\right),
\]
{\small{\[\mathbf{M}=\left(
  \begin{array}{ccc}
    4\omega_0^2{{(B-C)}}-l_{11}(t)-b_{11}(t)\!\!\!\!\!\!\!\!\!\!&-l_{12}(t)-b_{12}(t)\!\!\!\!\!\!\!\!\!\!&-l_{13}(t)- b_{13}(t) \\
    -l_{21}(t)- b_{21}(t)\!\!\!\!\!\!\!\!\!\! & 3\omega_0^2{{(A-C)}}-l_{22}(t)- b_{22}(t)\!\!\!\!\!\!\!\!\!\!& -l_{23}(t)- b_{23}(t)\\
    -l_{31}(t)- b_{31}(t)\!\!\!\!\!\!\!\!\!\!&-l_{32}(t)- b_{32}(t)\!\!\!\!\!\!\!\!\!\!&\omega_0^2(B-A)-l_{33}(t)- b_{33}(t)
  \end{array}\!\right) ,
\]}}
$\vec X$ is a vector with the components $X_j(t,\varphi,\theta,\psi)\ \ (j=\overline{1,3}),$ nonlinearly dependent on $\varphi,\theta,\psi.$

\section{Analysis of stabilization process at the orbits of small inclinations}

The components $m_{ij}(t)$ of matrix $\mathbf{M}$, depending, generally speaking, on the small parameter $\sin i$, we represent as the sum of the terms $m_{ij}^{(0)}(t)$, independent on $\sin(i)$, and the terms $m_{ij}^{(1)}(t)$, containing small parameter $\sin i$ as multiplier. Then $l_{ij}(t)=l_{ij}^{(0)}(t)+\sin i \cdot l_{ij}^{(1)}(t)\, ,$
$b_{ij}(t)=b_{ij}^{(0)}(t)+\sin i \cdot b_{ij}^{(1)}(t)\, $. Therefore
\begin{multline}\label{eq136}
\mathbf{M}=\mathbf{M}^{(0)}+\sin i \cdot \mathbf{M}^{(1)} = \mbox{diag}(4\omega_0^2(B-C),3\omega_0^2(A-C),\omega_0^2(B-A))\\
-[(l_{ij}^{(0)}(t))+(b_{ij}^{(0)}(t))]+\sin i \cdot[-(l_{ij}^{(1)}(t))-(b_{ij}^{(1)}(t))].
\end{multline}

Consider the system of linear approximation of equations (\ref{eq7}) at $i=0$ (the spacecraft orbit is equatorial):
\begin{equation}\label{eq97}
\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}^{(0)}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)=0.
\end{equation}

In this case the equalities $v_{C\xi}=R(\omega_0-\omega_3)$, \, $v_{C\eta}=~v_{C\zeta}=~0$ are valid on the basis of (\ref{eq33}). Relying on (\ref{eq151}) -- (\ref{eq153}) we have the following projections of vector $\vec B$ on the axes of orbital frame in octupole approximation at $i=0$:
\begin{multline}\label{eq105}
B_{\xi}=\left( \frac{R_{E}}{R}\right)^3\left(g^1_1\sin u_0-h^1_1\sin u_0 \right)\\
+\sqrt {3}\left(\frac{R_{E}}{R}\right)^4 \left(g_2^2\sin( 2u_0)-h_2^2\cos (2u_0)\right)+\frac{\sqrt {6}}{4}\left(\frac{R_{E}}{R}\right)^5\left(
\sqrt {15}g_3^3\sin(3u_0 )
\right.\\
\left.
-g_3^1\sin u_0
+h_3^1\cos u_0
-\sqrt {15}h_3^3\cos(3u_0)
\right),
\end{multline}
\begin{multline}\label{eq106}
 B_{\eta}=-\left( \frac{R_{E}}{R}\right)^3g^1_0
 -\sqrt {3}\left(\frac{R_{E}}{R}\right)^4
 \left(
g_2^1\cos u_0
+h_2^1\sin u_0
 \right)\\
 -\frac{1}{2}\left(\frac{R_{E}}{R}\right)^5
 \left(
\sqrt {15}g_3^2\cos(2u_0)
 +\sqrt {15}h_3^2\sin (2u_0)
 -3 g_3^0
 \right),
\end{multline}
\begin{multline}\label{eq107}
B_{\zeta}=2\left( \frac{R_{E}}{R}\right)^3
\left(
g_1^1\cos  u_0
+h_1^1\sin u_0
\right)+
\frac{3}{2}\left( \frac{R_{E}}{R}\right)^4
\left(
\sqrt {3}g_2^2\cos( 2u_0)
+\right.\\
\left.
+\sqrt {3}h_2^2\sin ( 2u_0)
-g_2^0
\right)+\sqrt {2}\left( \frac{R_{E}}{R}\right)^5
\left(
\sqrt {5}g_3^3\cos( 3u_0 )
\right.\\
\left.
+\sqrt {5}h_3^3\sin( 3u_0 )
-\sqrt {3}h_3^1\sin ( u_0 )
-\sqrt {3}g_3^1\cos( u_0 ),
\right)
\end{multline}
where $u_0=u-\phi$.
Then
\begin{align*}
l_{11}^{(0)}(t)&=-Qk_{L}v^2_{C\xi}\left(
B_{\eta}^2+B_{\zeta}^2\right), \
&l_{12}^{(0)}(t)&=l_{21}^{(0)}(t)=l_{13}^{(0)}(t)=l_{31}^{(0)}(t)=0,\\
l_{22}^{(0)}(t)&=-Qk_{L}v^2_{C\xi} B_{\eta}^2,\
&l_{23}^{(0)}(t)&=l_{32}^{(0)}(t)=-Qk_{L}v_{C\xi}^2
B_{\eta}B_{\zeta},\\
l_{33}^{(0)}(t)&=-Qk_{L}v_{C\xi}^2
B_{\zeta}^2.
\end{align*}

Let us represent $l_{ij}^{(0)}(t)$ and $b_{ij}^{(0)}(t)$ in the form
\begin{equation}\label{eq137}
l_{ij}^{(0)}(t)=l_{ij \,
cp}^{(0)}+\widetilde{l}_{ij}^{(0)}(t), \quad
b_{ij}^{(0)}(t)=b_{ij \,
cp}^{(0)}+\widetilde{b}_{ij}^{(0)}(t),
\end{equation}
where $l_{i j \, cp }^{(0) }={\langle l_{ij}^{(0)}(t)\rangle}_t$ and $b_{i j \, cp }^{
(0) }={\langle l_{ij}^{(0)}(t)\rangle}_t$ are the averaged values of functions $l_{ij}^{(0)}(t)$ and $b_{ij}^{(0)}(t)$ over  $t$ accordingly. Then
\begin{equation}\label{eq89}
\begin{array}{ll}
l_{11\, cp}^{(0)}=-Qk_{L}v_{C\xi}^2({\langle B_{\eta}^2\rangle}_t+{\langle B_{\zeta}^2\rangle}_t),&
l_{22\, cp}^{(0)} =-Qk_{L}v_{C\xi}^2{\langle B_{\eta}^2\rangle}_t,\\
l_{23\, cp}^{(0)}=l_{32\, cp}^{(0)}=-Qk_{L}v_{C\xi}^2{\langle B_{\eta}B_{\zeta}\rangle}_t,&
l_{33\, cp}^{(0)}=-Qk_{L}v_{C\xi}^2{\langle B_{\zeta}^2\rangle}_t, \\
\multicolumn{2}{c}{ l_{12\, cp}^{(0)}=l_{21\, cp}^{(0)}=l_{13\, cp}^{(0)}=l_{31\, cp}^{(0)}=0.}
\end{array}
\end{equation}
\begin{equation}\label{eq1371}
\begin{array}{ll}
b_{11cp}^{(0)}=-k_{M}({\langle B_{\eta}^2\rangle}_t+{\langle B_{\zeta}^2\rangle}_t),&
b_{12cp}^{(0)}=b_{21cp}^{(0)}=k_{M}{\langle B_{\xi}B_{\eta}\rangle }_t,\\ b_{13cp}^{(0)}=b_{31cp}^{(0)}=k_{M}{\langle
B_{\xi}B_{\zeta}\rangle }_t ,&
b_{22cp}^{(0)}=-k_{M}({\langle B_{\xi}^2\rangle}_t+{\langle B_{\zeta}^2\rangle}_t),\\
b_{23cp}^{(0)}=b_{32cp}^{(0)}=k_{M}{\langle B_{\eta}B_{\zeta}\rangle }_t,& b_{33cp}^{(0)}=-k_{M}({\langle
B_{\xi}^2\rangle}_t+{\langle B_{\eta}^2\rangle}_t) \, .
\end{array}
\end{equation}
Let us calculate the necessary time averaged functions:
\begin{eqnarray*}
{\langle B_{\eta }^{2}\rangle }_{t} &=&\left( \dfrac{{{R_{E}}}}{{R}}\right)
^{6}{g_{1}^{0}}^{2}+\dfrac{3}{2}\left( \dfrac{{{R_{E}}}}{{R}}\right)
^{8}\left( {g_{2}^{1}}^{2}+{h_{2}^{1}}^{2}-2g_{1}^{0}g_{3}^{0}\right)  \\
&&+\dfrac{3}{8}\left( \dfrac{{{R_{E}}}}{{R}}\right) ^{10}\left( 5\left( {%
g_{3}^{2}}^{2}+{h_{3}^{2}}^{2}\right) +6\,{g_{3}^{0}}^{2}\right) ,
\end{eqnarray*}%
\begin{multline*}%\label{eq91}
{\langle
B_{\zeta}^2\rangle}_t=
2\left(\dfrac {{{R_{E}}}}{{R}}\right)^{6}
\left( {g_1^1}^{2}+{h_1^1}^{2} \right) \\
-\dfrac{1}{8}
\left(\dfrac {{{R_{E}}}}{{R}}\right)^{8}
\left( 16\sqrt {6} \left( g_1^1 g_3^1 +h_1^1 h_3^1\right)-27
\left({g_2^2}^{2}+{h_2^2}^{2}\right)-18{g_2^0}^{2}  \right) \\
+\left(\dfrac {{{R_{E}}}}{{R}}\right)^{10}
 \left( 5\left({g_3^3}^{2}+{h_3^3}^{2}\right)+3\left({g_3^1}^{
2}+{h_3^1}^{2}\right) \right),
\end{multline*}
\begin{multline*}%\label{eq92}
{\langle B_{\eta}B_{\zeta}\rangle}_t=
\dfrac{1}{2}
\left(\dfrac {{{R_{E}}}}{{R}}\right)^{7}
\left( 3g_1^0 g_2^0-2\sqrt {3}\left(h_2^1 h_1^1
+g_2^1g_1^1\right) \right) \\
+\dfrac{3}{8}
\left(\dfrac {{{R_{E}}}}{{R}}\right)^{9}
\left( 4\sqrt {2}\left(g_2^1g_3^1+h_2^1h_3^1\right)
-3\sqrt {5}(g_3^2 g_2^2+h_3^2h_2^2)-6g_3^0g_2^0
 \right),
\end{multline*}
\begin{multline*}%\label{}
{\langle
B_{\xi}^2\rangle}_t=
\dfrac{1}{2}\left(\dfrac {{{R_{E}}}}{{R}}\right)^{6}
\left({g_1^1}^{2}+{h_1^1}^{2}\right)
\\
-\dfrac{\sqrt{6}}{4}
\left(\dfrac {{{R_{E}}}}{{R}}\right)^{8}
\left(g_3^1g_1^1+ h_3^1h_1^1-\sqrt{6}\,({g_2^2}^{2}
+{h_2^2}^{2}) \right)
\\
+\dfrac{3}{16}
\left(\dfrac {{{R_{E}}}}{{R}}\right)^{10}
\left( 15\,({g_3^3}^{2}+{h_3^3}^{2})+{g_3^1}^{2}+{h_3^1}^
{2} \right),
\end{multline*}
\begin{eqnarray*}
{\langle B_{\xi }B_{\eta }\rangle }_{t} &=&-\dfrac{\sqrt{3}}{2}\left( \dfrac{%
{{R_{E}}}}{{R}}\right) ^{7}\left(
g_{1}^{1}h_{2}^{1}-h_{1}^{1}g_{2}^{1}\right) -\dfrac{3\sqrt{2}}{8}\left(
\dfrac{{{R_{E}}}}{{R}}\right) ^{9}  \\
&&\times \left( \sqrt{10}%
(g_{2}^{2}h_{3}^{2}-h_{2}^{2}g_{3}^{2})-g_{3}^{1}h_{2}^{1}+h_{3}^{1}g_{2}^{1}\right) ,
\end{eqnarray*}%
\begin{equation*}
{\langle B_{\xi }B_{\zeta }\rangle }_{t}=-\dfrac{\sqrt{6}}{4}\left( \dfrac{{{%
R_{E}}}}{{R}}\right) ^{8}\left( h_{3}^{1}g_{1}^{1}-g_{3}^{1}h_{1}^{1}\right).
\end{equation*}

In accordance with such expansion the matrix $\mathbf{M}$ will have the form
\begin{equation}\label{eq138}
\mathbf{M}=\mathbf{M}^{(0)}+\sin
i \,\mathbf{M}^{(1)}(t)=\mathbf{M}_{cp}^{(0)}+\widetilde{\mathbf{M}}^{(0)}(t)+\sin
i \,\mathbf{M}^{(1)}(t),
\end{equation}
 where the matrix $\mathbf{M}_{cp}^{(0)}$ will look like
{\small {\[\mathbf{M}_{cp}^{(0)}=\left(
  \begin{array}{ccc}
    4\omega_0^2(B-C)-l_{11\, cp}^{(0)}-b_{11\, cp}^{(0)}\!\!\!\!\!\!\!\!& -b_{12\, cp}^{(0)}\!\!\!\!\!\!\!\!& -b_{13\, cp}^{(0)} \\
     -b_{21\, cp}^{(0)}\!\!\!\!\!\!\!\!&3\omega_0^2(A-C)-l_{22\, cp}^{(0)}-b_{22\, cp}^{(0)}\!\!\!\!\!\!\!\!&-l_{23\, cp}^{(0)}-b_{23\, cp}^{(0)} \\
      -b_{31\, cp}^{(0)}\!\!\!\!\!\!\!\!&-l_{32\, cp}^{(0)}-b_{32\, cp}^{(0)}\!\!\!\!\!\!\!\!&\omega_0^2(B-A)-l_{33\, cp}^{(0)}-b_{33\, cp}^{(0)}
  \end{array}\right).
\]}}

At the inequalities
\begin{equation}\label{eq139}
\left\{
\begin{array}{l}
4\omega_0^2(B-C)-l_{11\,cp}^{(0)}-b_{11\,cp}^{(0)}>0\\
(4\omega_0^2(B-C)-l_{11\,cp}^{(0)}-b_{11\,cp}^{(0)})(3\omega_0^2(A-C)-l_{22\,cp}^{(0)}-b_{22\,cp}^{(0)})-({b_{12\,cp}^{(0)}})^2>0\\
\det \mathbf{M}_{cp}^{(0)}>0
\end{array}\right.
\end{equation}
hold true the matrix $\mathbf{M}_{cp}^{(0)}$\ is positive--definite and hence the zero solution of the system
\[\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}^{(0)}_{cp}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)=0\] is asymptotically stable.

In assumption that control moments $\vec{M}_{L}$ and $\vec{M}_{M}$ are absent, the spacecraft is under influence of gravitational moment $\vec{M}_{G}$ only. In this case the inequalities (\ref{eq139}) degenerates into the well-known
Beletsky conditions \cite{Beletsky1975} of spacecraft gravitational stabilization: $B>A>C$.
The presence of control moments $\vec{M}_{L}$ and $\vec{M}_{M}$, as is shown with the help of computer modeling, lead to essentially increasing of the region corresponding to Beletsky conditions. For example the region of validity of inequalities
(\ref{eq139}), constructed on the plane of dimensionless inertial parameters $\delta=B/A$ and $\varepsilon=C/A$, is shown by horizontal shading in the Fig. \ref{Fig40} for the following parameters of spacecraft and its orbit:
\begin{eqnarray}
R &=&7\cdot 10^{6}\;\mbox{m,}\ i=0\,\mbox{rad,}\ A=10^{3}\;\mbox{kg}\cdot %
\mbox{m}^{2},\ Q=5\cdot 10^{-3}\;\mbox{C,}\ \   \notag \\
k_{L} &=&10\;\mbox{s}/\mbox{T},\ k_{M}=2\cdot 10^{6}\;\mbox{N}\cdot \mbox{m}/%
\mbox{T}^{2}  \label{eq140}
\end{eqnarray}

This region, evidently from the Fig. \ref{Fig40}, is rather wide. This is the evidence of efficiency of proposed method for such parameters $\delta$ and $\varepsilon$, which are beyond the bounds of Beletsky triangle ($\delta>1,\ \varepsilon<1$) \cite{Beletsky1975}, shown by vertical shading in the same figure.
\begin{figure}[!h]
\centering
\includegraphics[width=10cm]{f40.eps}
\caption{} \label{Fig40}
\end{figure}

This implies that this method is applicable to the spacecraft that cannot be stabilized by the moment of gravitational forces.
At the same time parameter $k_{M}$, taken in this example as equal to $2\cdot 10^6\; \mbox{N}\cdot\mbox{m}/\mbox{T}^2$, ensures the arising of sufficiently small control moment $\vec{M}_{M}$, not even exceeding the gravitational moment $\vec{M}_{G}$ in degree of value.

Comparing the parameters regions obtained in the case of joint use of the moments $\vec M_{M}$ and $\vec M_{L}$ (Fig. \ref{Fig40}) with the corresponding regions constructed for the cases when only one of the moments $\vec M_{M}$ or $\vec M_{L}$ is in use, we notice the significantly increasing of the mentioned region in conditions of joint use of the moments.
The results of computer modeling for parameters taken from (\ref{eq140}) except the cases when one of the moments $\vec M_{M}$ or $\vec M_{L}$ is absent and correspondingly $k_{M}=0$ or $ k_{L}=0$, are shown in Fig. \ref{Fig41}.
\begin{figure}[!h]
\centering
\includegraphics[width=10cm]{f41.eps}
\caption{} \label{Fig41}
\end{figure}

The horizontal shading in this figure mark the region corresponding to the use of only one control moment $\vec M_{L}$,
the vertical shading mark the region corresponding to the use of only one control moment $\vec M_{M}$, the grey color mark the region which arises only thanks to the joint use of the both control moments. The union of mentioned three regions give the resultant region of asymptotic stability, corresponding to the inequalities (\ref{eq139}).

Further we consider the matrix $\widetilde{\bold M}^{(0)}(t)= \bigl(\widetilde m_{ij}^{(0)}(t)\bigr)$ and estimate its norm.
\begin{multline*}
\Vert \widetilde{\mathbf{M}}^{(0)}(t)\Vert =\sum\limits_{i,j=1}^{3}|%
\widetilde{m}_{ij}^{(0)}(t)|=|m_{11}^{(0)}(t)-m_{11%
\,cp}^{(0)}|+2|m_{12}^{(0)}(t)-m_{12\,cp}^{(0)}|+2|m_{13}^{(0)}(t) \\
-m_{13\,cp}^{(0)}|+|m_{22}^{(0)}(t)-m_{22%
\,cp}^{(0)}|+2|m_{23}^{(0)}(t)-m_{23\,cp}^{(0)}|+|m_{33}^{(0)}(t)-m_{33%
\,cp}^{(0)}| \\
=|l_{11\,cp}^{(0)}+b_{11\,cp}^{(0)}-l_{11}^{(0)}(t)-b_{11}^{(0)}(t)|+2|b_{12%
\,cp}^{(0)}-b_{12}^{(0)}(t)|+2|b_{13\,cp}^{(0)}-b_{13}^{(0)}(t)| \\
+|l_{22\,cp}^{(0)}+b_{22\,cp}^{(0)}-l_{22}^{(0)}(t)-b_{22}^{(0)}(t)|+2|l_{23%
\,cp}^{(0)}+b_{23\,cp}^{(0)}-l_{23}^{(0)}(t)-b_{23}^{(0)}(t)|+|l_{33%
\,cp}^{(0)} \\
+b_{33\,cp}^{(0)}-l_{33}^{(0)}(t)-b_{33}^{(0)}(t)|==|(Qk_{L}R^{2}(\omega
_{0}-\omega _{E})^{2}+k_{M})(B_{\eta }^{2}+B_{\zeta }^{2}-{\langle B_{\eta
}^{2}\rangle }_{t} \\
-{\langle B_{\zeta }^{2}\rangle }_{t})|+2|k_{M}({\langle B_{\xi }B_{\eta
}\rangle }_{t}-B_{\xi }B_{\eta })|+2|k_{M}({\langle B_{\xi }B_{\zeta
}\rangle }_{t}-B_{\xi }B_{\zeta })| \\
+|Qk_{L}R^{2}(\omega _{0}-\omega _{E})^{2}(B_{\eta }^{2}-{\langle B_{\eta
}^{2}\rangle }_{t})+k_{M}(B_{\xi }^{2}+B_{\zeta }^{2}-{\langle B_{\xi
}^{2}\rangle }_{t}-{\langle B_{\zeta }^{2}\rangle }_{t})| \\
+2|(Qk_{L}R^{2}(\omega _{0}-\omega _{E})^{2}-k_{M})(B_{\eta }B_{\zeta }-{%
\langle B_{\eta }B_{\zeta }\rangle }_{t})| \\
+|Qk_{L}R^{2}(\omega _{0}-\omega _{E})^{2}(B_{\zeta }^{2}-{\langle B_{\zeta
}^{2}\rangle }_{t})+k_{M}(B_{\xi }^{2}+B_{\eta }^{2}-{\langle B_{\xi
}^{2}\rangle }_{t}-{\langle B_{\eta }^{2}\rangle }_{t})| \\
<6\left( R^{2}(\omega _{0}-\omega _{E})^{2}|Qk_{L}|+2k_{M}\right) \left(
B_{\eta }^{2}+{\langle B_{\eta }^{2}\rangle }_{t}\right)  \\
<108\left( (\omega _{0}-\omega _{E})^{2}|Qk_{L}|+2k_{M}\right) {g_{1}^{0}}%
^{2}\frac{R_{E}^{6}}{R^{4}}\equiv C_{1}.
\end{multline*}

Since $C_1=\text{const}$ depends only the matrix $\bold M_{cp}^{(0)}$, then the zero solution of differential system
(\ref{eq97}) is also asymptotically stable {\cite {BR 1954}}. At that this asymptotical stability is uniform since the coefficients in the system (\ref{eq97}) are almost-periodic in $t$. Furthermore the exponential asymptotic stability of the zero solution of the system (\ref{eq97}) implies the exponential asymptotic stability of the zero solution of the initial nonlinear system (\ref{eq7}) at $i=0$ {\cite {Malkin1966}}. This prove the possibility of spacecraft stabilization at
$i=0$ with the help of the moments $\vec{M}_{L}$ and $\vec{M}_{M}$.

In more general case i.e. in the case of slightly inclined orbits ($i\neq 0$ but $\sin i$ is small) it is appropriate to write the differential equations of perturbed motion (\ref{eq7}) in the form
\begin{equation}\label{eq141}
\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}^{(0)}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)+\sin i \cdot \mathbf{M}^{(1)}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)+\vec X=0.
\end{equation}

According to the theorem of total stability \cite {Malkin1966} the uniform asymptotic stability of the zero solution of the system (\ref{eq97}) is the sufficient condition of the total stability (the stability in the presence of constantly affecting perturbations) of this solution. As the constantly affecting perturbations one can consider small in their norm perturbations
\begin{equation}\label{eq142}
f=(f_1,f_2,f_3)^{\top}=\sin i\,\bold M^{(1)}(\varphi,\theta,\psi)^{\top} +
(X_1, X_2, X_3)^{\top}.
\end{equation}

Furthermore, the exponential asymptotic stability of the zero solution of the system (\ref{eq97}) according to the inequalities
\begin{equation}\label{eq143}
|f_j(t,\varphi,\psi,\theta)|<\sin i\,
(|\varphi|+|\theta|+|\psi|) \ \ (j=\overline{1,3})
\end{equation}
implies the asymptotic stability of the zero solution of the initial nonlinear system (\ref{eq7}) for sufficiently small values of $i$ {\cite {Malkin1966}}.

Thus, it was proved the total stability of direct equilibrium position and asymptotic stability at sufficiently small values of $i$ and therefore the adaptability of suggested stabilization method based on the joint use of the moments $\vec M_{M}$ and $\vec M_{L}$ is well founded for the orbits of small inclinations.

The computer modeling was implemented in order to confirm the efficiency of suggested stabilization method and approbate it for concrete values of parameters of spacecraft, its orbit and initial conditions of motion.
The results shown in Fig.\ref{Fig42}, \ref{Fig43}, \ref{Fig44} were obtained for the following parameters values:
\noindent $R=7\cdot10^6\; \mbox{m,}\ i=0\; \mbox{rad,}\
A=10^3\; \mbox{kg}\cdot\mbox{m}^2,\ Q=5\cdot10^{-3}\; \mbox{C,}\
\delta=2.2,\ \varepsilon=1.4,\ h_1=h_2=h_3=h=0.5A\omega_0\; \mbox{kg}\cdot\mbox{m}^2/\mbox{s} , \ k_{L}=10\; \mbox{s}/\mbox{T}, \ k_{M}=2\cdot 10^6\; \mbox{N}\cdot\mbox{m}/\mbox{T}^2 $.
The initial conditions of motion were chosen as follows $\varphi (0)=0.2\, \mbox{rad,}\ \psi
(0)=-0.2\, \mbox{rad,}\ \theta (0)=0.2\, \mbox{rad,}\
\omega_x(0)=0.1\omega_0\, ,$ $\ \ \omega_y(0)=1.1\omega_0\, , \ \
\omega_z(0)=0.1\omega_0\, .$
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f42.eps}
\caption{} \label{Fig42}
\end{figure}
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f43.eps}
\caption{} \label{Fig43}
\end{figure}
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f29.eps}
\caption{} \label{Fig29}
\end{figure}
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f44.eps}
\caption{} \label{Fig44}
\end{figure}

In Fig.~\ref{Fig42} the graphs of "airborne" angles via the argument of a latitude are plotted for spacecraft under the action of control moments $\vec M_{L}$ and $\vec M_{M}$.
In Fig.~\ref{Fig43} the graphs of quaternion components are plotted for the case of stabilized motion, and in Fig.~\ref{Fig29} are the same components in the absence of control moments $\vec M_{L}$ ($k_{L}=0$), $\vec M_{M}$ ($k_{M}=0$) and the damping moment $\vec M_{D}$ ($h=0$). The values of the moments $\vec{M}_{L}, \ \vec{M}_{G},\ \vec{M}_{D}$  and $ \vec{M}_{M}$, acting upon the spacecraft in the process of stabilization, are shown in Fig.~\ref{Fig44}. It is obvious that the mentioned values are in such limits that make evident the operability of suggested method of attitude control. The comparison of obtained results with similar results obtained with the use of only one of the moments $\vec M_{L}$ or $\vec M_{M}$ is made. It is revealed that the joint use of both control moments $\vec {M}_{L}$ and $\vec M_{M}$ in electrodynamical attitude stability system result in reduce of the transition time.

%\begin{figure}
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f43}
%\hfill
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f44}
%\\
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig43}}
%\hfill
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig44}}
%\end{figure}

\section{Analysis of stabilization process at the orbits of mean and large inclinations}

For substantiation of suggested attitude control method let us again turn to the dynamical Euler equations (\ref{eq7}) and at first consider the system of linear approximation of these equations:
\begin{equation}\label{eq99}
\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)=0.
\end{equation}

Represent the matrix $\mathbf{M}$ in the form $\mathbf{M}=\mathbf{M}_{cp}+\widetilde{\mathbf{M}}(t).$ Here
the matrix $\mathbf{M}_{cp}$ is the result of time componentwise averaging of the matrix $\mathbf{M}$, i.e.,
%\begin{equation}
{\small{\[\mathbf{M}_{cp}=\left(
  \begin{array}{ccc}
    4\omega_0^2{{(B-C)}}-l_{11cp} -b_{11cp} \!\!\!\!\!\!\!\!\!\!&-l_{12cp} -b_{12cp} \!\!\!\!\!\!\!\!\!\!&-l_{13cp} - b_{13cp}  \\
    -l_{21cp} - b_{21cp} \!\!\!\!\!\!\!\!\!\! & 3\omega_0^2{{(A-C)}}-l_{22cp} - b_{22cp} \!\!\!\!\!\!\!\!\!\!& -l_{23cp} - b_{23cp} \\
    -l_{31cp} - b_{31cp} \!\!\!\!\!\!\!\!\!\!&-l_{32cp} - b_{32cp} \!\!\!\!\!\!\!\!\!\!&\omega_0^2(B-A)-l_{33cp} - b_{33cp}
  \end{array}\!\right) ,
\]}}
%\end{equation}
where
\begin{equation}\label{eq145}
l_{11cp}=-Qk_{L}[{\langle v_{C\eta}^2B_{\xi}^2\rangle}_t-2 v_{C\xi}
{\langle v_{C\eta}B_{\xi}B_{\eta}\rangle }_t+v_{C\xi}^2({\langle
B_{\eta}^2\rangle}_t+ {\langle B_{\zeta}^2\rangle }_t)],
\end{equation}
\begin{equation}\label{eq146}
l_{12cp}=l_{21cp}=-Qk_{L}v_{C\xi}{\langle
v_{C\eta}B_{\zeta}^2\rangle }_t,
\end{equation}
\begin{equation}\label{eq147}
l_{13cp}=l_{31cp}=Qk_{L}(v_{C\xi}{\langle
v_{C\eta}B_{\eta}B_{\zeta}\rangle }_t -{\langle
v_{C\eta}^2B_{\xi}B_{\zeta}\rangle }_t),
\end{equation}
\begin{equation}\label{eq148}
l_{22cp}=-Qk_{L}[{\langle v_{C\eta}^2B_{\xi}^2\rangle}_t-2 v_{C\xi}
{\langle v_{C\eta}B_{\xi}B_{\eta}\rangle }_t+v_{C\xi}^2{\langle
B_{\eta}^2\rangle}_t+ {\langle v_{C\eta}^2B_{\zeta}^2\rangle
}_t],
\end{equation}
\begin{equation}\label{eq149}
l_{23cp}=l_{32cp}=-Qk_{L}v_{C\xi}(v_{C\xi}{\langle
B_{\eta}B_{\zeta}\rangle }_t -{\langle
v_{C\eta}B_{\xi}B_{\zeta}\rangle }_t),
\end{equation}
\begin{equation}\label{eq150}
l_{33cp}=-Qk_{L}(v_{C\xi}^2{ \langle B_{\zeta}^2\rangle }_t+{ \langle v_{C\eta}^2B_{\zeta}^2\rangle }_t) .
\end{equation}
\begin{equation*}%\label{}
\begin{array}{ll}
b_{11cp}=-k_{M}({\langle B_{\eta}^2\rangle}_t+{\langle B_{\zeta}^2\rangle}_t),&
b_{12cp}=b_{21cp}=k_{M}{\langle B_{\xi}B_{\eta}\rangle }_t,\\ b_{13cp}=b_{31cp}=k_{M}{\langle
B_{\xi}B_{\zeta}\rangle }_t ,&
b_{22cp}=-k_{M}({\langle B_{\xi}^2\rangle}_t+{\langle B_{\zeta}^2\rangle}_t),\\
b_{23cp}=b_{32cp}=k_{M}{\langle B_{\eta}B_{\zeta}\rangle }_t,& b_{33cp}=-k_{M}({\langle
B_{\xi}^2\rangle}_t+{\langle B_{\eta}^2\rangle}_t) \, ,
\end{array}
\end{equation*}

The averaged values appearing in $l_{ijcp}$ and $b_{ijcp}$ are calculated with the use of components of vector $\vec B$
in octupole approximation (\ref{eq151}) -- (\ref{eq153}):
\begin{multline}\label{eq112}
{\langle
B_{\eta}^2\rangle}_t=
\dfrac{1}{2}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{6}
\left( \left( 2\,{g_1^0}^{2}-{g_1^1}^{2}-{h_1^1}^{2} \right) \cos^{2}(i) +{h_1^1}^{2}+{g_1^1}^{2}
\right) \\
+\dfrac{3}{8}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{8}
\left(\left( -12\,{g_2^0}^{2}+8\,{g_2^1}^{2}+5\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -2\,{g_2^2}^{2}+8\,{h_2^1}^{2}
\right.\right.
\\
\left.\left.
-20\,g_1^0g_3^0-2\,{h_2^2}^{2} \right)\cos^{4} (i)+\left( 12 \left({g_2^0}^{2}+g_1^0 g_3^0\right)-6\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right)
\right.\right.
\\
\left.\left.
-6\,\left({g_2^1}^{2}+{h_2^1}^{2}\right) \right) \cos^{2}( i )+2({g_2^1}^{2}+{g_2^2}^{2}+{h_2^1}^{2}+{h_2^2}^{2})+\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right)
\right) \\
+\dfrac{3}{128}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{10}
\left(
 45\left( 20\,{g_3^0}^{2}-\left(\,{g_3^3}^{2}+\,{h_3^3}^{2}\right)-15\left({h_3^1}^{2}+{g_3^1}^{2}\right)
 \right.\right.
 \\
 \left.\left.
 +6\left({h_3^2}^{2}+{g_3^2}^{2}\right) \right)  \cos^{6} ( i )  + 5\left( 3\left({h_3^3}^{2}+{g_3^3}^{2}\right)+181\left({g_3^1}^{2}+{h_3^1}^{2}\right)
 \right.\right.
 \\
 \left.\left.
 -52\left({h_3^2}^{2}+{g_3^2}^{2}\right)-264\,{g_3^0}^{2} \right)\cos^{4} (i)+ 3\left( 172{g_3^0}^{2}+10\left({g_3^2}^{2}+{h_3^2}^{2}\right)
 \right.\right.
 \\
 \left.\left.
 -91\left({h_3^1}^{2}+{g_3^1}^{2}\right)-5\left({g_3^3}^{2}+{h_3^3}^{2}\right) \right) \cos^{2} (i)
\right.
 \\
\left.
 +40\left({h_3^2}^{2}+{g_3^2}^{2}\right)+45\left({g_3^3}^{2}+{h_3^3}^{2}\right)+43\left({h_3^1}^{2}+\,{g_3^1}^{2}\right)\right)
\end{multline}
\begin{multline}\label{eq113}
{\langle
v_{C\eta}^2B_{\xi}^2\rangle}_t=
-\dfrac{\omega^2_{E}\sin^2 i}{16}
\dfrac{R_{E}^6}{ R^4}
\left(
\left( 6\,{g_1^0}^{2}-3\,{g_1^1}^{2}-3\,{h_1^1}^{2}
\right)  \cos^{2} ( i ) -{g_1^1}^{2}-{
h_1^1}^{2}
\right.
\\
\left.
-6\,{g_1^0}^{2}\right)-
\dfrac{\omega^2_{E}\sin^2 i}{64}
\dfrac{R_{E}^8}{ R^6}
\left(
3\, \left( 8\left({h_2^1}^{2}+{g_2^1}^{2}\right)-2\left({g_2^2}^{2}+
\,{h_2^2}^{2}\right)
\right.\right.
\\
\left.\left.
+5\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -12\,{g_2^0}^{2}-20\,g_1^0g_3^0 \right)
 \cos^{4} ( i ) +2\, \left( 36\,{g_2^0}
^{2}+24\,g_1^0g_3^0
\right.\right.
\\
\left.\left.
-7\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -18\left({g_2^2}^{2}+\,{h_2^2}^{2}\right)
 \right)  \cos^{2}( i ) +7\,\sqrt {6}
 \left( h_3^1h_1^1+g_3^1g_1^1 \right)
\right.
\\
\left.
-24\left({h_2^1}
^{2}+\,{g_2^1}^{2}\right)-6\left({g_2^2}^{2}+\,{h_2^2}^{2}\right)+12\,g_1^0g_3^0-36\,{g_2^0}^{2}
\right)
\\
-\dfrac{3\omega^2_{E}\sin^2 i}{2048}
\dfrac{R_{E}^{10}}{ R^8}
\left(
45\, \left( 6\left({h_3^2}^{2}+\,{g_3^2}^{2}\right)-15\left({h_3^1}^{2}+\,{g_3^1}^{2}\right)-{h_3^3}^{2}-{g_3^3}^{2}
\right.\right.
\\
\left.\left.
+20\,{g_3^0}^{2} \right)  \cos^{6} ( i )  +5\, \left( 62\left({h_3^2}^{2}+\,{g_3^2}^{2}\right)+157\left({h_3^1}^{2}+\,{g_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-45\left({h_3^3}^{2}+{g_3^3}^{2}\right)-348\,{g_3^0}^{2} \right)   \cos^{4} ( i )+ \left( 247\left({h_3^1}^{2}+{g_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-615\left({h_3^3}^{2}+{g_3^3}^{2}\right)-310\left({g_3^2}^{2}+{h_3^2}^{2}\right)+1356\,{g_3^0}^{2} \right)   \cos^{2} ( i )
\right.
\\
\left.
-270\left({g_3^2}^{2}+\,{h_3^2}^{2}\right)-75\left({h_3^3}^{2}+\,{g_3^3}^{2}\right)-421\left({h_3^1}^{2}+{g_3^1}^{2}\right)-516\,{g_3^0}^{2}
\right)
\end{multline}
\begin{multline}\label{eq114}
{\langle
v_{C\eta}B_{\xi}B_{\eta}\rangle}_t=
\dfrac{\omega_{E}\sin^2 i\cos{i}}{4}
\dfrac{R_{E}^6}{ R^5}
\left(2\,{g_1^0}^{2}-{h_1^1}^{2}-{g_1^1}^{2}
\right)
\\
+\dfrac{\omega_{E}\sin^2 i\cos{i}}{32}
\dfrac{R_{E}^8}{ R^7}
\left(
3\, \left( 8\left({g_2^1}^{2}+{h_2^1}^{2}\right)-2\left({g_2^2}^{2}+{h_2^2}^{2}\right)
\right.\right.
\\
\left.\left.
+5\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -20\,g_1^0g_3^0-12\,{g_2^0}^{2} \right)  \cos^{3} ( i )   + \left( 36\,{g_2^0}^{2}+12\,g_1^0g_3^0
\right.\right.
\\
\left.\left.
-18\left({g_2^2}^{2}+{h_2^2}^{2}\right)-11\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) \right) \cos ( i )
\right)
\\
+\dfrac{3\omega_{E}\sin^2 i\cos{i}}{256}
\dfrac{R_{E}^{10}}{ R^9}
\left(
15\, \left( 6\left({h_3^2}^{2}+{g_3^2}^{2}\right)-15\left({h_3^1}^{2}+{g_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-{h_3^3}^{2}-{g_3^3}^{2}+20\,{g_3^0}^{2} \right)   \cos^{5} ( i )   +10\, \left( 2\left({h_3^2}^{2}+{g_3^2}^{2}\right)-3\left({h_3^3}^{2}+{g_3^3}^{2}\right)
\right.\right.
\\
\left.\left.
+19\left({h_3^1}^{2}+{g_3^1}^{2}\right)-36\,{g_3^0}^{2} \right)  \cos^{3} ( i )   +
3\, \left( 9\left({g_3^1}^{2}+{h_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-10\left({h_3^2}^{2}+\,{g_3^2}^{2}\right)-25\left({h_3^3}^{2}+{g_3^3}^{2}\right)+52\,{g_3^0}^{2} \right) \cos \left( i \right)
\right)
\end{multline}
\begin{multline}\label{eq115}
{\langle
B^2_{\zeta}\rangle}_t=
\left(\dfrac{R_{E}}{ R}\right)^6
\left(
 \left( {g_1^1}^{2}+{h_1^1}^{2}-2\,{g_1^0}^{2} \right)  \cos^{2} ( i )   +{h_1^1}^{2}+2\,{g_1^0}^{2}+{g_1^1}^{2}\right)
 \\
 -\dfrac{1}{64}\left(\dfrac{R_{E}}{ R}\right)^8
\left(
 3\, \left( 108\left({h_2^1}^{2}+{g_2^1}^{2}\right)-27\left({g_2^2}^{2}+\,{h_2^2}^{2}\right)
\right.\right.
\\
\left.\left.
+80\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -162\,{g_2^0}^{2}-320\,g_3^0g_1^0 \right)   \cos^{4} ( i )   +3\, \left( 180\,{g_2^0}^{2}
\right.\right.
\\
\left.\left.
+384\,g_3^0g_1^0-18\left({g_2^2}^{2}+\,{h_2^2}^{2}\right)-72\left({h_2^1}^{2}+{g_2^1}^{2}\right)
\right.\right.
\\
\left.\left.
-32\,\sqrt {6} \left(h_3^1h_1^1+g_3^1g_1^1 \right)  \right)   \cos^{2} ( i )  -81\left({g_2^2}^{2}+{h_2^2}^{2}\right)-16\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right)
\right.
\\
\left.
 -108\left({h_2^1}^{2}+{g_2^1}^{2}\right)-192\,g_3^0g_1^0-198\,{g_2^0}^{2}
\right)
-\dfrac{1}{16}\left(\dfrac{R_{E}}{ R}\right)^{10}
\left(
25 \left( 6\left({h_3^2}^{2}+{g_3^2}^{2}\right)
\right.\right.
\\
\left.\left.
-15\left({h_3^1}^{2}+{g_3^1}^{2}\right)-{h_3^3}^{2}-{g_3^3}^{2}+20\,{g_3^0}^{2} \right)   \cos^{6} ( i )   +15\, \left( 33\left({g_3^1}^{2}+{h_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-6\left({h_3^2}^{2}+{g_3^2}^{2}\right)-{g_3^3}^{2}-{h_3^3}^{2}-52\,{g_3^0}^{2} \right)  \cos^{4} ( i )   +3 \left( 116\,{g_3^0}^{2}-5\left({g_3^3}^{2}+\,{h_3^3}^{2}\right)
\right.\right.
\\
\left.\left.
-43\left({g_3^1}^{2}+{h_3^1}^{2}\right)-10\left({h_3^2}^{2}+{g_3^2}^{2}\right) \right)  \cos^{2} ( i )
\right.
\\
\left.
 -25\left({g_3^3}^{2}+{h_3^3}^{2}\right)-39\left({h_3^1}^{2}+{g_3^1}^{2}\right)-30\left({h_3^2}^{2}+{g_3^2}^{2}\right)
 -68\,{g_3^0}^{2}
\right)
\end{multline}
\begin{multline}\label{eq116}
{\langle
v^2_{C\eta}B^2_{\zeta}\rangle}_t=
-\dfrac{\omega^2_{E}\sin^2 i}{4}
\dfrac{R_{E}^6}{ R^4}
\left(
\left( 2\,{g_1^0}^{2}-{g_1^1}^{2}-{h_1^1}^{2} \right)
 \cos^{2}( i )
\right.
\\
\left.
 -3\,{h_1^1}^{2}-3\,{g_1^1}^{2}-2\,{g_1^0}^{2}
\right)-\dfrac{\omega^2_{E}\sin^2 i}{128}
\dfrac{R_{E}^8}{ R^6}
\left(
 \left( 108\left({g_2^1}^{2}+{h_2^1}^{2}\right)
\right.\right.
\\
\left.\left.
 +80\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -27\left({g_2^2}^{2}+{h_2^2}^{2}\right)-162\,{g_2^0}^{2}-320\,g_3^0g_1^0 \right)  \cos^4 (i)
\right.
\\
\left.
+2\, \left( 16\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right) -27\left({g_2^2}^{2}+{h_2^2}^{2}\right)+128\,g_3^0g_1^0+54\,{g_2^0}^{2} \right)   \cos^2 (i)
\right.
\\
\left.
-135\left({g_2^2}^{2}+{h_2^2}^{2}\right)-108\left({g_2^1}^{2}+{h_2^1}^{2}\right)+16\,\sqrt {6} \left( h_3^1h_1^1+g_3^1g_1^1 \right)
\right.
\\
\left.
+64\,g_3^0g_1^0-90\,{g_2^0}^{2}
 \right)-\dfrac{\omega^2_{E}\sin^2 i}{128}
\dfrac{R_{E}^{10}}{ R^8}
\left(
25\, \left( 6\left({h_3^2}^{2}+{g_3^2}^{2}\right)-15\left({h_3^1}^{2}+{g_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-{h_3^3}^{2}-{g_3^3}^{2}+20\,{g_3^0}^{2} \right)  \cos^6 (i)+15\, \left( 19\left({h_3^1}^{2}+{g_3^1}^{2}\right)-3\left({g_3^3}^{2}+{h_3^3}^{2}\right)
\right.\right.
\\
\left.\left.
+2\left({h_3^2}^{2}+{g_3^2}^{2}\right)-36\,{g_3^0}^{2} \right) \cos^4 (i)+3\, \left( 9\left({h_3^1}^{2}+{g_3^1}^{2}\right)-25\left({g_3^3}^{2}+{h_3^3}^{2}\right)
\right.\right.
\\
\left.\left.
-10\left({h_3^2}^{2}+{g_3^2}^{2}\right)+52\,{g_3^0}^{2} \right)  \cos^2 (i)-175\left({g_3^3}^{2}+{h_3^3}^{2}\right)
\right.
\\
\left.
-129\left({h_3^1}^{2}+{g_3^1}^{2}\right)-150\left({h_3^2}^{2}+{g_3^2}^{2}\right)-116\,{g_3^0}^{2}
\right)
\end{multline}
\begin{multline}\label{eq117}
{\langle
v_{C\eta}B_{\xi}B_{\zeta}\rangle}_t=
-\dfrac{\omega_{E}\sin^2 i}{16}
\dfrac{R_{E}^7}{ R^6}
\left(
7\, \left( \sqrt {3} \left( h_2^1h_1^1+g_2^1g_1^1 \right) -3\,g_2^0g_1^0 \right)  \cos^2 (i)
\right.
\\
\left.
+\sqrt {3} \left( h_2^1h_1^1+g_2^1g_1^1 \right) +9\,g_2^0g_1^0
\right)-\dfrac{3\omega_{E}\sin^2 i}{128}
\dfrac{R_{E}^9}{ R^8}
\left(
\left( 17\,\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right)
\right.\right.
\\
\left.\left.
-85\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right) +170\,g_3^0g_2^0 \right)  \cos^4 (i)+
 \left( 48\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right)
 \right.\right.
\\
\left.\left.
+6\,\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right) -148\,g_3^0g_2^0 \right)  \cos^2 (i)
\right.
\\
\left.
+\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right) +5\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right) +26\,g_3^0g_2^0
\right)
\end{multline}
\begin{multline}\label{eq118}
{\langle
B_{\eta}B_{\zeta}\rangle}_t=
-\dfrac{\cos{i}}{4}
\left(\dfrac{R_{E}}{ R}\right)^7
\left(7\, \left( \sqrt {3} \left( h_2^1h_1^1+g_2^1g_1^1 \right) -3\,g_2^0g_1^0 \right)  \cos^2 (i)
\right.
\\
\left.
+15\,g_2^0g_1^0-3\,\sqrt {3} \left( h_2^1h_1^1+g_2^1g_1^1 \right)
\right)-\dfrac{3\cos{i}}{64}
\left(\dfrac{R_{E}}{ R}\right)^9
\left(
 \left( 51\,\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right)
\right.\right.
\\
\left.\left.
-255\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right) +510\,g_3^0g_2^0 \right)  \cos^4 (i)+ \left( 266\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right)
\right.\right.
\\
\left.\left.
-22\,\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right) -636\,g_3^0g_2^0 \right)  \cos^2 (i)
\right.
\\
\left.
-5\,\sqrt {5} \left( g_3^2g_2^2+h_3^2h_2^2 \right) -43\,\sqrt {2} \left( h_2^1h_3^1+g_2^1g_3^1 \right) +174\,g_3^0g_2^0
 \right)
\end{multline}
\begin{equation}\label{eq119}
{\langle
v_{C\eta}B^2_{\zeta}\rangle}_t=0
\end{equation}
\begin{equation}\label{eq120}
{\langle
v^2_{C\eta}B_{\xi}B_{\zeta}\rangle}_t=
-\dfrac{\sqrt{6}\omega^2_{E}\sin^2 i \cos i}{32}
\dfrac{R^8_{E}}{ R^6}
\left( h_3^1g_1^1-h_1^1g_3^1 \right)  \left( 5\,
\cos^2 i -1 \right)
\end{equation}
\begin{equation}\label{eq121}
{\langle
v_{C\eta}B_{\eta}B_{\zeta}\rangle}_t=
\dfrac{\sqrt{6}\omega_{E}\sin^2 i }{32}
\dfrac{R^8_{E}}{R^7}
\left( h_3^1g_1^1-h_1^1g_3^1 \right)  \left( 5\,
 \cos^2 i -1 \right)
\end{equation}
\begin{multline*}%\label{}
{\langle
B_{\xi}^2\rangle}_t=
\dfrac{1}{4}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{6}
\left( \left( {h_1^1}^{2}-2\, {g_1^0}^2+{g_1^1}^{2} \right) \cos^{2}(i) +{g_1^1}^{2}+2\, {g_1^0}^2+{h_1^1}^{2}
\right) \\
+\dfrac{1}{32}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{8}
\left( 3
\left(
20\,g_1^0 g_3^0-5\,\sqrt {6}\left(h_3^1h_1^1+g_1^1g_3^1\right)+2\,\left({h_2^2}^{2}+{g_2^2}^{2}\right)+12\,{g_2^0}^{2}
\right.\right.
\\
\left.\left.
-8\,\left({g_2^1}^{2}+\,{h_2^1}^{2}\right)
\right)\cos^{4} (i)\right)
+\left( 6 \left(
6\,\left({g_2^2}^{2}+{h_2^2}^{2}\right)-12\,\left(g_3^0+\,{g_2^0}^{2}\right)
\right.\right.
\\
\left.\left.
+\sqrt {6}\left(g_1^1g_3^1+h_3^1h_1^1\right)
\right)\cos^{2}( i )+
\sqrt {6}\left(g_1^1g_3^1+h_3^1h_1^1\right)+6\left(2\,g_1^0g_3^0+3\,{g_2^0}^{2}
\right.\right.
\\
\left.\left.
+4\,{h_2^1}^{2}+4\,{g_2^1}^{2}+
{g_2^2}^{2}+{h_2^2}^{2}\right)
\right)
+\dfrac{3}{256}
\left(
\dfrac {{{R_{E}}}}{{R}}\right)^{10}
\left(15\left(15\,\left({h_3^1}^{2}+{g_3^1}^{2}\right)
\right.\right.
\\
\left.\left.
-6\,\left({h_3^2}^{2}+{g_3^2}^{2}\right)-20\,{g_3^0}^{2}
+{g_3^3}^{2}-{h_3^3}^{2}
\right)\cos^{6} ( i )+
5\left(21\,\left({g_3^3}^{2}+{h_3^3}^{2}\right)
\right.\right.
\\
\left.\left.
-53\,\left({g_3^1}^{2}+\,{h_3^1}^{2}\right)-34\,\left({h_3^2}^{2}+{g_3^2}^{2}\right)+132\,{g_3^0}^{2}
\right)\cos^{4} (i)
\right.
\\
\left.
+3\left(35\left(2\,{g_3^2}^{2}+2\,{h_3^2}^{2}+{g_3^3}^{2}+{h_3^3}^{2}\right)-19\left({h_3^1}^{2}
+\,{g_3^1}^{2}\right)-172\,{g_3^0}^{2}
\right)\cos^{2} (i)
 \right.
 \\
 \left.
+50\,\left({g_3^2}^{2}+{h_3^2}^{2}\right)+15\,\left({g_3^3}^{2}+{h_3^3}^{2}\right)+
113\,\left({h_3^1}^{2}+{g_3^1}^{2}\right)+156\,{g_3^0}^{2}
\right),
\end{multline*}
\begin{multline*}%\label{}
{\langle
B_{\xi}B_{\eta}\rangle}_t=
\dfrac{\sqrt{3}}{4}
\left(\dfrac{R_{E}}{ R}\right)^7\left(3\,\left(h_1^1g_2^1-\,h_2^1g_1^1\right)\cos^2 (i)
+\left(h_2^1g_1^1-h_1^1g_2^1\right)\right)
\\
+\dfrac{3}{64}
\left(\dfrac{R_{E}}{ R}\right)^9\left(5\left(2\,\sqrt {5}\left(g_3^2h_2^2-g_2^2h_3^2\right)+5\,\sqrt {2}\left( h_3^1g_2^1-h_2^1g_3^1\right)
\right)\cos^4 (i)
\right.
\\
\left.
+6\left(
2\,\sqrt {5}\left(g_3^2h_2^2-g_2^2h_3^2\right)+7\sqrt {2}\left(h_2^1g_3^1-h_3^1g_2^1\right)\right)\cos^2 (i)
\right.
\\
\left.
+3\left(
3\,\sqrt {2}\left(h_3^1g_2^1-h_2^1g_3^1\right)+2\sqrt {5}\,\left(g_2^2h_3^2-g_3^2h_2^2\right)
\right)
\right),
\end{multline*}
\begin{equation*}%\label{}
{\langle
B_{\xi}B_{\zeta}\rangle}_t=
\dfrac{\cos (i)\sqrt{6}}{8}
\left(\dfrac{R_{E}}{ R}\right)^8
\left(5\cos^{2}(i)-3\right)\left(h_1^1    g_3^1- h_3^1g_1^1\right).
\end{equation*}

Let us replace the matrix $\mathbf{M}$ by its averaged value $\mathbf{M}_{cp}$ in the system of equations (\ref{eq99}):
\begin{equation}\label{eq122}
\mathbf{J}\left(\begin{array}{c}\ddot{\varphi}\\
\ddot{\theta}\\\ddot{\psi}\end{array}\right)+\mathbf{H}\left(\begin{array}{c}\dot{\varphi}\\
\dot{\theta}\\\dot{\psi}\end{array}\right)+\mathbf{M}_{cp}\left(\begin{array}{c}{\varphi}\\
{\theta}\\{\psi}\end{array}\right)=0.
\end{equation}

By the change
\begin{equation*}
x_1=\dfrac{d{\varphi}}{du}, \; x_2=\dfrac{d{\theta}}{du}, \; x_3=\dfrac{d{\psi}}{du}, \;
x_4={\varphi}, \; x_5={\theta}, \; x_6={\psi}
\end{equation*}
the system (\ref{eq99}) can be reduced to the dimensionless form
\begin{equation}\label{eq101}
\frac{d\mathbf{x}}{du}=(\mathbf{N}+\widetilde{\mathbf{N}}(u))\mathbf{x},
\end{equation}
where the matrix $\mathbf{N}$ have the block structure:
\[\mathbf{N}=\left(
\begin{array}{c|c}
{-\dfrac{\mathbf{J}^{-1}\mathbf{H}}{\omega_0}}  & {-\dfrac{\mathbf{J}^{-1}\mathbf{M_{cp}}}{\omega_0^2}}  \\ \\[-0.3cm]
\hline \\ \\[-1.0cm]
 \begin{array}{ccc}1 & 0 & 0 \\
 0 & 1 & 0  \\
 0 & 0 & 1 \end{array} &
 \begin{array}{ccc} 0 & 0 & 0 \\
 0 & 0 & 0 \\
  0 & 0 & 0\end{array}
\end{array}\right)
\]

All solutions of the system $\dfrac{d\mathbf{x}}{du}={\mathbf{N}}\mathbf{x}$, corresponding to (\ref{eq122}), are asymptotically stable if and only if all the real parts of the roots of characteristic equation
$\mathrm{det}(\lambda I-\mathbf{N})=0$ are negative. The numerical computer analysis showed that there exists a domain of parameters of the Earth's artificial satellite and its orbit, in the presence of which this condition is fulfilled, and, hence, in this domain, real parts $\alpha_j$ of eigenvalues of the matrix $\mathbf{N}$ satisfy the inequalities $\alpha_j<-\mu$,
where $\ j=1..6,\ \mu=\mathrm{const}>0.$

For visualization of such domain it is advisable to construct its section fixing some parameters of spacecraft and its orbit.
For example the section of mentioned domain at the following fixed parameters
\begin{equation}\label{eq155}
\begin{array}{c}
R=7\cdot10^6\; \mbox{m,}\
A=10^3\; \mbox{kg}\cdot\mbox{m}^2,\ Q=5 \cdot 10^{-3}\; \mbox{C,}\ \\
\ h_1=h_2=h_3=h=0.5A\omega_0\; \mbox{kg}\cdot\mbox{m}^2/\mbox{s}, \ k_{L}=10\; \mbox{s}/\mbox{T}, \ k_{M}=2 \cdot 10^6\; \mbox{N}\cdot\mbox{m}/\mbox{T}^2,
\end{array}
\end{equation}
by the plain of inertial parameters $(\delta, \ \varepsilon)$ at $i=\pi/3\, \mbox{rad} \, $ is plotted in the Fig.~\ref{Fig45}.
\begin{figure}
\centering \includegraphics[width=8cm]{f45.eps}
\caption{} \label{Fig45}
\end{figure}

Moreover, the comparison of the given domain of parameters, obtained with the joint use of the moments $\vec M_{M}$ and $\vec M_{L}$ (Fig.~\ref{Fig45}), with corresponding domains, constructed for the cases, when only one of these moments is in use, was fulfilled. All calculations were performed for parameters of spacecraft and its orbit, cited in equalities
(\ref{eq155}), except the cases of absence of one of control moments $\vec M_{M}$ or $\vec M_{L}$, for which $k_{M}=0$ or $ k_{L}=0$ accordingly. The results are shown in the Fig.~\ref{Fig46}.
\begin{figure}[!h]
\centering
\includegraphics[width=8cm]{f46.eps}
\caption{} \label{Fig46}
\end{figure}

The horizontal shading mark the region corresponding to the use of only one control moment $\vec M_{L}$,
the vertical shading mark the region corresponding to the use of only one control moment $\vec M_{M}$, the grey color mark the region which arises only thanks to the joint use of the both control moments.

According to \cite{Martinuk1989}, let us consider the function
\begin{equation}\label{eq124}
\chi(u)=\exp\left(-\frac{\mu}{2}u\right)\sum\limits_{k=0}^{5}\frac{\left(12\max\limits_{i,j}|n_{i,j}|\right)^k}{k!}u^k\
.
\end{equation}

The maximum value of this function on semiaxis [$0,\,+\infty$) we denote as $D$. If
\begin{equation}\label{eq102}
\|\widetilde{\mathbf{N}}(u)\|=\max\limits_j\sum\limits_{i=1}^6|\widetilde{{n}}_{ij}|<\dfrac{\mu}{2}D,
\end{equation}
then the zero solution of the system (\ref{eq101}) (the same system in dimension form is (\ref{eq99})) is asymp\-totically stable \cite{Martinuk1989}. Hence, the zero solution of initial nonlinear system (\ref{eq6}) will be total stable \cite{Malkin1966}.

It is easy to verify that in domain of parameters of spacecraft and its orbit where the real parts of all eigenvalues of the matrix $\mathbf{N}$ are negative, the inequality (\ref{eq102}) is valid.
As is obvious from the foregoing, the suggested method of spacecraft attitude stabilization is operable in the mentioned domain and this is confirmed by computer modeling.

\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f47.eps}
\caption{} \label{Fig47}
\end{figure}
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f48.eps}
\caption{} \label{Fig48}
\end{figure}
\nopagebreak
\begin{figure}[!h]
\centering \includegraphics[width=10cm]{f49.eps}
\caption{} \label{Fig49}
\end{figure}

In Fig. \ref{Fig47} --- Fig. \ref{Fig49} are represented the computational results of spacecraft stabilized attitude motion in the presence of the following values of parameters of the spacecraft and its orbit:
$R=7\cdot10^6\; \mbox{m,}\ i=\pi/3\; \mbox{rad,}\ A=10^3\; \mbox{kg}\cdot\mbox{m}^2,\ Q=5\cdot10^{-3}\; \mbox{C,}\
\delta=0.7,\ \varepsilon=0.7,\ h_1=h_2=h_3=h=0.5A\omega_0\; \mbox{kg}\cdot\mbox{m}^2/\mbox{s}, \ k_{L}=10\; \mbox{s}/\mbox{T}, \ k_{M}=2\cdot 10^6 \; \mbox{N}\cdot\mbox{m}/\mbox{T}^2$. At such choice of parameter $\delta$ the gravitational moment $\vec M_{G}$ is disturbing from the point of view of stabilization process in the direct equilibrium position. As initial values,  $\varphi (0)=0.2\, \mbox{rad}$,\; $\psi (0)=-0.2\, \mbox{rad}$,\; $\theta (0)=0.2\, \mbox{rad}$,\; $\omega_x(0)=0.1\omega_0\, \mbox{rad/s}$,\; $\omega_y(0)=1.1\omega_0\, \mbox{rad/s}$, \; $\omega_z(0)=0.1\omega_0\,  \mbox{rad/s}.$
%\begin{figure}
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f47}
%\hfill
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f48}
%\\
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig47}}
%\hfill
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig48}}
%\end{figure}
%
%\begin{figure}
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f49}
%\hfill
%\includegraphics[width=8.0cm,%
%keepaspectratio]{f29}
%\\
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig49}}
%\hfill
%\parbox[t]{7.5cm}{\vspace{-0.1cm}\caption{} \label{Fig29}}
%\end{figure}

In Fig.~\ref{Fig47} the graphs of "airborne" angles via the argument of a latitude are plotted for spacecraft under the action of control moments $\vec M_{L}$ and $\vec M_{M}$.
In Fig.~\ref{Fig48} the graphs of quaternion components are plotted for the case of stabilized motion, and in Fig.~\ref{Fig29} are the same components in the absence of control moments $\vec M_{L}$ ($k_{L}=0$), $\vec M_{M}$ ($k_{M}=0$) and the damping moment $\vec M_{D}$ ($h=0$).
The comparison of Fig.~\ref{Fig47} and Fig.~\ref{Fig48}, Fig.~\ref{Fig29} make it clear that the suggested stabilization system permits to reach the prescribed spacecraft attitude position, and the stabilization process is rather fast as in the case of the equatorial orbit ($i=0$).
The values of the moments $\vec{M}_{L}, \  \vec{M}_{M},\ \vec{M}_{G},\ \vec{M}_{D}$  acting upon the spacecraft in the process of stabilization, are shown in Fig.~\ref{Fig49}. It is obvious that the mentioned values are in such limits that make evident the operability of suggested method of attitude control. The comparison of obtained results with similar results obtained with the use of only one of the moments $\vec M_{L}$ or $\vec M_{M}$ is made. It is revealed that the joint use of both control moments $\vec {M}_{L}$ and $\vec M_{M}$ in electrodynamical attitude stability system result in reduce of the transition time.

CONCLUSION

Thus, the possibility of spacecraft attitude stabilization on the orbits with arbitrary inclinations with the help of control moments $\vec M_{L}$, $\vec M_{M}$ and $\vec M_{D}$ is proved. It is proved and verified by computer modeling that the suggested attitude control permits to reach the stabilized motion in the short time. It should be noted that all computations
were carried out with the use of the octupole approximation of geomagnetic field in order to achieve the necessary calculation accuracy. It was revealed that the use of more simple models of geomagnetic field (dipole or quadrupole approximations) may give rise to incorrect results.

%REFERENCES

\begin{thebibliography}{99}

%\footnotesize

%\itemsep 2pt

\bibitem{Kov.1975}
 Kovalenko A.P.\, Magnitnye sistemy upravleniya kosmicheskimi letatel'nymi apparatami (Magnetic Systems for Spaceship Control), Mashinostroenie, Moscow, 1975.

\bibitem{Alp.1978}
Alpatov A.P., Dranovskii V.I., Saltykov Yu.D., Khoroshilov V.S.\, Dinamika kosmicheskikh apparatov s magnitnymi sistemami upravleniya (Dynamics of Spacecraft with Magnetic Control Systems), Mashino\-stroenie, Moscow, 1978.

\bibitem{KRes2003}
Tikhonov A.A. \, {\it A Method of Semipassive Attitude Stabilization of a Spacecraft in the Geomagnetic Field},\, Cosmic Research, 41, no.  1, (2003), 63--73.

\bibitem{AutRemCon2007en}
Antipov K.A. and Tikhonov A.A. \, {\it Parametric Control in the Problem of Spacecraft Stabilization in the Geomagnetic Field},\, Automation and Remote Control, 68, no. 8, (2007), 1333--1345.

\bibitem{Vestnik1999}
Petrov K.G. and Tikhonov A.A. \, {\it The moment of Lorentz forces, acting upon the charged satellite in the geomagnetic field. Part 1: The Earth's magnetic field strenth in the orbital coordinate system},\, Vestnik of St.-Petersburg State University, Ser.1, 1~(no.~1), (1999), 92--100.

\bibitem{KRes2002}
Tikhonov A.A. and Petrov K.G. \, {\it Multipole Models of the Earth's Magnetic Field},\,
 Cosmic Research, 40, no. 3, (2002), 219--229.

\bibitem{Vest2006}
Antipov K.A., Petrov K.G., Tikhonov A.A. \,
{\it The electrodynamical method for attitude stabilization of a dynamically symmetric spacecraft}, Vestnik of
St.-Petersburg State University. Ser.1, 1~(no.~1), (2006), 75--86.

\bibitem{Beletsky1975}
Beletskii V.V. \, Dvizhenie sputnika otnositel'no tsentra mass v gravitatsionnom pole (Satellite Motion with Respect to the Center of Mass in the Gravitational Field),Mosk. Gos. Univ., Moscow, 1975.

\bibitem{Martinuk1989}
Martynyuk, A.A., Lakshmikantam, V., and Lila, S. \, Ustoichivost' dvizheniya: metod integral'nykh neravenstv (Stability of Motion: The Integral Inequality Method), Naukova Dumka, Kiev, 1989.

\bibitem{BR 1954}
Bellmann R. \, The stability theory of differential equations solutions, 1954.

\bibitem{Malkin1966}
Malkin I.G. \, Teoriya ustoichivosti dvizheniya (Theory of Motion Stability), Nauka, Moscow, 1966.
\end{thebibliography}
%\newpage

\vspace*{1.5cm}

\begin{center}
{\normalsize Submitted in July 2011.}
\end{center}

\clearpage

\vspace*{1cm}

\begin{center}
\textbf{\large O stabilizaciji letilice u orbitalnom koordinatnom sistemu}
\end{center}

\medskip

\noindent

U radu se razmatra letilica u kru\v znoj orbiti bliskoj zemlji. Letilica medjudejstvuje sa geomagnetnim poljem momentima Lorencove i magnetne sile. Prihva\' cena je oktupolna aproksimacija zemljinog magnetnog polja. Elektromagnetni parametri letilice, naime, elektrostati\v cki moment naelektrisanja prvog reda kao i sopstveni magnetni moment su kontrolisane kvaziperiodi\v cne funkcije. Upravljaju\' ci algoritmi za elektromagnetne parametre letilice koji dozvoljavaju stabilizovanje polo\v zaja letilice u orbitalnom koordinatnm sistemu su dobijeni. Stabilnost stabilizovane orijentacije letilice je dokazana kako analiti\v cki tako i PC izra\v cunavanjima.

\vspace*{8cm}

\begin{center}
\rule[5mm]{12.5cm}{0.2mm}\\[-4.5mm]
% \rule[15mm]{12cm}{0.15mm}\\[-10mm]

% \vspace*{-3mm}

{\scriptsize doi:10.2298/TAM1202127A \qquad \qquad \qquad
\qquad \qquad Math.Subj.Class.: 70E50, 70Q05, 34H15}
\end{center}
\end{document}

70E50 Stability problems
70Q05 Control of mechanical systems
34H15 Stabilization

