# Local positioning system as a classic alternative to atomic navigation.

B. Dubetsky  Independent Researcher, 1849 S Ocean Dr, Apt 207, Hallandale
Florida 33009, United States<EMAIL_ADDRESS>
###### Abstract

A local positional system (LPS) is proposed, in which particles are launched
at given velocities, and a sensor system measures the trajectory of particles
in the platform frame. These measurements allow us to restore the position and
orientation of the platform in the frame of the rotating Earth, without
solving navigation equations. When the platform velocity is known and if the
platform orientation stays the same, the LPS-technique allows a navigational
accuracy of 100$\mu$ per one hour to be achieved . In this case, the LPS
technique is insensitive to the type of platform trajectory. If there are also
velocimeters installed on the platform, then one can restore the velocity and
angular rate of the platform rotation in respect to the Earth. Instead of
navigational equations, it is necessary to obtain the classical trajectory of
a particle in the field of a rotating gravity source. Taking into account the
gravity-gradient, Coriolis, and centrifugal forces, the exact expression for
this trajectory is derived, which can be widely used in atomic interferometry.
A new iterative method for restoring the orientation of the platform without
using gyroscopes is developed. The simulation allowed us to determine the
conditions under which the LPS-navigation error per hour is about $10$m.

###### pacs:

03.75.Dg; 37.25.+k; 04.80.-y

##  I Introduction

Since its birth about 40 years ago c1 , the field of atom interferometry has
matured significantly. The current state and prospects in this area are
presented, for example, in the reviews c1.1  and the proposals c1.2 ; c1.3 ;
c1.4 ; c1.5 ; c1.5.1 ; c1.5.2 .

Among other applications, atom interferometers (AIs) are now proposed to use
as inertial measurement units (IMUs) in navigation.

The precision accuracy of AI allows us to hope for significant progress,
improving the accuracy and duration of the navigation. These applications are
based on the possibility of using AI as an accelerometer c2 . If the atoms
move with acceleration $a$, then the AI phase is given by

$\phi=\mathbf{k}\cdot\mathbf{a}T^{2},$ (1)

where $\mathbf{k}$ is the effective wave vector of the Raman pulses, $T$ is
the delay time between these pulses. The error of the acceleration measurement
is

$\delta a=\phi_{err}/kT^{2},$ (2)

where $\phi_{err}$ is the accuracy of the AI phase measurement. To date, for
example, the accuracy of differential acceleration measurements has reached

$\delta a=1.4ppt$ (3)

at the delay time

$T=955ms.$ (4)

For navigation applications, a system of 6 AIs was created and tested c4 .
However, the direct use of precise AI as an IMU is difficult. Matter of fact
is that the IMUs data should be substituted into the navigation equations, and
for an accurate numerical solution of these equations, the IMUs measurements
should be repeated with a time step of at least $\tau\sim 1$ms c5 . To do
this, the interrogation time of the AI must be no more than $\tau$, and
therefore $T\lesssim\tau/2$. Thus, with the same accuracy of the phase
measurement $\phi_{err}$ and the effective wave vector $k$, the acceleration
measurement error increases by 6 orders of magnitude to the level of $5$ppm.
With such an inaccuracy, atomic interferometers will be even an order of
magnitude worse than conventional accelerometers c6  and will not result in
any progress in navigation.

Two methods have been proposed to avoid this difficulty. In one method, AIs is
used in a hybrid mode together with conventional IMUs c7 ; c8 ; c9 ; c9.1 ;
c10 . Another option c11  is related to the fact that one can restore the
position of the platform directly from the AI phase. without measuring the
atomic acceleration and, in principle, without any use of conventional IMUs.
To understand why this is possible, one notes that Eq. (1) is valid only for
the uniform accelerated motion of the atom. If the atomic cloud moves in the
platform frame, then, ignoring the quantum effects, for the phase instead of
(1) one gets

$\displaystyle\phi$ $\displaystyle=$ $\displaystyle\mathbf{k}\cdot\mathbf{p},$
(5a) $\displaystyle\mathbf{p}$ $\displaystyle=$
$\displaystyle\mathbf{x}\left(t_{3}\right)-2\mathbf{x}\left(t_{2}\right)+\mathbf{x}\left(t_{1}\right),$
(5b) where $\mathbf{x}(t)$ is the classical trajectory of the cloud in the
platform frame. Knowing the trajectory at the previous moments $t_{1}$ and
$t_{2}=t_{1}+T$ and measuring the phase of AI, one restores the trajectory at
the moment $t_{3}=t_{1}+2T$. At the same time, there is no need to solve
navigation equations, and as a result, the restriction on the feasible delay
time $T$ is removed. One sees that when using AI as an IMU, one actually does
double work: first, by measuring the trajectory,
$\mathbf{x}\left(t_{3}\right)$, gets the acceleration, and then from the
navigation equations gets the trajectory.

The error of measuring the position $\mathbf{x}\left(t_{3}\right)$ is given by

$\delta_{x}\sim\phi_{err}/k\sim 1.3\text{pm,}$ (6)

where to estimate the parameter $\phi_{err}/k$, one can use the Eqs. (3-4).
Despite the fact that at each step the error (6) does not depend on $T$, the
total error $\Delta X$ during navigation time $t_{n}$ accumulates after
repeated measurement of phases at times multiple to $T$, . For the given
$t_{n}$, the number of steps

$n\sim t_{n}/T.$ (7)

The use of AIs with a longer interrogation time leads to a reduction in the
number of steps and, consequently, to a reduction in the accumulated error.
Thus, one sees that, proposed in c11 , the rejection of the use of AIs as
accelerometers, allows the use of an unprecedently high accuracy of trajectory
measurement (6) for an unprecedently accurate navigation.

Atom interference is an exclusively quantum phenomenon. Nevertheless, the
expression (5) for the phase is purely classical. This means that instead of
quantum objects, atoms, one can use classical material points, which in the
future we will call particles. We do not specify the kind of particles.
However, we assume that the mass of the particle $M$ is sufficiently large to
neglect the recoil velocity $\hbar k/M$, where, in our case, $\hbar k$ is the
change in momentum that accompanies the measurement of the particle position.

Like in the patent c11 , in this article, we propose to abandon the use of
accelerometers, gyroscopes, and navigation equations, and instead measure the
trajectories of particles that are launched in the platform system from
certain points, at certain speeds, and with a time interval $T$. We will
develop a theory of such navigation here and perform numerical simulations.
The obvious advantages of classical objects are a dramatic mitigation in the
requirements of high vacuum, disappearing the necessity to solve the problem
of the phase multiple ambiguity, and the absence of a recoil effect. The
obvious disadvantage of our method is that for classical objects, the accuracy
of the measurement (6) of their position is hardly achievable at the present
time. We have, however, developed methods here that will significantly relax
the requirements for $\delta_{x}$ The ultimate goal of numerical simulation is
to answer the question of what the accuracy of measuring the coordinates and
velocities of particles should be in order to achieve a given accuracy of
navigation $\Delta X$ for a given period of time $t_{n}$.Specifically, we will
find the answer to this question for

$\displaystyle t_{n}$ $\displaystyle=$ $\displaystyle 1\text{h,}$ (8a)
$\displaystyle\Delta X$ $\displaystyle=$ $\displaystyle 10\text{m.}$ (8b)

The accelerometer in conventional navigation measures the acceleration of
particle $\mathbf{a}$ in the platform system. For $\mathbf{a}$ in the absence
of rotation one has

$\mathbf{a}=\mathbf{g}-\mathbf{A},$ (9)

where $\mathbf{g}$ and $\mathbf{A}$ are the gravitational field of the Earth
and the acceleration of the platform. Knowing $\mathbf{g}$ and measuring
$\mathbf{a}$, one obtains the acceleration of platform $\mathbf{A}$,
substituting it into the navigation equation

$\mathbf{\ddot{X}}\left(t\right)=\mathbf{A},$ (10)

to determine the trajectory of the platform $\mathbf{X}\left(t\right)$.

In Ref. c11 , it is actually proposed to replace the accelerations
$\mathbf{a}$ and $\mathbf{A}$, with the trajectories $\mathbf{x}(t)$ and
$\mathbf{X}\left(t\right)$. The procedure for determining the position of the
platform by measuring the position of the particle in the absence of rotation
of the platform is illustrated in Fig. 1.

Figure 1: Exploiting the measurement of the particle position in the platform
frame for. navigation. The trajectories of the platform and the particle are
shown with solid and dotted curves.

We are interested in the trajectory of the platform $\mathbf{X}\left(t\right)$
in the coordinate system $\left\\{X,Y,Z\right\\}$ rotating with the Earth. On
the platform, one installs a local frame $\left\\{x,y,z\right\\}$ with the
origin at the point $\mathbf{X}\left(t\right)$ and sensors that measure the
position of the particle $\mathbf{x}(t)$ in that local frame. At the moment of
time $t_{1}$ one launches the particle. After the launch, the particle is
completely decoupled from the platform and moves along the trajectory
$\mathbf{X}_{E}\left(t\right)$ in a vacuum under the influence of only the
rotating Earth. Knowing the potential of the Earth’s gravitational field and
the velocity of the Earth rotation $\vec{\Omega}_{E}\left(t\right),$ one can
obtain the trajectory $\mathbf{X}_{E}\left(t\right)$. At time $t_{2}$ one
measures the position of the particle $\mathbf{x}$. Then from Fig. 1 one gets
for the platform trajectory

$\mathbf{X}\left(t_{2}\right)=\mathbf{X}_{E}\left(t_{2}\right)-\mathbf{x}.$
(11)

The position of the particle at the time of launch in the system
$\left\\{X,Y,Z\right\\}$, obviously coincides with the position of the
platform, $\mathbf{X}_{E}\left(t_{1}\right)=\mathbf{X}\left(t_{1}\right)$. To
determine the trajectory $\mathbf{X}_{E}\left(t\right)$, it is also necessary
to get the initial velocity of the particle,
$\mathbf{v}_{E}\left(t_{1}\right)$, and it depends on the velocity of the
platform $\vec{V}\left(t_{1}\right)$, the velocity of the particle
$\mathbf{v}$ in the platform frame, the initial orientation of the platform,
the rotation matrix $R\left(t_{1}\right)$, and the initial rate of the
platform rotation $\mathbf{\Omega}\left(t_{1}\right)$. It is also necessary to
take into account that at the moment $t_{2}$ the orientation of the system
changes , the rotation matrix $R\left(t_{2}\right)\not=R\left(t_{1}\right)$
and one must rotate the vector $\mathbf{x}$ from the frame
$\left\\{x,y,z\right\\}$ in which it is measured to the frame
$\left\\{X,Y,Z\right\\}$. Below we will derive an obvious generalization of
the formula (11) taking into account all these factors. But here we we would
like to note that for successful navigation, restoring the position of the
platform, it is also necessary to restore the parameters

$Z=\left\\{\mathbf{V}\left(t\right),\mathbf{\Omega}\left(t\right),R\left(t\right)\right\\}.$
(12)

However, in order to demonstrate the power of the navigation method we are
considering here, let’s assume that these parameters are known exactly. We
also assume that the particles are launched with a step $T=300$ms, and the
position of the particle is measured with accuracy [standard deviation (SD)]

$\delta_{x}=1\mu.$ (13)

Then the time dependence of the navigation error is shown in Fig.2.

Figure 2: The green curves are the time dependences of the navigation error
$\Delta X\left(t\right)$ for an ensemble of 30 randomly generated platform
trajectories $\mathbf{X}\left(t\right)$. The red curve is the average error in
the ensemble $\left\langle\Delta X\left(t\right)\right\rangle$, the black
curve is SD $\sqrt{\left\langle\left[\Delta
X\left(t\right)\right]^{2}\right\rangle-\left\langle\Delta
X\left(t\right)\right\rangle^{2}}.$

We assume that the accuracy of measuring the trajectory of a classical
particle (13) is 6 orders of magnitude worse than the accuracy of measuring
(6) the trajectory of a quantum atomic cloud, achievable by atomic
interferometry methods c1 . And nevertheless, ideally, the accuracy of
navigation in $1$ hour does not exceed $150{\mu}$ and practically remains the
same for different trajectories.

Following c6 , in the future we will assume that the north-east-down (NED)
local frame is installed on the platform. Likewise GPS, we propose to install
in NED-frame a sensor system to measure the distance to the particle, which we
will call the local positioning system (LPS), see Fig. 3.

Figure 3: LPS - measurement of the particle’s coordinates. From the distances
to three red satellites, one determines two positions of the particle, true
white and false red. Similarly, for green satellites, one gets true white and
false green. The true solution is obviously the same for the red and green
satellites.

As in GPS, we will call these sensors satellites. Our code can be used for any
satellite location. However, in the simulation, we use only satellites located
on the frame axes. When three satellites are located on the positive semi-axes
of the system, the vector connecting the particle and the satellites do not
lie in the same plane, and the coordinates of the vector $\mathbf{x}$ can be
obtained from the distances to the satellites. The equations for $\mathbf{x}$
have two solutions, one true, the other false. To overcome this ambiguity, one
can use an additional system of satellites located, for example, symmetrically
on negative semi-axes and choose as true the solution that is the same for
both satellite systems. During the simulation, pseudo-random corrections were
made to the exact distances using the normal random number generator c12  with
SD $\delta_{x}$.

At the launching time, we need to restore the values of the parameters (12).
Errors in restoring these parameters degrade navigation accuracy by several
orders of magnitude. These parameters can be restored using conventional IMUs.
In this case, we come to another kind of hybrid. Our code can be included in
the code of operation of such a hybrid. In this paper, however, we will
consider the possibilities of reconstructing $Z$ by launching particles with
different velocities and simultaneously measuring their coordinates and
velocities through time T.

The most significant factor affecting the accuracy of navigation is the
platform rotation. In the case of atomic navigation, it was proposed to use
the gimbal c13  to reduce the effect of rotation. In this paper, we also
assume the presence of gimbal, which we define as reference systems that
perform only noisy rotations at small angles, with zero mean value and SD
$\delta_{g}$. In the future, we will determine the quality of the gimbals, SD
$\delta_{g}$, to achieve the ultimate goal of navigation (8).

An elementary loop of the navigation is the restoration of the position,
velocity, orientation and rate of rotation of the platform at the time
$t_{2}=t_{1}+T$ from their values at the time of the particle launch $t_{1}$.
Then the restored parameters at the time $t_{2}$ are used for navigation at
the time $t_{3}=t_{2}+T$ and so on. We will call this method first order
navigation. Due to the large Doppler broadening of the spectral lines, one
could not use both the atomic interferometry c1  and the Ramsey fringes c14
for this method. We propose installing in the NED-frame velocimeters together
with satellites. Measuring the velocity of a particle in the platform frame
will restore the velocity of the platform $V\left(t_{2}\right),$ while
measuring the distance difference between particles launched at different
velocities will restore the orientation of the platform $R\left(t_{2}\right),$
and finally, measuring the velocity difference of these particles will restore
the rotation rate of the platform $\mathbf{\Omega}\left(t_{2}\right).$

The severe requirements on the platform velocity recovery can be mitigated in
the case of the 2nd order navigation, when one measures particle’s positions
at three times $\left\\{t_{1},t_{2}=t_{1}+T,t_{3}=t_{1}+2T\right\\}$ and
employs the 2nd order difference (5b) to restore platform position
$\mathbf{X}\left(t_{3}\right)$. The reason here is that even with an arbitrary
trajectory of the platform, the difference (5b) does not depend on the initial
velocity of the particle, if the platform does not rotate, and the Earth’s
gravitational field is uniform. The weak inhomogeneity of the field and the
noisy, small rotation of the gimbal lead to such a weak dependence of
$\mathbf{p}$ on the platform velocity that there is no need in velocimeters.

The requirements for the quality of the gimbal can also be relaxed if one
launches simultaneously 3 particles with different velocities and restores the
orientation of the gimbal from the differences in their trajectories. In the
case of the atom interferometry, iterative method of this kind restoration was
developed in c13 . Here we generalize this method to the case of LPS-
navigation (LPS-N).

We offer here 5 varieties of LPS-N:

  1. 1.1

1st order LPS-N, the quality of the gimbal is high enough to neglect its
rotation;

  2. 1.2

1st order LPS-N, the quality of the gimbal is sufficient to neglect its
rotation, but one does not neglect the rate of rotation, and restores it
applying data of velocimeters;

  3. 1.3

1st order LPS-N, one restores the orientation and rotation rate of the gimbal
applying satellites’ and velocimeters’ data;

  4. 2.1

2nd order LPS-N, the quality of the gimbal is high enough to neglect its
rotation;

  5. 2.2

2nd order LPS-N, one restores the orientation of the gimbal from satellite
data.

The article is arranged as follows. The exact expression for the trajectory of
a particle in the field of the rotating Earth, when the potential of this
field is expanded to the 1st order gravity-gradient terms, is obtained in c11
. For the completeness of this study, we derive this expression in the next
section. The technique of restoring the orientation, rotation rate, velocity,
and then the position of the platform is described in Sec. III for LPS-N of
the 1st order and in Sec. IV for LPS-N of the 2nd order. The results of the
simulation and their discussion were carried out in the Sec. V. In the
appendix Sec. C we describe the model of the Earth’s gravitational field
accepted in our code.

##  II Particle trajectory

If the source of gravity rotates at a constant velocity $\mathbf{\Omega}_{E}$,
then the trajectory of the particle obeys the following equation of motion

$\overset{{}_{\bullet\bullet}}{\mathbf{X}}_{E}=-\mathbf{\triangledown}U\left(\mathbf{X}_{E}\right)-\mathbf{\Omega}_{E}\times\left(\mathbf{\Omega}_{E}\times\mathbf{X}_{E}\right)-2\mathbf{\Omega}_{E}\times\mathbf{v}_{E},$
(14)

where $\mathbf{v}_{E}=\overset{{}_{\bullet}}{\mathbf{X}}_{E}$ is the velocity
of the particle in the rotating Earth frame, $U\left(\mathbf{X}\right)$ is the
potential of the gravitational field. If the motion occurs in a small vicinity
of the point $\mathbf{X}_{c}$, then expanding the potential to the second
order terms,

$U\left(\mathbf{X}_{E}\right)\approx
U\left(\mathbf{X}_{c}\right)-\mathbf{g}_{E}\left(\mathbf{X}_{c}\right)\cdot\mathbf{x}_{E}-\dfrac{1}{2}\mathbf{x}_{E}^{T}\underline{\Gamma}\left(\mathbf{X}_{c}\right)\mathbf{x}_{E},$
(15)

where

$\mathbf{x}_{E}=\mathbf{X}_{E}-\mathbf{X}_{c},$ (16)

and

$\displaystyle\mathbf{g}_{E}\left(\mathbf{X}\right)$ $\displaystyle=$
$\displaystyle-\vec{\triangledown}U\left(\mathbf{X}\right),$ (17a)
$\displaystyle\underline{\Gamma}\left(\mathbf{X}\right)$ $\displaystyle=$
$\displaystyle-\mathbf{\partial}_{\mathbf{X}}\mathbf{\partial}_{\mathbf{X}}^{T}U\left(\mathbf{X}\right)$
(17b) are the gravitational acceleration and the 1st order gravity-gradient
tensor, one comes to the equation of motion
$\overset{{}_{\bullet\bullet}}{\mathbf{x}}_{E}=\mathbf{g}\left(\mathbf{X}_{c}\right)-2\mathbf{\Omega}_{E}\times\mathbf{v}_{E}+\underline{q}\left(\mathbf{X}_{c}\right)\mathbf{x}_{E},$
(18)

where

$\displaystyle\mathbf{g}\left(\mathbf{X}_{c}\right)$ $\displaystyle=$
$\displaystyle\mathbf{g}_{E}\left(\mathbf{X}_{c}\right)-\mathbf{\Omega}_{E}\times\left(\mathbf{\Omega}_{E}\times\mathbf{X}_{c}\right),$
(19a) $\displaystyle\underline{q}\left(\mathbf{X}_{c}\right)$ $\displaystyle=$
$\displaystyle\underline{\Gamma}\left(\mathbf{X}_{c}\right)-\underline{\Omega}_{E}^{2},$
(19b) and the tensor $\underline{\Omega}_{E}$ is given by
$\underline{\Omega}_{Eik}=-\varepsilon_{ikm}\Omega_{Em}.$ (20)

Despite the fact that Eq. (18) is a linear system of equations with permanent
coefficients, most of the papers used expansions of
$\mathbf{x}_{E}\left(t\right)$ into power series over time. In the Refs. c14.1
; c14.2 , the terms up to the 4th order were calculated, the 5th and 6th order
terms were obtained in Ref. c14.3 , the expansion to the 7th order was
obtained in the article c14.4 . The exact solution of the equation (18) is
obtained in the patent c11 . Here we derive exact solutions for arbitrarily
large values of the tensor $\underline{q}\left(\mathbf{X}_{c}\right)$ and the
rate of rotation of the gravity field source, $\mathbf{\Omega}_{E}$. In this
case, we are using the method proposed in c15 . Consider first the homogeneous
equation $\left[\mathbf{g}\left(\mathbf{X}_{c}\right)=0\right]$. Its solution
is $\mathbf{x}_{E}\left(t\right)=\mathbf{x}\left(\omega\right)e^{-i\omega t}$,
where $\omega$ and $\mathbf{x}\left(\omega\right)$ are the eigenvalue and
eigenvector of the matrix

$\underline{A}\left(\omega\right)=\underline{q}\left(\mathbf{X}_{c}\right)+\omega^{2}I+2i\omega\underline{\Omega}_{E},$
(21)

where $I$ is the unity matrix. One notices that owing the symmetry of the
tensor $\underline{q}\left(\mathbf{X}_{c}\right),$ the determinant
$\left|\underline{A}\left(\omega\right)\right|$ is an even function $\omega$,
$\left|\underline{A}\left(\omega\right)\right|=\left|\underline{A}\left(-\omega\right)\right|$,
and, therefore, the characteristic.equation

$\left|\underline{A}\left(\omega\right)\right|=0,$ (22)

is reduced to the 3rd order equation for $\omega^{2}$

$\omega^{6}+a_{2}\omega^{4}+a_{1}\omega^{2}+a_{0}=0,$ (23)

where

$\displaystyle a_{0}=\left|q\left(\mathbf{X}_{c}\right)\right|,$ (24a)
$\displaystyle
a_{1}=\dfrac{1}{2}\left\\{Tr^{2}\left[\underline{q}\left(\mathbf{X}_{c}\right)\right]-Tr\left[\underline{q}^{2}\left(\mathbf{X}_{c}\right)\right]\right\\}-4\mathbf{\Omega}_{E}^{T}\underline{q}\Omega_{En},$
(24b) $\displaystyle
a_{2}=Tr\left[\underline{q\left(\mathbf{X}_{c}\right)}\right]-4\Omega_{E}^{2}.$
(24c) Eq. (23) has 6 roots  $\displaystyle\omega_{i}$ $\displaystyle=$
$\displaystyle\sqrt{z_{i}-\dfrac{a_{2}}{3}},~{}\text{for }i=1,2,3$ (25a)
$\displaystyle\left\\{\omega_{4},\omega_{5},\omega_{6}\right\\}$
$\displaystyle=$
$\displaystyle-\left\\{\omega_{1},\omega_{2},\omega_{3}\right\\},$ (25b) where
$z_{i}$ is the root of the equation $z^{3}+pz+s=0$ and  $\displaystyle p$
$\displaystyle=$ $\displaystyle a_{1}-\dfrac{a_{2}^{2}}{3},$ (26a)
$\displaystyle s$ $\displaystyle=$ $\displaystyle
a_{0}-\dfrac{a_{2}a_{1}}{3}+\dfrac{2a_{2}^{3}}{27}.$ (26b) The eigenvectors
are given by
$\mathbf{x}\left(\omega\right)=\left(\begin{array}[]{c}A_{22}\left(\omega\right)A_{33}\left(\omega\right)-A_{23}\left(\omega\right)A_{32}\left(\omega\right)\\\
A_{23}\left(\omega\right)A_{31}\left(\omega\right)-A_{21}\left(\omega\right)A_{33}\left(\omega\right)\\\
A_{21}\left(\omega\right)A_{32}\left(\omega\right)-A_{22}\left(\omega\right)A_{31}\left(\omega\right)\end{array}\right).$
(27)

Let’s return now to the equation (18). We introduce a 6-dimensional vector

$\mathbf{y}=\left(\begin{array}[]{c}\mathbf{x}_{E}\\\
\mathbf{v}_{E}\end{array}\right),$ (28)

which obeys the equation

$\mathbf{\dot{y}}=\left(\begin{array}[]{cc}0&I\\\
\underline{q}&-2\underline{\Omega}_{E}\end{array}\right)\mathbf{y}+\left(\begin{array}[]{c}0\\\
\mathbf{g}\end{array}\right).$ (29)

One gets that solution is given by

$\displaystyle\mathbf{y}\left(t\right)$ $\displaystyle=$
$\displaystyle\underline{\Phi}\left(t\right)\mathbf{y}\left(0\right)+\int_{0}^{t}dt^{\prime}\underline{\Phi}\left(t-t^{\prime}\right)\left(\begin{array}[]{c}0\\\
\mathbf{g}\end{array}\right),$ (32)
$\displaystyle\underline{\Phi}\left(t\right)$ $\displaystyle=$
$\displaystyle\underline{Y}\left(t\right)\underline{Y}^{-1}\left(0\right),$
(33)

where $6\times 6$ matrix $\underline{Y}\left(t\right)$ is composed of
eigenvectors,

$\underline{Y}\left(t\right)=\left\\{\mathbf{y}_{1}e^{-i\omega_{1}t},\ldots,\mathbf{y}_{6}e^{-i\omega_{6}t}\right\\}.$
(34)

Since

$\mathbf{y}_{m}=\left(\begin{array}[]{c}\mathbf{x}_{m}\\\
-i\omega_{m}\mathbf{x}_{m}\end{array}\right),$ (35)

one can represent the matrix (34) as

$\underline{Y}\left(t\right)=\left(\begin{array}[]{cc}\underline{x}_{+}e^{-i\underline{\tilde{\omega}}t}&\underline{x}_{-}e^{i\underline{\tilde{\omega}}t}\\\
-i\underline{x}_{+}\underline{\tilde{\omega}}e^{-i\underline{\tilde{\omega}}t}&i\underline{x}_{-}\underline{\tilde{\omega}}e^{i\underline{\tilde{\omega}}t}\end{array}\right),$
(36)

where $3\times 3-$matrices $\underline{x}$± and $\omega$ are composed of
eigenvectors (27) and eigenfrequencies

$\displaystyle\underline{x}_{\pm}$ $\displaystyle\equiv$
$\displaystyle\left\\{\mathbf{x}\left(\pm\omega_{1}\right),\mathbf{x}\left(\pm\omega_{2}\right),\mathbf{x}\left(\pm\omega_{3}\right)\right\\},$
(37a) $\displaystyle\underline{\tilde{\omega}}$ $\displaystyle=$
$\displaystyle diag\left(\omega_{1},\omega_{2},\omega_{3}\right).$ (37b) To
inverse the matrix $\underline{Y}\left(0\right)$, one can use the method
developed in c16 . It is convenient to represent the $6\times 6-$matrices
$\underline{\Phi}\left(t\right)$ and $\underline{Y}^{-1}\left(0\right)$ as
$2\times 2-$matrices  $\displaystyle\underline{\Phi}\left(t\right)$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\underline{x}_{x}\left(t\right)&\underline{x}_{v}\left(t\right)\\\
\underline{\dot{x}}_{x}\left(t\right)&\underline{\dot{x}}_{v}\left(t\right)\end{array}\right),$
(38c) $\displaystyle\underline{Y}^{-1}\left(0\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}b_{11}&b_{12}\\\
b_{21}&b_{22}\end{array}\right),$ (38f) in which the matrix elements are
$3\times 3-$matrices. Then one gets for $b_{mn}$
$\displaystyle\underline{x}_{+}b_{11}+\underline{x}_{-}b_{21}$
$\displaystyle=$ $\displaystyle I,$ (39a)
$\displaystyle\underline{x}_{+}\underline{\tilde{\omega}}b_{11}-\underline{x}_{-}\underline{\tilde{\omega}}b_{21}$
$\displaystyle=$ $\displaystyle 0,$ (39b)
$\displaystyle\underline{x}_{+}\underline{\tilde{\omega}}b_{12}-\underline{x}_{-}\underline{\tilde{\omega}}b_{22}$
$\displaystyle=$ $\displaystyle i,$ (39c)
$\displaystyle\underline{x}_{+}b_{12}+\underline{x}_{-}b_{22}$
$\displaystyle=$ $\displaystyle 0.$ (39d) Solving the Eqs. (39) and
substituting their solutions consequently in the Eqs. (33) and (32), after
integrating into the Eq. (32), one arrives at the next exact solution of the
Eq. (18)  $\displaystyle\mathbf{x}_{E}\left(t\right)$
$\displaystyle=\underline{x}_{x}\left(t\right)\mathbf{x}_{E}\left(0\right)+\underline{x}_{v}\left(t\right)\mathbf{v}_{E}\left(0\right)+\underline{x}_{g}\left(t\right)\mathbf{g}\left(\mathbf{X}_{c}\right),$
(40a) $\displaystyle\underline{x}_{x}\left(t\right)$
$\displaystyle=\left(\underline{x}_{+}e^{-i\underline{\tilde{\omega}}t}+\underline{x}_{-}e^{i\underline{\tilde{\omega}}t}\underline{\tilde{\omega}}^{-1}\underline{x}_{-}^{-1}\underline{x}_{+}\underline{\tilde{\omega}}\right)\underline{x}_{i},$
(40b) $\displaystyle\underline{x}_{v}\left(t\right)$
$\displaystyle=i\left(\underline{x}_{+}e^{-i\underline{\tilde{\omega}}t}\underline{x}_{+}^{-1}\underline{x}_{-}-\underline{x}_{-}e^{i\underline{\tilde{\omega}}t}\right)\underline{v}_{i}$
(40c) $\displaystyle\underline{x}_{g}\left(t\right)$
$\displaystyle=\left(\underline{x}_{+}\dfrac{e^{-i\underline{\tilde{\omega}}t}-1}{\underline{\tilde{\omega}}}\underline{x}_{+}^{-1}\underline{x}_{-}+\underline{x}_{-}\dfrac{e^{i\underline{\tilde{\omega}}t}-1}{\underline{\tilde{\omega}}}\right)\underline{v}_{i},$
(40d) where  $\displaystyle\underline{x}_{i}$ $\displaystyle=$
$\displaystyle\left(x_{+}+x_{-}\tilde{\omega}^{-1}x_{-}^{-1}x_{+}\tilde{\omega}\right)^{-1},$
(41a) $\displaystyle\underline{v}_{i}$ $\displaystyle=$
$\displaystyle\left(x_{-}\tilde{\omega}+x_{+}\tilde{\omega}x_{+}^{-1}x_{-}\right)^{-1}.$
(41b)

##  III 1st order LPS-N

Suppose that in the platform frame, one launches a particle at the time $t=0$
from the point $\mathbf{x}\left(0\right)$ at the velocity
$\mathbf{v}\left(0\right)$. After the interrogation time $T$, one measures the
position and velocity of the particle
$\left\\{\mathbf{x}\left(T\right),\mathbf{v}\left(T\right)\right\\}$ using the
local system of satellites and velocimeters and applying Eqs. (78, 81). The
navigation purpose is to restore the position, orientation, velocity, and
rotation rate of the platform,
$\left\\{\mathbf{X}\left(T\right),R\left(T\right),\mathbf{V}\left(T\right),\mathbf{\Omega}\left(T\right)\right\\},$
from these measurements.

After the launch, the particles are decoupled from the platform and move into
the Earth gravitational and inertial fields.. Its position in respect to the
Earth origin is given by

$\mathbf{X}_{E}\left(t\right)=\mathbf{X}_{c}+\mathbf{x}_{E}\left(t\right),$
(42)

where $\mathbf{x}_{E}\left(t\right)$ one obtains from Eq. (40a). If the origin
of the platform frame moves along the trajectory $\mathbf{X}\left(T\right)$,
then

$\mathbf{x}\left(T\right)=R\left(T\right)\left[\mathbf{x}_{E}\left(T\right)+\mathbf{X}_{c}-\mathbf{X}\left(T\right)\right].$
(43)

Since the derivative of the rotation matrix is determined by the rotation
rate,
$\dot{R}\left(T\right)=-\underline{\Omega}\left(T\right)R\left(T\right),$ one
obtains for the particle velocity

$\mathbf{v}\left(T\right)=-\mathbf{\Omega}\left(T\right)\times\mathbf{x}\left(T\right)+R\left(T\right)\left[\mathbf{v}_{E}\left(T\right)-\mathbf{V}\left(T\right)\right].$
(44)

The particle is detached from the platform only after launch, at $T>0$. If
$T=0$, one gets from Eqs. (43, 44)

$\displaystyle\mathbf{x}_{E}\left(0\right)$ $\displaystyle=$
$\displaystyle\mathbf{X}\left(0\right)-\mathbf{X}_{c}+R^{T}\left(0\right)\mathbf{x}\left(0\right),$
(45a) $\displaystyle\mathbf{v}_{E}\left(0\right)$ $\displaystyle=$
$\displaystyle\mathbf{V}\left(0\right)+R^{T}\left(0\right)\left[\mathbf{v}\left(0\right)+\mathbf{\Omega}\times\mathbf{x}\left(0\right)\right].$
(45b) Then one gets  $\displaystyle\mathbf{x}\left(T\right)$ $\displaystyle=$
$\displaystyle
R\left(T\right)\left\\{\underline{x}_{x}\left(T\right)\left[\mathbf{X}\left(0\right)-\mathbf{X}_{c}+R^{T}\left(0\right)\mathbf{x}\left(0\right)\right]\right.$
(46a)
$\displaystyle+\underline{x}_{v}\left(T\right)\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\left(\mathbf{v}\left(0\right)+\mathbf{\Omega}\left(0\right)\times\mathbf{x}\left(0\right)\right)\right]$
$\displaystyle\left.+\underline{x}_{g}\left(T\right)\mathbf{g}\left(\mathbf{X}_{c}\right)+\mathbf{X}_{c}-\mathbf{X}\left(T\right)\right\\};$
$\displaystyle\mathbf{v}\left(T\right)$ $\displaystyle=$
$\displaystyle-\mathbf{\Omega}\left(T\right)\times\mathbf{x}\left(T\right)+R\left(T\right)$
(46b)
$\displaystyle\times\left\\{\underline{\dot{x}}_{x}\left(T\right)\left[\mathbf{X}\left(0\right)-\mathbf{X}_{c}+R^{T}\left(0\right)\mathbf{x}\left(0\right)\right]\right.$
$\displaystyle+\underline{\dot{x}}_{v}\left(T\right)\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\left(\mathbf{v}\left(0\right)+\mathbf{\Omega}\left(0\right)\times\mathbf{x}\left(0\right)\right)\right]$
$\displaystyle\left.+\underline{\dot{x}}_{g}\left(T\right)\mathbf{g}\left(\mathbf{X}_{c}\right)-\mathbf{V}\left(T\right)\right\\}.$
The equation (46a) gives the exact trajectory of the particle in the platform
frame, which can, for example, be used to exact calculation of the AI phase,
instead of the approximate calculations in the articles c14.1 ; c14.2 ; c14.3
; c14.4 . Like in the previous studies c11 ; c13 , here we are going to use
Eq. (46) to restore the position and velocity of the platform. For them, one
has
$\displaystyle\mathbf{X}\left(T\right)=\mathbf{X}_{c}+\underline{x}_{x}\left(T\right)\left[\mathbf{X}\left(0\right)-\mathbf{X}_{c}+R^{T}\left(0\right)\mathbf{x}\left(0\right)\right]$
$\displaystyle+\underline{x}_{v}\left(T\right)\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\left(\mathbf{v}\left(0\right)+\mathbf{\Omega}\left(0\right)\times\mathbf{x}\left(0\right)\right)\right]$
$\displaystyle+\underline{x}_{g}\left(T\right)\mathbf{g}\left(\mathbf{X}_{c}\right)-R^{T}\left(T\right)\mathbf{x}\left(T\right),$
(47)
$\displaystyle\mathbf{V}\left(T\right)=\underline{\dot{x}}_{x}\left(T\right)\left[\mathbf{X}\left(0\right)-\mathbf{X}_{c}+R^{T}\left(0\right)\mathbf{x}\left(0\right)\right]$
$\displaystyle+\underline{\dot{x}}_{v}\left(t\right)\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\left(\mathbf{v}\left(0\right)+\mathbf{\Omega}\times\mathbf{x}\left(0\right)\right)\right]$
$\displaystyle+\underline{\dot{x}}_{g}\left(T\right)\mathbf{g}\left(\mathbf{X}_{c}\right)-R^{T}\left(T\right)\left[\mathbf{v}\left(T\right)+\mathbf{\Omega}\left(T\right)\times\mathbf{x}\left(T\right)\right].$
(48)

###  III.1 Rotation rate restore

One sees that in order to obtain the position and velocity of the platform, it
is necessary to initially restore the rotation rate
$\mathbf{\Omega}\left(T\right)$ and the rotation matrix $R\left(T\right)$. For
these restorations, one can use a differential technique. If from the same
point one launches 3 particles with velocities $\mathbf{v},$
$\mathbf{v}^{\prime}$ and $\mathbf{v}^{\prime\prime}$, then using Eq.. (46b)
to get the velocity differences at the time $T$

$\displaystyle\delta\mathbf{v}_{1}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{v}\left(T\right)-\mathbf{v}^{\prime}\left(T\right),$
(49a) $\displaystyle\delta\mathbf{v}_{2}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{v}\left(T\right)-\mathbf{v}^{\prime\prime}\left(T\right).$
(49b) one arrives to the following equations for
$\mathbf{\Omega}\left(T\right)$
$\mathbf{a}_{i}\times\mathbf{\Omega}\left(T\right)=\mathbf{b}_{i},$ (50)

where

$\displaystyle\mathbf{a}_{i}$ $\displaystyle\equiv$
$\displaystyle\delta\mathbf{x}_{i}\left(T\right),$ (51a)
$\displaystyle\mathbf{b}_{i}$ $\displaystyle\equiv$
$\displaystyle\delta\mathbf{v}_{i}\left(T\right)-R\left(T\right)\dot{x}_{v}\left(T\right)R^{T}\left(0\right)\delta\mathbf{v}_{i},$
(51b) $\displaystyle\delta\mathbf{v}_{i}$ $\displaystyle\equiv$
$\displaystyle\delta\mathbf{v}_{i}\left(0\right).$ (51c) One sees that one
velocity difference,for example $\delta\mathbf{v}_{1},$ is not sufficient, it
allows us to restore only the component of the rotation rate perpendicular to
$\mathbf{a}_{1}.$ That is why one starts 3 particles and considers 2 velocity
differences. Decomposing $\mathbf{\Omega}\left(T\right)$ as
$\mathbf{\Omega}\left(T\right)=\mathbf{\Omega}_{\perp}\left(T\right)+\Omega_{\parallel}\left(T\right)\mathbf{a}_{1},$
(52)

one gets from Eq. (50) for $i=1$

$\mathbf{\Omega}_{\perp}\left(T\right)=\dfrac{\mathbf{b}_{1}\times\mathbf{a}_{1}}{a_{1}^{2}}$
(53)

Substituting Eq.. (53) in Eq. (52) and then in Eq. (50) for i-2, one obtains
the equation for $\Omega_{\parallel}\left(T\right).$ After solving it, one
arrives to the following result

$\displaystyle\mathbf{\Omega}\left(t\right)=\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}\right),$
(54a)
$\displaystyle\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}\right)=\dfrac{\mathbf{b}_{1}\times\mathbf{a}_{1}}{a_{1}^{2}}$
$\displaystyle+\dfrac{\mathbf{a}_{1}\cdot\left\\{\left[\mathbf{b}_{2}a_{1}^{2}-\mathbf{b}_{1}\left(\mathbf{a}_{1}\cdot\mathbf{a}_{2}\right)\right]\times\mathbf{a}_{2}\right\\}}{\left(\mathbf{a}_{2}\times\mathbf{a}_{1}\right)^{2}a_{1}^{2}}\mathbf{a}_{1}.$
(54b)

###  III.2 Rotation matrix restore

Here, also, one should use the differential technique c11 ; c13 . For the same
3 particles with launch velocities $\mathbf{v},$ $\mathbf{v}^{\prime}$ and
$\mathbf{v}^{\prime\prime}$ instead of the velocities’ difference (49), one
now measures the positions’ differences

$\displaystyle\delta\mathbf{x}_{1}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{x}\left(T\right)-\mathbf{x}^{\prime}\left(T\right),$
(55a) $\displaystyle\delta\mathbf{x}_{2}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{x}\left(T\right)-\mathbf{x}^{\prime\prime}\left(T\right),$
(55b) which depend only on the rotation matrices,
$\displaystyle\delta\mathbf{x}_{i}\left(T\right)$ $\displaystyle=$
$\displaystyle r\mathbf{a}_{i},$ (56a) $\displaystyle r$ $\displaystyle=$
$\displaystyle R\left(T\right)R^{T}\left(0\right),$ (56b)
$\displaystyle\mathbf{a}_{i}$ $\displaystyle=$ $\displaystyle
R\left(0\right)\underline{x}_{v}\left(T\right)R^{T}\left(0\right)\delta\mathbf{v}_{i},$
(56c) where $\delta\mathbf{v}_{i}$ is defined in Eq. (51c). The matrix $r$ in
the Rodriguez representation c16.1  is given by
$r=\cos\delta+\dfrac{\left(1-\cos\delta\right)}{\delta^{2}}\mathbf{\delta\delta}^{T}-\dfrac{\sin\delta}{\delta}\underline{\delta},$
(57)

where the vector $\mathbf{\delta}$ contains 3 parameters that define the
rotation. In the Ref. c13 , an iterative method was proposed for calculating
the angle of rotation $\mathbf{\delta}$ from the phase differences of $6$ AIs,
in which the atoms were launched at different speeds, and the wave vectors did
not lie in the same plane. Here we adapt this method for the case of LPS.

In the case of rotation at a small angle, $\delta\ll 1$, when $r\approx
1-$$\delta$, and matrix $\delta$ defined as

$\underline{\delta}_{ik}=-\varepsilon_{ikm}\delta_{m},$ (58)

one gets that

$\mathbf{a}_{i}\times\mathbf{\delta}=\mathbf{b}_{i},$ (59)

where

$\mathbf{b}_{i}=\delta\mathbf{x}_{i}\left(T\right)-\mathbf{a}_{i},$ (60)

and $\left|b_{i}\right|\ll\left|a_{i}\right|$. This equation coincides with
Eq. (50), and hence,

$\mathbf{\delta}=\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}\right),$ (61)

where $\mathbf{f}$ is the vector function defined in Eq. (54b).

Let us now turn to the general case. Then, from Eq. (56a), one gets

$\sin\delta\mathbf{a}_{i}\times\mathbf{n}+\left(1-\cos\delta\right)\left[\mathbf{n}\left(\mathbf{na}_{i}\right)-\mathbf{a}_{i}\right]-\mathbf{b}_{i}=0,$
(62)

where $\mathbf{n}=\mathbf{\delta}/\delta$. One can introduce the parameter

$\sigma=\left|\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}\right)\right|$
(63)

and expand the solutions of Eq. (62) in the series

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}_{s=1}^{\infty}\delta^{\left(s\right)}\sigma^{s},$
(64a) $\displaystyle\mathbf{n}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}_{s=0}^{\infty}\mathbf{n}^{\left(s\right)}\sigma^{s}$
(64b) The lower order terms in these series are equal
$\displaystyle\delta^{\left(1\right)}$ $\displaystyle=$ $\displaystyle 1,$
(65a) $\displaystyle\mathbf{n}^{\left(0\right)}$ $\displaystyle=$
$\displaystyle\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}\right)/\sigma$
(65b) In Appendix B , we have derived the recurrence relations for the
coefficients
$\left\\{\mathbf{n}^{\left(s-1\right)},\delta^{\left(s\right)}\right\\}$.

##  IV 2nd order LPS-N

In this case, one restores the trajectory and orientation of the platform
$\left\\{\mathbf{X}\left(2T\right),R\left(2T\right)\right\\}$ from their
previous values at moments $t=0$ and $t=T$, using a second-order difference

$\mathbf{p}=\mathbf{x}\left(2T\right)-2x\left(T\right)+\mathbf{x}\left(0\right),$
(66)

in which the particle positions in the platform frame,
$\mathbf{x}\left(T\right)$ and $x\left(2T\right)$ are measured using
satellites.

An important advantage of such navigation is a tremendous decrease of the
navigation sensitivity to the initial velocity of the platform. For successful
navigation, it is sufficient to put c11 ; c13

$\mathbf{V}\left(0\right)=\dfrac{\mathbf{X}\left(T\right)-\mathbf{X}\left(-T\right)}{2T},$
(67)

where, for a given loop of the navigation cycle, $t\in\left[0,2T\right]$, the
positions of the platform $\mathbf{X}\left(-T\right)$ and
$\mathbf{X}\left(T\right)$ are restored in the previous loops.

If all the particles are launched from the center of the platform
$\left[\mathbf{x}\left(0\right)=0\right],$ then the lever=arm term,
$\mathbf{\Omega}\left(0\right)\times\mathbf{x}\left(0\right),$ disappears, and
simultaneously the need to restore the angular velocity of the platform
rotation $\mathbf{\Omega}\left(0\right)$ disappears also. Let’s also assume
for simplicity that $\mathbf{X}_{c}=\mathbf{X}\left(0\right)$. Then

$\displaystyle\mathbf{p}$ $\displaystyle=$
$\displaystyle\left[R\left(2T\right)\underline{x}_{v}\left(2T\right)-2R\left(T\right)\underline{x}_{v}\left(T\right)\right]$
(68)
$\displaystyle\times\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\mathbf{v}\left(0\right)\right]$
$\displaystyle+\left[R\left(2T\right)\underline{x}_{g}\left(2T\right)-2R\left(T\right)\underline{x}_{g}\left(T\right)\right]\mathbf{g}\left(\mathbf{X}\left(0\right)\right)$
$\displaystyle-R\left(2T\right)\left[\mathbf{X}\left(2T\right)-\mathbf{X}\left(0\right)\right]+2R\left(T\right)\left[\mathbf{X}\left(T\right)-\mathbf{X}\left(0\right)\right].$

In the absence of rotation $\left[R\left(t\right)=I\right]$ and for a uniform,
non-rotating gravitational field

$\underline{x}_{v}\left(t\right)=t;$ (69)

the first term in Eq. (68) disappears and the difference $\mathbf{p}$ does not
depend on the initial velocity of the platform. The simulation showed that
despite the weak inhomogeneity of the field, the presence of inertial forces
due to the rotation of the Earth and the small rotations of the gimbal, the
first term remains so small that the error in restoring the velocity (67) has
a much smaller impact on the accuracy of navigation than other factors, and
one may not take it into account. Then the coordinate of the platform
$\mathbf{X}\left(2T\right)$ is restored as

$\displaystyle\mathbf{X}\left(2T\right)$ $\displaystyle=$ $\displaystyle
R^{T}\left(2T\right)\left\\{-\mathbf{p}+\left[R\left(2T\right)\underline{x}_{v}\left(2T\right)-2R\left(T\right)\underline{x}_{v}\left(T\right)\right]\right.$
(70)
$\displaystyle\times\left[\mathbf{V}\left(0\right)+R^{T}\left(0\right)\mathbf{v}\left(0\right)\right]$
$\displaystyle+\left[R\left(2T\right)\underline{x}_{g}\left(2T\right)-2R\left(T\right)\underline{x}_{g}\left(T\right)\right]\mathbf{g}\left(\mathbf{X}\left(0\right)\right)$
$\displaystyle\left.+2R\left(T\right)\mathbf{X}\left(T\right)+\left[R\left(2T\right)-2R\left(T\right)\right]\mathbf{X}\left(0\right)\right\\}.$

If one neglects the rotation of the gimbal, assuming
$R\left(2T\right)=R\left(T\right)=R\left(0\right)$, then no further steps are
required. Otherwise, like in Sec. III.2, one launches simultaneously from the
origin 3 particles with velocities $\mathbf{v}\left(0\right)=\mathbf{v},$
$\mathbf{v}^{\prime},$ and $\mathbf{v}^{\prime\prime}$ and measures the double
differences

$\displaystyle\delta\mathbf{p}_{1}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{p}\left(T\right)-\mathbf{p}^{\prime}\left(T\right),$
(71a) $\displaystyle\delta\mathbf{p}_{2}\left(T\right)$ $\displaystyle=$
$\displaystyle\mathbf{p}\left(T\right)-\mathbf{p}^{\prime\prime}\left(T\right),$
(71b) which lead to the equations for the rotation matrix $R\left(2T\right)$.
In this case, one comes to the equation
$\delta\mathbf{p}_{i}+2R\left(T\right)\underline{x}_{v}\left(T\right)R^{T}\left(0\right)\delta\mathbf{v}_{i}=ra_{i}$
(72)

where

$\displaystyle r$ $\displaystyle=$ $\displaystyle
R\left(2T\right)R^{T}\left(0\right),$ (73a) $\displaystyle\mathbf{a}_{i}$
$\displaystyle=$ $\displaystyle
R\left(0\right)\underline{x}_{v}\left(2T\right)R^{T}\left(0\right)\delta\mathbf{v}_{i},$
(73b) where $\delta\mathbf{v}_{i}$ is defined in Eq. (51c). After that one
follows the method of solving Eq. (56a), described above and in appendix B,
with the only difference that now $\mathbf{a}_{i}$ is given by Eq. (73b) and
$\mathbf{b}_{i}=\delta\mathbf{p}_{i}+2R\left(T\right)\underline{x}_{v}\left(T\right)R^{T}\left(0\right)\delta\mathbf{v}_{i}-\mathbf{a}_{i}.$
(74)

##  V Simulation results

We perform testing the LPS-N on an ensemble of 30 randomly generated platform
trajectories, shown in Fig. 4.

Figure 4: Platform trajectories. Platform moves along each trajectory during
the time $t_{n}=1$hour. The trajectories are chosen so that the platform
velocity $V\left(0\right)=V\left(t_{n}\right)=0\ $and it reaches a maximum at
$t=t_{n}/2$. Only those trajectories were selected for which
$V\left(t_{n}/2\right)\in\left[200\text{m/s,}240\text{m/s}\right]$. (a)
Trajectories latitudes and longitudes. (b) Time dependences of the platform
trajectories’ altitudes. (The trajectories were randomly generated in the
western part of the Earth northern hemisphere for the sole purpose that they
were visible in Fig. 4a)

Regarding the orientation, we assume that the North-East-Down (NED) frame has
been installed initially on the platform, i.e. the initial gimbal rotation
matrix is given by

$R\left(0\right)=\left(\begin{array}[]{ccc}-\sin\phi\cos\lambda&-\sin\phi\sin\lambda&\cos\phi\\\
-\sin\lambda\cos\lambda&\cos\lambda&0\\\
-\cos\phi\cos\lambda&-\cos\phi\sin\lambda&-\sin\phi\end{array}\right),$ (75)

where $\phi$ and $\lambda$ are initial latitudes and longitude of the
platform. Later on gimbal commits only random noisy rotations at small angles
of the order of gimbal stabilization level $\delta_{g}$, in the frequency band
$f_{g}=100$Hz.

Our navigation code consists of the following points

  1. 1.

One calculates the exact values of the particle’s coordinate and velocity
using the Eqs. (46, 40b \- 40d) and exact expressions for the trajectory and
rotation matrix of the platform and the Earth’s gravitational field (see
appendix C below).

  2. 2.

Coordinates and velocities provide us accurate ”readings” of satellites and
velocimeters. One exploits the pseudo-random normal number generators c12
with SD $\delta_{x}$. and $\delta_{v}$ to make corrections to the sensors’
readings

  3. 3.

The sensor data are the input parameters for restoring the coordinates and
velocities in the platform frame using Eqs. (78, 81).

  4. 4.

The orientation of the gimbal is either assumed unchanged (for methods 1.1,
1.2, 2.1), or restored (for methods 1.3 and 2.2). In the last case, one solves
Eq. (62) using the iterative method (see appendix B).

  5. 5.

If necessary, one restores the instantaneous rotation rate of the gimbal
$\mathbf{\Omega}\left(T\right)$ [according to the formula (54a)].

  6. 6.

Finally, after restoring the velocity of the platform [see Eqs. (48) or (67)],
one obtains the desired coordinates of the platform from the Eqs. (47) or (70)

With the exception of the first, in the each loop of the navigation cycle, one
uses the initial position, speed, orientation, and rotation speed of the
platform restored in the previous loop. As a result, the navigation error
$\Delta X$ accumulates and increases on 7 orders of magnitude, from $1{\mu}$
to $10$m at the end of the cycle.

The accuracy of navigation increases with more accurate measurement of the
position of the particles in the platform frame. To achieve it, we propose
using not one set of satellites and velocimeters, but ensembles of these
sensors. Specifically, during the simulation, we assumed that there are
$n_{x}$ satellites located in the platform on each of the semi-axes. Then, it
is obvious that the ensemble of measured coordinates consists of

$N_{x}=n_{x}^{6}$ (76)

elements. At the same time, on each of the axes, we have arranged
$n_{v}=n_{x}^{2}$ velocimeters so that the ensemble of velocity measurements
consists of the same number of elements

$N_{v}=n_{v}^{3}=N_{x}.$ (77)

The simulation showed that the navigation error on the one side is not
sensitive to the position of satellites and velocimeters, and on the other
side depends on the values of particles’ launch velocities. We have not been
able to find any qualitative or quantitative criteria to minimize this
dependence. We emphasize that we are talking about a 9-dimensional velocity
space $\left\\{\mathbf{v,v}^{\prime}\mathbf{,v}^{\prime\prime}\right\\}$.
Therefore, we randomly selected the velocity components in a given range,
gradually increasing the size of this range until, at $n_{x}=1$, the
navigation accuracy drops to the value of $\Delta X\left(t_{n}\right)\sim
100$m. We expected that with an increase in the number of satellites along
each of the semi-axes $n_{x}$ and a simultaneous increase in the number of
velocimeters $n_{v}=n_{x}^{2}$, the error of the navigation will drop to the
desired level of $\Delta X\left(t_{n}\right)\sim 10$m. The simulation
confirmed this expectation (see the red curves in the fourth column in Fig.
6). The values of launching velocities chosen for the simulation are given in
the table 1.

We assume that a navigation cuboid with dimensions
$\left\\{L_{N},L_{E},L_{D}\right\\}$ is installed in the (NED)-coordinate
system of the platform, which is shielded from all fields and in which a
vacuum is maintained sufficient for one to believe that particles move freely
only under the action of the Earth gravitational field. The minimal sizes of
the cuboid at which there will be no collision of particles with the cuboid
walls are also given in the table 1. Since in the 2nd order LPS-N one observes
the particle motion during the time $2T$, then the sizes of the cuboid are
larger than for the 1st order LPS-N at the given $T.$

Table 1: Launching velocities chosen for navigation simulations and corresponding sizes of the navigation cuboids. | M
---
e
t
h
o
d
| Launching
---
velocities
[m/s]
| Minimal
---
cuboid
size
[m]
1.1 | $\begin{array}[]{c}v_{N}=-0.40\\\ v_{E}=-0.40\\\ v_{D}=-0.99\end{array}$ | $\begin{array}[]{c}L_{N}=0.12\\\ L_{E}=0.12\\\ L_{D}=0.20\end{array}$
1.2 | $\begin{array}[]{ccc}v_{N}=-0.40&v_{N}^{\prime}=-0.40&v_{N}^{\prime\prime}=0.40\\\ v_{E}=-0.40&v_{E}^{\prime}=0.40&v_{E}^{\prime\prime}=-0.40\\\ v_{D}=-1.55&v_{D}^{\prime}=-3.04&v_{D}^{\prime\prime}=-1.79\end{array}$ | $\begin{array}[]{c}L_{N}=0.12\\\ L_{E}=0.12\\\ L_{D}=0.49\end{array}$
1.3 | $\begin{array}[]{ccc}v_{N}=2.06&v_{N}^{\prime}=-1.96&v_{N}^{\prime\prime}=0.13\\\ v_{E}=-0.29&v_{E}^{\prime}=0.76&v_{E}^{\prime\prime}=0.67\\\ v_{D}=-4.54&v_{D}^{\prime}=-2.00&v_{D}^{\prime\prime}=1.15\end{array}$ | $\begin{array}[]{c}L_{N}=1.2\\\ L_{E}=0.3\\\ L_{D}=1.7\end{array}$
2.1 | $\begin{array}[]{ccc}v_{N}=1.77&v_{N}=1.77&v_{N}=-1.77\\\ v_{E}=1.77&v_{E}^{\prime}=-1.77&v_{E}^{\prime}=1.77\\\ v_{D}=2.37&v_{D}^{\prime}=2.05&v_{D}^{\prime\prime}=1.51\end{array}$ | $\begin{array}[]{c}L_{N}=2.1\\\ L_{E}=2.1\\\ L_{D}=2.4\end{array}$
2.2 | $\begin{array}[]{ccc}v_{N}=2.94&v_{N}^{\prime}=2.94&v_{N}^{\prime\prime}=-2.94\\\ v_{E}=2.94&v_{E}^{\prime}=-2.94&v_{E}^{\prime\prime}=2.94\\\ v_{D}=0.99&v_{D}^{\prime}=-5.61&v_{D}^{\prime\prime}=-2.97\end{array}$ | $\begin{array}[]{c}L_{N}=3.6\\\ L_{E}=3.6\\\ L_{D}=3.95\end{array}$

To illustrate the LPS-N, the build up in time of the navigation errors when
the platform moves along the trajectory 20 (see Fig. 4) is shown in Fig. 5.
Although the various plots on this and next figures may be difficult to read
at standard magnification, they can be read easily using the zoom feature when
read online in PDF format. The navigation method 1.3 is used here. In this
method, the largest number of platform parameters are restored, the
orientation of the gimbal, then sequentially the rotation rate of the gimbal,
the velocity and, finally, the position of the platform. At the same time,
only an error in restoring orientation is caused with the accuracy of
measuring the difference in the coordinates of particles launched at different
velocities. Errors in restoring the subsequent parameter are related both to
the quality of satellites and velocimeters, and to the accuracy of restoring
previous parameters. It is for this reason that the error in Fig.
5$b^{\prime}$ is 5 orders of magnitude larger than the error in Fig. 2.

Figure 5: The time dependence of the parameters of the platform motion and
rotation (plots $i$) and the errors of their restoring.(plots $i^{\prime}$),
where $i=a,...,e$. Trajectory 20 (see fig. 4). LPS-N method 1.3. Plots
$\left\\{a,a^{\prime}\right\\}$ for the altitude of the platform origine
trajectory, plots $\left\\{b,b^{\prime}\right\\}$ for the latitude and
longitude of this trajectory, $\left\\{c,c^{\prime}\right\\}$ for the platform
velocity $\mathbf{V}\left(t\right)$, $\left\\{d,d^{\prime}\right\\}$ for the
matrix elements of the rotation matrix
$r\left(t\right)=R\left(t\right)R^{T}\left(0\right)$,
$\left\\{e,e^{\prime}\right\\}$ for the platform rotation rate
$\mathbf{\Omega}\left(t\right).$

The simulation results are summarized in the table in Fig. 6.

Figure 6: Five LPS-N methods. Results of simulations. For the same 30
trajectories as in FIG. 1 one generated a navigation process. This was done
for 11 sets of pseudo-random errors of satellites and velocimeters. As a
result, one has ensemble of 330 $\Delta X\left(t\right)-$dependences. The
average over this ensemble errors for $t_{n}=1$hour of the navigation are
shown by red dots in the 4th column. In the same plots, the blue, green and
black dots show the maximum and minimum values of the error and SD for the
entire ensemble, respectively. The $\Delta X\left(t\right)-$dependences for 30
platform trajectories with a given set of sensors’ errors are shown in the 2nd
and 3rd columns. They contain plots for those sets of sensor errors at which
one of the $\Delta X\left(t\right)$ curves reaches the minimum (best) or
maximum (worst) value of $\Delta X\left(t_{n}\right)$.

The five rows of this table correspond to three 1st order LPS-N methods (1.1,
1.2, and 1.3) and two 2nd order LPS-N methods (2.1 and 2.2). One sees that in
method 1.1, when it is assumed that the quality of the gimbal is sufficient to
neglect the rotation at all, in order to achieve the ultimate goal of
navigation (8), it is necessary to use a gimbal with a rms of
$\delta_{g}=10^{-7}rad$ and velocimeters with an accuracy of
$\delta_{v}=20\mu/s$. The requirements can be mitigated for the gimbal by an
order of magnitude, and for the velocimeters by a factor of 2.5 if one uses
method 1.2, when only the angular rate of the gimbal rotation is restored.
Moreover, one can weaken the requirement for a gimbal by an order of magnitude
in the method 1.3, when both the rate of rotation and the orientation of the
gimbal are restored, but satellites with an rms error of $\delta_{x}=100$nm
are needed. We have verified that with another random choice of trajectories
and orientations, the navigation error $\Delta X\left(t\right)$ remains the
same in the order of its magnitude.

If the 2nd order difference (5b) is used, then, due to its low sensitivity to
the platform velocity, one may not use the velocimeters at all. For both 2nd
order LPS-N methods, 2.1 and 2.2, we verified that even the use of ideal
velocimeters, with zero error of velocity measurement, has no effect on the
navigation error $\Delta X\left(t\right)$. In addition, the lever-arm term
disappears, when one launches all three particles from the same point [see Eq.
(68)], and therefore the platform rotation rate may not be restored. Finally,
with the 2nd order LPS-N , the requirements for the quality of gimbal are
weakened by an order of magnitude (compare the 4th and 1st rows in the table
in Fig. 6).

In conclusion, the simulation showed that the LPS can provide navigation
accuracy in 1 hour of the order of 10m, if a gimbal with rms
$\delta_{g}=10^{-6}$rad and satellites, distance sensors, with an error
$\delta_{x}=1\mu$ are used.

However, the reverse formulation of the question is also possible. If the
accuracy of the satellites and velocimeters, the rms of gimbal, and the
navigation box sizes are given, then our code will allow us to determine the
expected accuracy of LPS-navigation $\Delta X$ for a given time $t_{n}.$

###### Acknowledgements.

The author is grateful to Drs. K. Tintschurin, B. Young, I. Teper, and A. Zorn
for the fruitful discussions.

##  Appendix A Particle’s position and velocity measurements in the platform
frame

###  A.1 LPS - measurement of the position

Let’s assume that there are 3 satellites installed on the platform at the
points $\mathbf{s}_{i}=\left(x_{i},y_{i},z_{i}\right)^{T}$, which measure the
distances to a given point $\mathbf{x}=\left(x,y,z\right)^{T},$
$r_{i}=\left|\mathbf{x}-\mathbf{s}_{i}\right|$. Then for the point coordinates
one has c17-1

$\displaystyle z=\left(-b\pm\sqrt{b^{2}-ac}\right)/a,$ (78a) $\displaystyle
x=b_{0x}+b_{1x}z,$ (78b) $\displaystyle y=b_{0y}+b_{1y}z,$ (78c)
$\displaystyle b=b_{0x}b_{1x}+b_{0y}b_{1y}-b_{1x}x_{1}-b_{1y}y_{1}-z_{1},$
(78d) $\displaystyle a=1+b_{1x}^{2}+b_{1y}^{2},$ (78e) $\displaystyle
c=\rho_{1}^{2}-r_{1}^{2}+b_{0x}^{2}+b_{0y}^{2}-2\left(b_{0x}x_{1}+b_{0y}y_{1}\right),$
(78f) $\displaystyle b_{0x}=\left(a_{2}y_{31}-a_{3}y_{21}\right)/\Delta,$
(78g) $\displaystyle b_{1x}=\left(z_{31}y_{21}-z_{21}y_{31}\right)/\Delta,$
(78h) $\displaystyle b_{0y}=\left(x_{21}a_{3}-x_{31}a_{2}\right)/\Delta,$
(78i) $\displaystyle b_{1y}=\left(x_{31}z_{21}-x_{21}z_{31}\right)/\Delta,$
(78j) $\displaystyle\Delta=x_{21}y_{31}-y_{21}x_{31},$ (78k)
$\displaystyle\rho_{i}=\left|\mathbf{s}_{i}\right|,$ (78l) $\displaystyle
a_{i}=\dfrac{1}{2}\left(\rho_{i}^{2}-\rho_{1}^{2}-r_{i}^{2}+r_{1}^{2}\right),$
(78m)
$\displaystyle\left(x_{ik},y_{ik},z_{ik}\right)=\left(\mathbf{s}_{i}-\mathbf{s}_{k}\right)^{T}.$
(78n)

###  A.2 Velocity measurement

For the 1st order LPS-N, one must measure the velocity $\mathbf{v}$ of a
particle in the platform frame at a given point $\mathbf{x}$. If 3
velocimeters, located at points $\mathbf{s}_{vi}$, measure the projection of
the particle velocity

$v^{\left(i\right)}=\mathbf{n}_{i}\cdot\mathbf{v},$ (79)

where

$\mathbf{n}_{i}=\dfrac{\mathbf{x}-\mathbf{s}_{vi}}{\left|\mathbf{x}-\mathbf{s}_{vi}\right|}.$
(80a) then for the particle velocity one gets
$\mathbf{v}=\underline{n}^{-1}\left(\begin{array}[]{c}v^{\left(1\right)}\\\
v^{\left(2\right)}\\\ v^{\left(3\right)}\end{array}\right),$ (81)

where the matrix $\underline{n}$ is given by

$\underline{n}=\left(\begin{array}[]{c}\mathbf{n}_{1}^{T}\\\
\mathbf{n}_{2}^{T}\\\ \mathbf{n}_{3}^{T}\end{array}\right).$ (82)

##  Appendix B Recurrent relations

To decompose

$\sin\delta=\mathop{\displaystyle\sum}_{s=0}^{\infty}\dfrac{\left(-1\right)^{s}\sigma^{2s+1}}{\left(2s+1\right)!}\left(\dfrac{\delta}{\sigma}\right)^{2s+1}$
(83)

one needs to get a series for a degree

$\left(\dfrac{\delta}{\sigma}\right)^{p}=p!\sum_{v=0}^{\infty}\delta_{v}^{\left(p\right)}\sigma^{v}$
(84)

To calculate the coefficient $\delta_{v}^{\left(p\right)},$ one truncates the
row (64) with the first $\left(\nu+1\right)$ terms

$\delta=\sigma\sum_{i=1}^{\nu+1}\delta^{\left(i\right)}\sigma^{i-1}.$ (85)

Then

$\displaystyle\left(\delta/\sigma\right)^{p}$ $\displaystyle=$ $\displaystyle
p!\mathop{\displaystyle\sum}_{r_{1}=0}^{p}\ldots\mathop{\displaystyle\sum}_{r_{\nu+1}=0}^{p}\delta_{p,\mathop{\displaystyle\sum}_{i=1}^{\nu+1}r_{i}}$
(86)
$\displaystyle\times\sigma^{\mathop{\displaystyle\sum}_{i=1}^{\nu+1}r_{i}\left(i-1\right)}\mathop{\displaystyle\prod}\limits_{i=1}^{\nu+1}\dfrac{\left[\delta^{\left(i\right)}\right]^{r_{i}}}{r_{i}!},$

where $\delta_{ik}$ is the Kronecker symbol, and so

$\displaystyle\delta_{v}^{\left(p\right)}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}_{r_{1}=0}^{p}\ldots\mathop{\displaystyle\sum}_{r_{\nu+1}=0}^{p}\delta_{p,\mathop{\displaystyle\sum}_{i=1}^{\nu+1}r_{i}}$
(87)
$\displaystyle\times\delta_{\nu,\mathop{\displaystyle\sum}_{i=1}^{\nu+1}r_{i}\left(i-1\right)}\mathop{\displaystyle\prod}\limits_{i=1}^{\nu+1}\dfrac{\left[\delta^{\left(i\right)}\right]^{r_{i}}}{r_{i}!}.$

Substituting the decomposition (84) in Eq. (83), one finds that

$\displaystyle\sin\delta$ $\displaystyle\equiv$
$\displaystyle\mathop{\displaystyle\sum}_{s=0}^{\infty}u_{s}\sigma^{s+1},$
(88a) $\displaystyle u_{s}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}_{s^{\prime}=0}^{\left[s/2\right]}\left(-1\right)^{s^{\prime}}\delta_{s-2s^{\prime}}^{\left(2s^{\prime}+1\right)},$
(88b) where $\left[a\right]$ is the Entire part of a. Similarly, one gets

$\displaystyle
1-\cos\delta=\mathop{\displaystyle\sum}_{s=0}^{\infty}c_{s}\sigma^{s+2},$
(89a) $\displaystyle
c_{s}=\left\\{\begin{array}[]{c}\mathop{\displaystyle\sum}_{s^{\prime}=0}^{\left[s/2\right]}\left(-1\right)^{s^{\prime}}\delta_{s-2s^{\prime}}^{\left(2s^{\prime}+2\right)},\text{
for }s\geq 0\\\ 0,\text{ for }s<0\end{array}\right..$ (89d) Multiplying the
power series in Eq. (62) one comes to the following equation for the
coefficients

$\mathop{\displaystyle\sum}_{s_{1}=0}^{s-1}\mathbf{w}_{is_{1}}-\mathbf{a}_{i}c_{s-2}+\mathop{\displaystyle\sum}_{s_{1}=0}^{s-2}c_{s_{1}}\mathop{\displaystyle\sum}_{s_{2}=0}^{s-s_{1}-2}\mathbf{n}^{\left(s_{2}\right)}\left(\mathbf{n}^{\left(s-s_{1}-2-s_{2}\right)}\mathbf{a}_{i}\right)=\dfrac{\mathbf{b}_{i}}{\sigma}\delta_{s1}.$
(90)

where

$\mathbf{w}_{is_{1}}=u_{s_{1}}\mathbf{a}_{i}\times\mathbf{n}^{\left(s-s_{1}-1\right)}$
(91)

It follows from Eqs. (88b, 87, 65a) that

$u_{0}=1$ (92)

and for $s=1$ one gets the equation

$\mathbf{a}_{i}\times\mathbf{n}^{\left(0\right)}=\dfrac{\mathbf{b}_{i}}{\sigma}$
(93)

the solution of which, taking into account the normalization condition,
$\left|\mathbf{n}^{\left(0\right)}\right|=1$, coincides with (65).

The coefficients of the our interest,
$\left\\{\mathbf{n}^{\left(s-1\right)},\delta^{\left(s\right)}\right\\}$ are
contained only in the terms $\mathbf{w}_{i0}$ and
$\mathbf{w}_{i\left(s-1\right)}$. For $\mathbf{w}_{i0}$ from Eq. (92) one has,

$\mathbf{w}_{i0}=\mathbf{a}_{i}\times\mathbf{n}^{\left(s-1\right)},$ (94)

and from Eqs. (88b, 93) one gets

$\mathbf{w}_{i\left(s-1\right)}=\left(\delta_{s-1}^{\left(1\right)}+q_{s}\right)\dfrac{\mathbf{b}_{i}}{\sigma}$
(95)

where

$q_{s}\equiv\mathop{\displaystyle\sum}_{s^{\prime}=1}^{\left[\left(s-1\right)/2\right]}\left(-1\right)^{s^{\prime}}\delta_{s-1-2s^{\prime}}^{\left(2s^{\prime}+1\right)}.$
(96)

From Eq. (87) for $p=1$ and $\nu=s-1$, one sees that
$\delta_{s-1}^{\left(1\right)}$ contains only one term at $\left\\{r_{1}\ldots
r_{s}\right\\}=\left\\{0\ldots 01\right\\}$, which is equal to

$\delta_{s-1}^{\left(1\right)}=\delta^{\left(s\right)}.$ (97)

Then, with $s>1$, one can rewrite Eq. (90) as

$\displaystyle\mathbf{a}_{i}\times\mathbf{n}^{\left(s-1\right)}=-\left(\delta^{\left(s\right)}+q_{s}\right)\dfrac{\mathbf{b}_{i}}{\sigma}+\mathbf{b}_{i}^{\left(s\right)},$
(98a)
$\displaystyle\mathbf{b}_{i}^{\left(s\right)}=\mathbf{a}_{i}c_{s-2}-\mathop{\displaystyle\sum}_{s_{1}=1}^{s-2}\mathbf{w}_{is_{1}}$
$\displaystyle-\mathop{\displaystyle\sum}_{s_{1}=0}^{s-2}c_{s_{1}}\mathop{\displaystyle\sum}_{s_{2}=0}^{s-s_{1}-2}\mathbf{n}^{\left(s_{2}\right)}\left(\mathbf{n}^{\left(s-s_{1}-2-s_{2}\right)}\cdot\mathbf{a}_{i}\right).$
(98b)

This is an equation of the same type as Eq. (50) for the speed of rotation and
Eq. (59) for the angle of rotation $\mathbf{\delta}$. By analogy with those
equations, the solution of Eq. (19a) is

$\mathbf{n}^{\left(s-1\right)}=-\left(\delta^{\left(s\right)}+q_{s}\right)\mathbf{n}^{\left(0\right)}+\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(s\right)}\right),$
(99)

where the function $\mathbf{f}$ is given by Eq. (54b). Let us consider the
normalization condition $\mathbf{n}^{2}=1.$ It means that

$\displaystyle
2\left[n^{\left(0\right)}\cdot\mathbf{n}^{\left(s-1\right)}+\mathop{\displaystyle\sum}_{s_{1}=1}^{\left[s/2\right]-1}\mathbf{n}^{s_{1}}\mathbf{n}^{s-s_{1}-1}\right]$
(100) $\displaystyle=$
$\displaystyle-\delta_{s,2\left[s/2\right]+1}\left[\mathbf{n}^{\left(\left(s-1\right)/2\right)}\right]^{2}$

Substituting here. Eq. (99) one arrives to the expression for the coefficient
$\delta^{\left(s\right)},$

$\displaystyle\delta^{\left(s\right)}=-q_{s}+\mathbf{n}^{\left(0\right)}\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(s\right)}\right)+$
$\displaystyle\mathop{\displaystyle\sum}_{s_{1}=1}^{\left[s/2\right]-1}\mathbf{n}^{s_{1}}\mathbf{n}^{s-s_{1}-1}+\dfrac{1}{2}\delta_{s,2\left[s/2\right]+1}\left[\mathbf{n}^{\left(\left(s-1\right)/2\right)}\right]^{2}.$
(101)

The expressions (99, 101) are the desired recurrent relations.

Consider them for the lowest orders. If $s=2,$ then from the Eqs..(89a, 96) it
is obvious that $c_{0}=\dfrac{1}{2}$ and $q_{2}=0,$ which means that

$\mathbf{b}_{i}^{\left(2\right)}=\dfrac{1}{2}\left[\mathbf{a}_{i}-\mathbf{n}^{\left(0\right)}\left(\mathbf{n}^{\left(0\right)}\cdot\mathbf{a}_{i}\right)\right].$
(102)

and, therefore,

$\displaystyle\delta^{\left(2\right)}$ $\displaystyle=$
$\displaystyle\mathbf{n}^{\left(0\right)}\cdot\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(2\right)}\right),$
(103a) $\displaystyle\mathbf{n}^{\left(1\right)}$ $\displaystyle=$
$\displaystyle-\delta^{\left(2\right)}\mathbf{n}^{\left(0\right)}+\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(2\right)}\right).$
(103b)

For $s=3$, one consequently finds that

$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle
u_{1}=\delta^{\left(2\right)},$ (104a) $\displaystyle\mathbf{w}_{i1}$
$\displaystyle=$
$\displaystyle\delta^{\left(2\right)}\mathbf{a}_{i}\times\mathbf{n}^{\left(1\right)},$
(104b) after that  $\displaystyle\mathbf{b}_{i}^{\left(3\right)}$
$\displaystyle=$
$\displaystyle\delta^{\left(2\right)}\left[\mathbf{a}_{i}-\mathbf{a}_{i}\times\mathbf{n}^{\left(1\right)}-\mathbf{n}^{\left(0\right)}\left(\mathbf{n}^{\left(0\right)}\mathbf{a}_{i}\right)\right]$
(105)
$\displaystyle-\dfrac{1}{2}\left(\mathbf{n}^{\left(0\right)}\left(\mathbf{n}^{\left(1\right)}\cdot\mathbf{a}_{i}\right)+\mathbf{n}^{\left(1\right)}\left(\mathbf{n}^{\left(0\right)}\cdot\mathbf{a}_{i}\right)\right).$

And since $q_{3}=-\dfrac{1}{6}$, then

$\delta^{\left(3\right)}=\dfrac{1}{6}+\mathbf{n}^{\left(0\right)}\cdot\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(3\right)}\right)+\dfrac{1}{2}\mathbf{n}^{\left(1\right)2},$
(106)

and, finally,

$\mathbf{n}^{\left(2\right)}=\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(3\right)}\right)-\mathbf{n}^{\left(0\right)}\left[\mathbf{n}^{\left(0\right)}\cdot\mathbf{f}\left(\mathbf{a}_{i},\mathbf{b}_{i}^{\left(3\right)}\right)+\dfrac{1}{2}\mathbf{n}^{\left(1\right)2}\right].$
(107)

##  Appendix C Earth gravitational field

We consider navigation in the gravitational field of a rotating Earth. This
field consists of normal and anomaly parts

$U\left(\mathbf{X}\right)=U_{n}\left(\mathbf{X}\right)+U_{a}\left(\mathbf{X}\right)$
(108)

For the potential of the normal component, one has c17

$\displaystyle U_{n}\left(\mathbf{X}\right)$
$\displaystyle=-\dfrac{GM}{X}\left(1-\sum_{m=1}^{\infty}J_{2m}\left(\dfrac{a_{E}}{X}\right)^{2m}P_{2m}\left(\cos\theta\right)\right),$
(109a) $\displaystyle J_{2m}$
$\displaystyle=\left(-1\right)^{m+1}\dfrac{3e^{2m}}{\left(2m+1\right)\left(2m+3\right)}\left(1-m\left(1-\dfrac{5J_{2}}{e^{2}}\right)\right),$
(109b) where $GM=3.986004418\cdot 10^{14}$m${}^{\text{3}}$s${}^{\text{-2}}$ is
the geocentric gravitational constant, $a=6.37813659\cdot 10^{6}$m is the
semimajor axis of the geoid, $e=\sqrt{f\left(2-f\right)}$ is the first
eccentricity, $f$ is polar flattening $\left(1/f=298.25765\right)$,
$J_{2}=1.0826267\cdot 10^{-3}$ is a dynamic form factor of the Earth c18 ,
$P_{n}\left(x\right)$ is the Legendre polynomial, $\theta$ is the polar angle
of the vector $\mathbf{X}$. From here one derives the following expressions
for the normal components of the gravitational acceleration and the gravity-
gradient tensor of the Earth
$\displaystyle\mathbf{g}_{nE}\left(\mathbf{X}\right)$ $\displaystyle=$
$\displaystyle\mathbf{g}=-\dfrac{GM}{x^{2}}\left(b\mathbf{n}+b_{1}\mathbf{z}\right);,$
(110a) $\displaystyle\underline{\Gamma}_{n}\left(\mathbf{X}\right)$
$\displaystyle=$
$\displaystyle\dfrac{GM}{x^{3}}\left(-b\delta_{ik}+b_{2}n_{i}n_{k}+b_{3}\left(n_{i}z_{k}+z_{i}n_{k}\right)\right.$
(110b) $\displaystyle\left.+b_{4}z_{i}z_{k}\right),$ where  $\displaystyle
b=1-\dfrac{1}{\sin^{2}\theta}\mathop{\displaystyle\sum}_{m=1}^{\infty}J_{2m}\left(\dfrac{a}{x}\right)^{2m}\left[P_{2m}\left(\cos\theta\right)\left(2m+1\right.\right.$
$\displaystyle\left.\left.-\left(4m+1\right)\cos^{2}\theta\right)+2m\cos\theta
P_{2m-1}\left(\cos\theta\right)\right],$ (111a) $\displaystyle
b_{1}=\dfrac{1}{\sin^{2}\theta}\mathop{\displaystyle\sum}_{m=1}^{\infty}J_{2m}\left(\dfrac{a}{x}\right)^{2m}2m\left[P_{2m-1}\left(\cos\theta\right)\right.$
$\displaystyle\left.-\cos\theta P_{2m}\left(\cos\theta\right)\right],$ (111b)
$\displaystyle
b_{2}=3-\dfrac{1}{\sin^{4}\theta}\mathop{\displaystyle\sum}_{m=1}^{\infty}J_{2m}\left(\dfrac{a}{x}\right)^{2m}\left\\{P_{2m}\left(\cos\theta\right)\left[\left(3+2m\right)\right.\right.$
$\displaystyle\left.\times\left(2m+1\right)-\cos^{2}\theta\left(20m^{2}+28m+6\right)\right]$
$\displaystyle\left.+\cos^{4}\theta\left(16m^{2}+16m+3\right)\right]+2m\cos\theta$
$\displaystyle\left.P_{2m-1}\left(\cos\theta\right)\left[5+4m-\left(3+4m\right)\cos^{2}\theta\right]\right\\},$
(111c) $\displaystyle
b_{3}=\dfrac{1}{\sin^{4}\theta}\mathop{\displaystyle\sum}_{m=1}^{\infty}J_{2m}\left(\dfrac{a}{x}\right)^{2m}2m\left\\{2P_{2m-1}\left(\cos\theta\right)\left(m+1\right.\right.$
$\displaystyle\left.\left.-m\cos^{2}\theta\right)-\cos\theta
P_{2m}\left(\cos\theta\right)\left[4m+3-\left(4m+1\right)\cos^{2}\theta\right]\right\\},$
(111d) $\displaystyle
b_{4}=\dfrac{1}{\sin^{4}\theta}\mathop{\displaystyle\sum}_{m=1}^{\infty}J_{2m}\left(\dfrac{a}{x}\right)^{2m}2m$
$\displaystyle\times\left\\{\left[2m+1-\left(2m-1\right)\cos^{2}\theta\right]P_{2m}\left(\cos\theta\right)\right.$
$\displaystyle\left.-2\cos\theta P_{2m-1}\left(\cos\theta\right)\right\\}$
(111e) $\mathbf{z}$ is a unit vector along the Z-axis.

For the anomaly part of the potential $U_{a}\left(\mathbf{X}\right),$ we could
not find an analytical expression. Nevertheless, it is known c19  that the rms
of anomalous acceleration

$g_{a}\sim 30\text{mGal}$ (112)

and from Plate 5 in Ref. c20  one can conclude that the rms of the anomaly
part of the gravity-gradient tensor $\Gamma_{a}\sim 100$E. It follows that the
typical size of the anomaly field oscillations in space is

$L_{a}\sim g_{a}/\left|\Gamma_{a}\right|\sim 3\text{km.}$ (113)

Unfortunately, the known analytical c21  or numerical c21.1 ; c22  approaches
to the calculation of the gravitational field in the space around the Earth,
including its anomalous part, cannot be used in our calculations, since they
assume a known field on the Earth’s surface. Therefore, for our code, we
assumed that the anomalous potential consists of the sum of sinusoidal terms

$U_{a}\left(\mathbf{X}\right)=\sum_{j=1}^{n_{a}}h_{aj}\cos\left(\mathbf{k}_{aj}\cdot\mathbf{X}\right),$
(114)

in which the wave vectors were generated by the formula

$k_{aji}=\pm 2\pi/L_{a}(1+2r),$ (115)

where $r$ is a pseudo-random number uniformly distributed over the interval
$\left[0,1\right]$, and for amplitudes $h_{aj}$, a generator of normal pseudo-
random numbers was used, after which they were normalized by the condition

$g_{a}^{2}=\dfrac{1}{2}\sum_{j=1}^{n_{a}}h_{aj}^{2}k_{aj}^{2}.$ (116)

Of course, Eq. (114) has nothing common with reality [it follows, for example,
that particles move in a hypothetical medium with a density of
$\rho=-4\pi\triangle U_{a}\left(\vec{r}\right)/G\not=0,$ which can even be
negative. Despite this, we accept it, since at $n_{a}\gg 1$, the anomalous
field has a rms (112) and a spatial period (113). In addition, the
representation of the anomalous potential in the form of sums of harmonic
terms can qualitatively correspond to the topographic undulations observed by
the GRACE and GOCE satellites.

## References

  * (1) B Ya Dubetskii, A P. Kazantsev, V P. Chebotayev, V P. Yakovlev, Interference of atoms in separated optical fields, Pis’ma Zh. Eksp. Teor. Fiz. 39, 531 (1984) [JETP Lett. 39, 649 (1984)].
  * (2) G M Tino and M A Kasevich(ed) 2014 Atom interferometry Proc. of the Int. School of Physics ‘Enrico Fermi’ vol 188 (IOS Press)
  * (3) B. Canuel et al, Technologies for the ELGAR large scale atom interferometer array, arXiv:2007.04014v1 [physics.atom-ph]
  * (4) Ming-Sheng Zhan et al, ZAIGA: Zhaoshan long-baseline atom interferometer gravitation antenna, Int. J. Mod. Phys. D 29, 194005 (2020).
  * (5) L. Badurina et al, AION: an atom interferometer observatory and network, J. Cosmol. Astropart. Phys. JCAP05(2020)011
  * (6) B. Battelier et al, Exploring the Foundations of the Universe with Space Tests of the Equivalence Principle, arXiv:1908.11785v3 [physics.space-ph]
  * (7) Y. A. El-Neaj et al, AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space,EPJ Quantum Technology  7, 6 (2020).
  * (8) G M. Tino et al, SAGE: A proposal for a space atomic gravity explorer, Eur. Phys. J. D. 73, 228 (2019).
  * (9) M. Abe et al, Matter-wave Atomic Gradiometer Interferometric Sensor (MAGIS-100), arXiv:2104.02835 [physics.atom-ph].
  * (10) M. Kasevich, S. Chu, Atomic interferometry using stimulated Raman transitions, Phys. Rev. Lett. 67, 181 (1991)..
  * (11) P. Asenbaum, C. Overstreet, M. Kim, J. Curti, and M. A. Kasevich, Atom-Interferometric Test of the Equivalence Principle at the $10^{-12}$ Level, Phys. Rev. Lett. 125, 191101 (2020).
  * (12) B. Canuel, F. Leduc, D. Holleville, A. Gauguet, J. Fils, A. Virdis, A. Clairon, N. Dimarcq, Ch. J. Borde, A. Landragin, 6-axis inertial sensor using cold-atom interferometry, Phys. Rev. Lett. 97, 010402 (2006).
  * (13) M. S. Ahmed, D. V. Cuk, Comparison of Different Computation Methods for Strapdown Inertial Navigation Systems, Sci. Tech. Rev. LV, 22 (2005).
  * (14) C. Jekeli, ”Inertial navigation systems with geodetic applications,” de Gruyter, Berlin, New York, 2001.
  * (15) Aaron Canciani. Integration of cold atom interferometry ins with other sensors. Master’s thesis, Second Lieutenant, Air Force Institute of Technology (USAF), 2012.
  * (16) J. Lautier, L. Volodimer, T. Hardin, S. Merlet, M. Lours, F. Pereira Dos Santos, and A. Landragin, Hybridizing matter-wave and classical accelerometers, Appl. Phys. Lett. 105, 144102 (2014).
  * (17) P. Cheiney, L. Fouché, S. Templier, F. Napolitano, B. Battelier, P. Bouyer, and B. Barrett. Navigation-compatible hybrid quantum accelerometer using a Kalman filter. Phys. Rev. Applied 10, 034030 (2018).
  * (18) Y. Wu, J. Guo, X. Feng, L. Q. Chen, C-H. Yuan, W. Zhang, Atom-light hybrid quantum gyroscope, Phys. Rev. Applied 14, 064023 (2020).
  * (19) X. Wang, A. Kealy, C. Gilliam, S. Haine, J. Close, B. Moran, K. Talbot, S. Williams, K. Hardman, C. Freier, P. Wigley, A. White, S. Szigeti, S. Legge, Enhancing Inertial Navigation Performance via Fusion of Classical and Quantum Accelerometers, arXiv:2103.09378v1 [quant-ph].
  * (20) M. A. Kasevich and B. Dubetsky, Kinematic Sensors Employing Atom Interferometer Phases, US Patent 7,317,184 (2005).
  * (21) G.E.P. Box, , M.E. Muller, A note on the generation of random normal deviates, Annals Math. Stat, 29, 610 (1958).
  * (22) M. A. Kasevich and B. Dubetsky, The phase of an atom interferometer as a direct source for precise navigation, private communication.
  * (23) N. F. Ramsey, A new molecular beam resonance method, Phys. Rev. 76, 996 (1949).
  * (24) B. Dubetsky and M. A. Kasevich, Atom interferometer as a selective sensor of rotation or gravity, Phys. Rev. A 74, 023615 (2006).
  * (25) K. Bongs, R. Launay, and M. A. Kasevich, High-order inertial phase shifts for time-domain atom interferometers, Appl. Phys. B 84, 599 (2006).
  * (26) J. M. Hogan, D. M. S. Johnson, M. A. Kasevich, Light-pulse atom interferometry, arXiv:0806.3261 [physics.atom-ph], appear in the Proceedings of the International Summer School of Physics ”Enrico Fermi” on Atom Optics and Space Physics (Varenna, July 2007).
  * (27) S. Kleinerta, E. Kajaria, A. Rouraa, W. P. Schleich, Representation-free description of light-pulse atom interferometry including non-inertial effects, Phys. Rept. 605, 1 (2015).
  * (28) G. L. Kotkin and V. G. Serbo, ”Collection of problems in classical mechanics,” Oxford, New York, Pergamon Press (1971), problem 6.23.
  * (29) C. K. Chui, G. Chen, ”Kalman filtering with real-time applications,” Springer, Berlin, Heidelberg (1999), Sec. 1.1
  * (30) O. Rodrigues, ”Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace, et de la variation des coordonnées provenant de ces déplacements considérés indépendants des causes qui peuvent les produire”, Journal de Mathématiques Pures et Appliquées 5, 380 (1840).
  * (31) V. S. Shebshaevich, P. P. Dmitriev, N. V. Ivancevich, A. V. Kalugin, E. G. Kovalevsky, I. V. Kudryavtsev, V. Yu. Kutikov, Yu. B. Molchanov, Yu. A. Maksyutenko, Setevye sputnikovye radionavigatsionnye sistemy, RADIO I SVYAZ’, Moscow 1993, p. 220 (in Russian).
  * (32) W.A. Heiskanen and H. Moritz, ”Physical Geodesy,” W.H. Freeman and Co. San Francisco (1967), Sec. 2-9
  * (33) E. Groten, ”Parameters of Common Relevance of Astronomy, Geodesy, and Geodynamics,” pp. 134-140 in H. Moritz, Geodetic reference system 1980, Journal of Geodesy, 74, 128 (2000).
  * (34) Committee on Earth Gravity from Space, in Satellite Gravity and the Geosphere: Contributions to the Study of the Solid Earth and Its Fluid Envelopes (National Academies Press, Washington, DC, 1997), p. 13.
  * (35) Y. M. Wang, GSFCO0 mean sea surface, gravity anomaly, and vertical gravity gradient from satellite altimeter data, J. Geophys. Res. 106, 31167 (2001).
  * (36) J. Yu, C. Jekeli & M. Zhu, Analytical solutions of the Dirichlet and Neumann boundary-value problems with an ellipsoidal boundary, Journal of Geodesy 76, 653–667 (2003).
  * (37) C. Jekeli, J. K. Lee, J. H. Kwon, Modeling errors in upward continuation for INS gravity compensation, Journal of Geodesy, 81, 297 (2007).
  * (38) Y. M. Wang, Geodetic Boundary Value Problems, In Encyclopedia of Geodesy, E.W. Grafarend (ed.), Springer International Publishing Switzerland (Outside the USA) 2016.