2.1. Calculation of the equatorial and horizontal coordinates of navigation celestial bodies

The main problem in the development of the Astronomy Package was the development and software implementation of the mathematical apparatus for calculating the equatorial coordinates (right ascension and declination) of navigation celestial bodies: the Sun, the Moon, the brightest (navigation) stars and navigation planets (Venus, Mars, Jupiter and Saturn).

The equatorial coordinates of the stars are not constant. They change due to the body’s motion along the celestial sphere, precession and nutation of the Earth’s axis, as well as an aberration. Based on this, the calculation of the equatorial coordinates of celestial bodies in a general form can be represented by the following formulae:

(1) $$\left\{\matrix{{\alpha = \tilde \alpha {\;_m} + \Delta {\alpha _{AB}},} \cr {\delta = \tilde \delta {\;_m} + \Delta {\delta _{AB}},} \cr} \right.$$

where ${{\tilde \alpha }}_{m},{{\tilde \delta }}_{m}$ are mean places of the stars for the current epoch, considering precession and nutation; and $\mathrm{\Delta }{\alpha }_{AB},\mathrm{\Delta }{\delta }_{AB}$ are corrections for aberration. Here, ${\alpha }_{m},{\delta }_{m}$ are mean places of the celestial bodies at the current epoch without considering precession and nutation, for different types of celestial bodies calculations carried out in different ways. Calculation of mean places, considering precession, nutation and aberration was considered in the article by Kozik et al. (Reference Kozik, Denisova and Alcybeev2020).

Having calculated the values $\alpha $ and $\delta $ , we can calculate the horizontal coordinates for a given geographic latitude $\varphi$ : $h$ is the altitude of the celestial body above the horizon and $A$ is the azimuth. Altitude $h$ is expressed from the following formula:

(2) $${\rm sin}\ h={\rm sin}\ \varphi\ {\rm sin}\ \delta +{\rm cos}\ \varphi\ {\rm cos}\ \delta\ {\rm cos}\ {t}_{\rm{local}},$$

obtained from the navigation triangle. In equation (2), ${t}_{\rm{local}}$ is the local hour angle.

Several formulae are known for the calculation of $A$ , one of which is obtained from the formula

(3) $${\rm cot}\ A\ {\rm sin}\ {t}_{\rm{local}}={\rm cot}\ {\Delta }\ {\rm sin}\ \vartheta -{\rm cos}\ \vartheta\ {\rm cos}\ {t}_{\rm{local}},$$

where $\mathrm{\Delta }=90^\circ -\delta$ (polar distance), $\vartheta =90^\circ -\varphi$ . From equation (3), we obtain

(4) $${\rm{cot}}\ A = {{{\rm{tan}}\ \delta\ {\rm{cos}}\ \varphi } \over {{\rm{sin}}\ {t_{\rm{local}}}}} - {{{\rm{sin}}\ \varphi } \over {{\rm{tan}}\ {t_{\rm{local}}}}}.$$

The advantage of equation (4) is that all the necessary variables for the calculation are known, but it also has a drawback: the resulting value must be factored by a quarter since the function ${\rm{cot}}\left(\cdot \right)$ has a period $\tau =\pi $ .

There is a formula obtained from the cosine theorem,

(5) $${\rm{cos}}\ A = {{{\rm{sin}}\ \delta } \over {{\rm{cos}}\ h\ {\rm{cos}}\ \varphi }} - {\rm{tan}}\ h\ {\rm{tan}}\ \varphi .$$

When using equation (5), all the variables necessary for the calculation are known, since it is calculated by equation (2) in the previous step. In addition, since the function ${\rm{cos}}\left(\cdot \right)$ has a period $\tau =2\pi $ , it is necessary to look only at for which side of the horizon the value $A$ is in the west or in the east (unlike equation (4), where, as noted earlier, it is necessary to determine the correct quadrant of $A$ ), it depends on the value in equation (2), if ${t}_{\rm{local}}\lt 180^\circ$ , then $A:=360^{\circ}-A$ , if ${t}_{\rm{local}}\gt 180^\circ$ , then $A:=A$ .

2.2. Calculation of the time of the apparent sunrise (sunset) and solar illumination

In practice, it is important to know the time of sunrise and sunset for orientation during the daylight and nighttime hours. At the moment of apparent sunrise, the centre of the Sun is below the true horizon, and its altitude can be computed using the formula

(6) $$h=-d+p-{R}_{\underset{}{\odot }}-\rho, $$

where $d$ is the value of the dip of the horizon, $d=1.76\sqrt{{\tilde e}}$ , where ${\tilde e}$ is the height of the observer’s eye above sea level in metres;

$$p = {{180} \over \pi } \cdot 60 \cdot {\rm{arcsin}}{{{a_e}} \over {AU \cdot r}}$$

is the parallax of the Sun in arc minutes, where ${a}_{e}=6378140$ metres is the equatorial radius of the Earth, $AU$ is astronomical unit in metres, $r$ is the distance to the Sun in the astronomical units, which is calculated according to the theory Variations Séculaires des Orbites Planétaires (VSOP) (Bretagnon and Francou, Reference Bretagnon and Francou1988; Guo and Xu, Reference Guo and Xu2015);

$${R_{\mathop \odot \limits_{} }} = {{180} \over \pi } \cdot 60 \cdot {\rm{arcsin}}{{{a_{\mathop \odot \limits_{} }}} \over {AU \cdot r}}$$

is the apparent angular radius of the Sun in minutes, where ${a}_{{\odot }_{}}=6.95700\times 1{0}^{8}$ metres is the radius of the Sun; $\rho $ is astronomical refraction, raising the image of the Sun.

Astronomical refraction for different apparent altitudes ${h}_{a}$ in navigation practice is usually calculated by finding a value for normal conditions (temperature $t^{\circ}=10\, \mathrm{}^\circ \rm{C}$ , atmospheric pressure $B=760\,\mathrm{}\rm{mmHg}$ ) and then considering corrections for different temperatures and atmospheric pressure using special tables

$$\rho ={\rho }_{0}+\mathrm{\Delta }{\rho }_{t^{\circ}}+\mathrm{\Delta }{\rho }_{B}.$$

Since, in this case, it is necessary to calculate the value of astronomical refraction at the time of sunrise (sunset) of the Sun, at which the apparent altitude is equal to the dip with the opposite sign $-d$ , the value ${\rho }_{0}$ is selected from Table 1 or calculated by the formula

$${\rho }_{0}=34.4+0.002{d}^{2}+0.2d.$$

Table 1.

Correspondence of astronomical refraction for normal conditions ${\rho }_{0}$ to the dip $d$

Temperature correction values $\mathrm{\Delta }{\rho }_{t^\circ }$ are selected from the Table 2 or calculated by the formula

$$\Delta {\rho _{{t^{\circ} }}} = \Delta \rho _{{t^{\circ} }}^{0}\left( {1 - {d \over {10}}} \right) + \Delta \rho _{{t^{\circ} }}^{1}{d \over {10}},$$

where

$$\matrix{{{\rm{\Delta }}\rho _{t^\circ }^0 = 0.0011{t^{{\circ}2}} - 0.2355t^{\circ} + 2.2,{\rm{\;\;\;\;\Delta }}\rho _{t^{\circ}}^1 = 0.0013{t^{{\circ}2}} - 0.2638t^{\circ} + 2.4.} \hfill \cr }$$

Table 2.

Correspondence of temperature correction values $\Delta {\rho }_{{t}^{\circ}}$ to the dip $d$ at the corresponding temperature $t^{\circ}$

Pressure correction values $\mathrm{\Delta }{\rho }_{B}$ are selected from the Table 3 or calculated by the formula

$$\Delta {\rho _B} = \left\{ {\matrix{ {0.049\left( {B - 760} \right),\ \ \ d \lt 5,} \hfill \cr {0.053\left( {B - 760} \right),\ \ \ d \gt 5.} \hfill \cr } } \right.$$

Table 3.

Correspondence of pressure correction values $\Delta {\rho }_{B}$ to the dip $d$ at the corresponding temperature $t^{\circ}$

The hour angle is directly related to time, thus, by calculating the value of the local hour angle ${t}_{\rm{local}}$ at the time of sunrise (sunset), we can calculate the time of its rise (sunset).

From equation (2), we have

(7) $${\rm{cos}}\ {t_{\rm{local}}} = {{{\rm{sin}}\ h - {\rm{sin}}\ \varphi\ {\rm{sin}}\ \delta } \over {{\rm{cos}}\ \varphi\ {\rm{cos}}\ \delta }}.$$

Since the moment of the sunrise (sunset) is not known and therefore the declination of the Sun $\delta $ at the moment of the rise (sunset) is unknown, based on this, it will not be possible to directly use equation (7). Then, we will search ${t}_{\rm{local}}$ by the method of successive approximations. Finding an approximate value ${t}_{\rm{local}}$ is carried out according to the following algorithm.

1. 1. At the moment ${T}_{\rm{gr}}^{i}$ , the values $\alpha $ and $\delta $ and the value of the Greenwich hour angle of the first point of Aries are calculated from equation (1). Following the resolutions XVII and XVIII of the General Assemblies of the IAU, the Greenwich Mean Sidereal Time ${S}_{0}^{m}$ in ${0}^{\rm{h}}$ UT1 is determined by the following formula in seconds of time:

   (8) $${S}_{0}^{m}=24110.54841 + 8640184.812866T + 0.093104{T}^{2}-6.2\times 1{0}^{-6}{T}^{3},$$

   where $T$ is the time counted in Julian centuries by 36 525 days in UT1 from the ${1}^{\rm{st}}$ Jan. 2000 epoch, $1{2}^{\rm{h}}$ UT1 (JD 2 451 545.0). Considering the value ${S}_{0}^{m}$ from equation (8), the value is calculated by the formula

   

   where $\mathrm{\Delta }\psi $ is the nutation of the first point of Aries in ecliptic longitude, and ${T}^{d}$ is the time from the beginning of the day, expressed in fractions of a day.

   As an initial approximation ${T}_{\rm{gr}}^{0}$ , it is assumed

   $${T}_{\rm{gr}}^{0}={0}^{\rm{h}}-N,$$

   where $N$ is the time zone number.

2. 2. According to equation (7), one can calculate ${t}_{\rm{local}}$ at the time of sunrise ${t}_{\rm{local}}:={t}_{\rm{local}}$ , and at the time of sunset ${t}_{\rm{local}}:=360^{\circ}-{t}_{\rm{local}}$ .

3. 3. The Greenwich hour angle of the Sun is calculated by ${t}_{\rm{gr}}={t}_{\rm{local}}-\lambda$ .

4. 4. The updated Greenwich hour angle of the first point of Aries is calculated using .

5. 5. Use the calculated value $\mathrm{\Delta }{T}^{\rm{h}}$ to approximate value ${T}_{\rm{gr}}^{i+1}$ . The value $\mathrm{\Delta }{T}^{\rm{h}}$ is found by dividing the difference between the hour angles by the rate of change. Since the hour angle of the first point of Aries changes for $15^{\circ}02.5^{\prime}$ by an hour, the formula for calculating the value $\mathrm{\Delta }{T}^{\rm{h}}$ has the following form:

   

6. 6. A new time approximation ${T}_{\rm{gr}}^{i+1}={T}_{\rm{gr}}^{i}+\mathrm{\Delta }{T}^{\rm{h}}$ is calculated. If

   (9) $$\left|\mathrm{\Delta }{T}^{\rm{h}}\right|\gt \varepsilon, $$

   return to step 1 and calculate the values $\alpha $ , $\delta $ and at the moment ${T}_{\rm{gr}}^{i+1}$ ; in the case where

   (10) $$\left|\mathrm{\Delta }{T}^{\rm{h}}\right|\lt \varepsilon, $$

   we exit the loop and accept ${T}_{\rm{gr}}$ as the last value of the time. In conditions (9) and (10), $\varepsilon$ is some set value characterising the permissible error.

The calculation of solar illumination, rising (setting) of planets, stars and the Moon is carried out according to the same algorithm, the only difference is in how the altitude $h$ is calculated at the time of the event.

To calculate the time of the beginning of morning civil twilight and the end of evening civil twilight: $h=-6^\circ $ , to calculate the time of the beginning of morning nautical twilight and the end of evening nautical twilight: $h=-12^\circ $ .

When calculating the altitude of the Moon at the time of the apparent rising of its upper limb, the same equation (6) is applied, but using the following values: the mean parallax of the Moon ${p}_{\rm{Moon}}=57^{\prime}$ (the exact parallax can be calculated using different methods (Chapront, Reference Chapront1995; Chapront and Francou, Reference Chapront and Francou2003; Escalona-Buendia and Pina, Reference Escalona-Buendía and Pina2002)) and ${R}_{\rm{Moon}}=15.5^{\prime}$ .

Since the radius of the planets is zero, the formula for calculating the altitude is simplified to the form

(11) $$h=-\rho +p.$$

In the case of calculating the altitudes of the planets, equation (11) is simplified to $h=-\rho $ , but since the rising of planets and stars above the apparent horizon is an almost imperceptible event, in practice, one can take the altitude $h=0$ .

2.3. Identification of nautical planets and nautical stars

To identify stars and planets by horizontal coordinates, the inverse problem is solved. By horizontal coordinates $h$ and $A,$ equatorial coordinates $\alpha $ and $\delta $ are calculated. Value $\alpha $ is calculated by the formula

where is the Greenwich hour angle of the first point of Aries, which is calculated at the time of observations; ${\lambda }_{W}^{E}$ is reckonable longitude, in accordance with the sign. The value ${t}_{\rm{local}}$ is expressed from the ratio

$${\rm{cos}}\ {t_{\rm{local}}} = {{{\rm{sin}}\ h} \over {{\rm{cos}}\ \varphi\ {\rm{cos}}\ \delta }} - {\rm{tan}}\ \varphi\ {\rm{tan}}\ \delta, $$

obtained similarly to equation (5).

The value $\delta $ is expressed from the ratio

$${\rm {cos}}\ \delta ={\rm sin}\ h\ {\rm sin}\ \varphi +{\rm cos}\ h\ {\rm cos}\ \varphi\ {\rm cos}\ A,$$

similar to the calculation of the value by equation (2). Based on the found values of the equatorial coordinates, one can identify the star, for example, from the Yale Bright Star Catalog (2017). The identified celestial body can be a planet, therefore, to obtain accurate information, it is necessary to calculate the equatorial coordinates of all visible planets and check for coincidences.

As a rule, the reckoning position has an error, which means that it is impossible to get an exact match. That inaccuracy initiated a new technique through which it is possible to estimate the probability of finding a star within the specified error ellipse. Based on the expected errors in altitude and azimuth, an ellipse is set, where the semi-minor axis is the error in the altitude $b=\mathrm{\Delta }h$ . The semi-minor axis $a=\mathrm{\Delta }A\cdot \rm{cos}\ h$ . The probability of identifying a star within the ellipse is calculated by the distance of the star from the centre of the ellipse using the Laplace function

$$\Phi \left( x \right) = {2 \over {\sqrt {2\pi } }} \int _0^x{\mkern 1mu} {e^{ - {{{z^2}} \over 2}}}dz$$

as a Taylor series

$$\Phi \left( x \right) = {2 \over {\sqrt {2\pi } }}\sum\limits_{q = 0}^\infty {\mkern 1mu} {{{{\left( { - 1} \right)}^q}} \over {{2^q} \cdot q!}} \cdot {{{x^{2q + 1}}} \over {2q + 1}} = {2 \over {\sqrt {2\pi } }}\left( {x - {{{x^3}} \over 6} + {{{x^5}} \over {40}} - {{{x^7}} \over {336}} + \ldots } \right).$$

2.4. Position determination by lines of position

In the case of automated positioning based on observations with respect to reference points, the most probable coordinates are determined by the least squares (LS) method. In the presented work, to determine the location, the LS method was used, implemented by solving equations expressed in matrix form.

Let $n$ be the number of observed navigation parameters (or lines of position) used for determining the ship’s position. The required corrections to the reckoning position are calculated by solving the matrix equation

$${\cal X}={\left({\cal{A}}^{T}{\cal{K}}^{-1}\cal{A}\right)}^{-1}{\cal{A}}^{T}{\cal{K}}^{-1}{\cal L}={\left({\cal{A}}^{T}{\cal C}{\cal A}\right)}^{-1}{{\cal A}}^{T}{\cal C}\cal{L},$$

where $\cal{X}$ is the matrix of the searched corrections to the reckoning position

$${\cal X} = \left[ {\matrix{ {\Delta \varphi } \hfill \cr {\Delta w} \hfill \cr } } \right],$$

and $\Delta\varphi, \Delta w$ are the desired corrections to the calculated coordinates, matrix $\cal{A}$ of size $n\times 2$ is the coefficients matrix of equitations of lines of position

$${\cal A} = {\left[ {\matrix{ {{a_1}} \hfill & 38; {{a_2}} \hfill & 38; \ldots \hfill & 38; {{a_n}} \hfill \cr {{b_1}} \hfill & 38; {{b_2}} \hfill & 38; \ldots \hfill & 38; {{b_n}} \hfill \cr } } \right]^T},$$

matrix $\cal{K}$ of size $n\times n$ is the correlation matrix of error of navigation parameters

$${\cal K} = \left[ {\matrix{ {m_1^2} \hfill & 38; {{K_{12}}} \hfill & 38; {{K_{13}}} \hfill & 38; \ldots \hfill & 38; {{K_{1n}}} \hfill \cr {{K_{21}}} \hfill & 38; {m_2^2} \hfill & 38; {{K_{23}}} \hfill & 38; \ldots \hfill & 38; {{K_{2n}}} \hfill \cr \ldots \hfill & 38; \ldots \hfill & 38; \ldots \hfill & 38; \ldots \hfill & 38; \ldots \hfill \cr {{K_{n1}}} \hfill & 38; {{K_{n2}}} \hfill & 38; {{K_{n3}}} \hfill & 38; \ldots \hfill & 38; {m_n^2} \hfill \cr } } \right].$$

Its diagonal elements are the total dispersion of errors of navigational parameters ${m}_{i}^{2}$ , and non-diagonal elements are correlation points ${K}_{ij}$ , $i$ -th and $j$ -th navigational parameters (for navigational parameters of one type, ${K}_{ij}={m}_{0}^{2}$ ), $\cal{C}$ is the invertible matrix of the correlation matrix $\cal{K}$ , ${\cal C}={\cal{K}}^{-1}$ .

The multiplication of matrices ${\cal{A}}^{T}\cal{C}\cal{A}$ is a matrix of the absolute term of the normal equation of the LS method. Invertible matrix ${\left({\cal{A}}^{T}\cal{C}\cal{A}\right)}^{-1}$ gives correlation matrix of errors of calculated corrections is matrix ${\tilde {\cal K}}$

$${{\cal A}^T}{\cal C}{\cal A} = \left[ {\matrix{ {{A_1}} \hfill & 38; {{B_1}} \hfill \cr {{A_2}} \hfill & 38; {{B_2}} \hfill \cr } } \right] = \left[ {\matrix{ {\mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{a_j}{a_i}} \hfill & 38; {\mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{a_j}{b_i}} \hfill \cr {\mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{b_j}{a_i}} \hfill & 38; {\mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{b_j}{b_i}} \hfill \cr } } \right],$$

$$\tilde {\cal K} = {\left( {{{\cal A}^T}{\cal C}{\cal A}} \right)^{ - 1}} = {1 \over {{\rm{det}}\left( {{{\cal A}^T}{\cal C}{\cal A}} \right)}}\left[ {\matrix{ {{B_2}} \hfill & 38; { - {A_2}} \hfill \cr { - {B_1}} \hfill & 38; {{A_1}} \hfill \cr } } \right] = \left[ {\matrix{ {m_\varphi ^2} \hfill & 38; {{K_{\varphi w}}} \hfill \cr {{K_{\varphi w}}} \hfill & 38; {m_w^2} \hfill \cr } } \right].$$

Matrix $\cal{L}$ is a matrix of absolute terms ${l}_{i}$ , $i=1,2,\ldots,n$ , of equations of position lines

$${\cal L} = {\left[ {\matrix{ {{l_1}} \hfill & 38; {{l_2}} \hfill & 38; \ldots \hfill & 38; {{l_n}} \hfill \cr } } \right]^T},$$

thus, the multiplications of matrices ${\cal{A}}^{T}\cal{C}\cal{A}$ is a matrix of absolute terms of normal equation of the LS method

$${{\cal A}^T}{\cal C}{\cal L} = \left[ {\matrix{ {\mathop \sum \limits_{j = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{a_j}{l_i}} \hfill \cr {\mathop \sum \limits_{i = 1}^n \mathop \sum \limits_{i = 1}^n {c_{ij}}{b_j}{l_i}} \hfill \cr } } \right].$$

With mutually independent navigation parameters, the matrices $\cal{K}$ and $\cal{C}$ become diagonal matrices, while the diagonal elements of the matrix $\cal{K}$ become equal to the partial variances, and the diagonal elements of the matrix $\cal{C}$ become equal to the weights of the corresponding navigation parameters. All non-diagonal elements vanish in this case (Gruzdev, Reference Gruzdev1979, Reference Gruzdev2005).

Thus, the automated positioning and assessing its accuracy required not only observed navigation parameters, but also the enumerable navigation parameters, gradients (their modules, directions and projections), and a priori estimates of the accuracy of the navigation parameters.

Reckoned navigation parameters and their gradients when using the LS method are found by using the inverse geodesics computation problem (Gruzdev, Reference Gruzdev2005).