Software package for nautical astronomy

Abstract

The article is devoted to the mathematical theories and algorithms necessary for the implementation of a software package that fully automates the calculations necessary in Nautical Astronomy. The article describes a method for calculating the equatorial and horizontal coordinates of the celestial bodies at any moment of time. The authors describe the calculation of the time of the apparent rising (setting) of the Sun, solar illumination and events of other celestial bodies. A formula for calculating astronomical refraction is proposed. A matrix method for implementing the method of least squares for determining the coordinates of a place along the lines of position is described. An algorithm for identifying navigational planets is also described and a method for estimating the error for it is proposed. Based on this, the results of the development of the software package ‘Astronomy Package’ for Nautical Astronomy are presented.

1. Introduction

Nautical Astronomy (Celestial Navigation) includes a set of practical methods necessary to determine the coordinates of a ship’s position in the sea by the measured altitudes of celestial bodies. Despite the active development of Satellite Navigation, Nautical Astronomy is also developing (Pierros, 2018; Lušić, 2018). There has also been the development of digital sextants (Kozik et al., 2020, 2021; Li et al., 2021). Nautical Astronomy, by decision of the International Maritime Organization, is included in the International Convention on Standards of Training, Certification and Watchkeeping for Seafarers (IMO, 2012), as a necessary competence for all watchkeeping officers (Gagarskiy, 2015). This decision is because the means of Nautical Astronomy can be used as a backup method for determining the coordinates of a ship’s position, especially in emergency and force majeure situations (Leonov et al., 2015). Additionally, Nautical Astronomy has several advantages: widespread availability, complete autonomy and sufficient accuracy for the high seas. However, the use of celestial bodies as navigational reference points is possible only if there is a catalogue with the coordinates of celestial bodies on board. For example, the Marine Astronomical Yearbook published by the Institute of Applied Astronomy of the Russian Academy of Sciences, The Nautical Almanac published jointly by the US Naval Observatory and the UK Hydrographic Office, Ephemerides Nautiques published by the French Bureau des Longitudes, and others. In addition to specialised catalogues, calculations require a significant amount of time.

Section B-II/1, paragraph 20 of the International Convention on the Training, Certification and Watchkeeping of Seafarers, 1978, as amended in 2010 (Manila Amendments), states: ‘Training in Celestial Navigation may include the use of electronic nautical almanac and Celestial Navigation calculation software’ (IMO, 2012).

Nevertheless, there are currently no software products that can adequately automate positioning by celestial bodies, but active development is in progress in this area. We might note the development of programs by the Institute of Applied Astronomy of the Russian Academy of Sciences (Lukashova et al., 2018; Sveshnikov et al., 2016), as well as mobile applications for the Android operating system, which cannot be approved for use on ships due to the peculiarities of international certifications. There are other good examples, which however are not suitable for solving problems on the ship (Sadiq et al., 2017; Bensky, 2018, Kozik and Sibilev, 2021). Proceeding from this idea, we have decided to develop a software package that can adequately automate the calculation of all possible tasks of Nautical Astronomy.

The article presents the main mathematical methods that were used as the basis of the developed software package. Methods include: calculation of the equatorial and horizontal coordinates of navigation celestial bodies; calculation of the time of the apparent rising (setting) of navigation celestial bodies, solar illumination and events of others celestial bodies; positioning along the lines of position; and an algorithm for identifying navigational planets and stars.

The article proposes the results of the development of the software package ‘Astronomy Package’ for Nautical Astronomy. The software package is developed in C++ for the Windows XP/7/8/8.1/10 operating system. The package includes the following programs: Astronomical almanac, StarFinder, Position, Astroblank, Sunrise, Tables of Altitudes and Azimuths, Calendarium (this article does not describe this program), and Auxiliary celestial sphere.

2. Methods

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. (2020).

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, 1988; Guo and Xu, 2015);

$${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, 1995; Chapront and Francou, 2003; Escalona-Buendia and Pina, 2002)) 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, 1979, 2005).

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, 2005).

3. Results

Astronomy Package software was developed in C++ for the Windows XP/7/8/8.1/10 operating systems and designed to automate computational processes during determining a position using Nautical Astronomy. The proposed software package includes the following software components: Astronomical almanac, Position, Astroblank, Sunrise, StarFinder, Tables of Altitudes and Azimuths, Calendarium, and Auxiliary celestial sphere.

3.1. Astronomical almanac

The Astronomical almanac is designed to calculate the equatorial and horizontal coordinates of celestial bodies and other auxiliary values required in Nautical Astronomy. The program completely replaces Nautical Almanacs (The Nautical Almanac, Brown’s Nautical Almanac, etc.). Calculations in the program are carried out according to the algorithms and methods from Sections 2.1 and 2.2. The program user interface is shown in Figure 1.

Figure 1. Example of the user interface of the Astronomical almanac. Tab with parameters for planets. Example with parameters for Venus on November 11, 2018, at 17:19:07, for latitude $59^{\circ} 59.0^{\prime}\, \rm{N}$ and longitude $29^{\circ} 46.0^{\prime}\, \rm{E}$ .

The Astronomical almanac program allows you to calculate the parameters for Sun; Moon; navigation planets (Venus, Mars, Jupiter, Saturn), in addition to navigation planets, there is also the possibility of calculating parameters for Mercury, Neptune, Uranus and the dwarf planet Pluto; 9 096 bright stars from the Yale Bright Star Catalog (2017), including 160 navigation stars from the Nautical Astronomical Yearbook and 173 navigation stars from The Nautical Almanac or Brown’s Nautical Almanac.

The working area of the program consists of four tabs: Sun, Moon, Planets, Stars. In each tab, a function is implemented for calculating the parameters for the corresponding celestial bodies. In the program, the parameters are indicated: RA ( $\alpha $ ) is right ascension, Dec ( $\delta $ ) is declination, SHA ( ${\tau }^{\mathrm{*}}=360^{\circ}-\alpha$ ) is sidereal hour angle, GHA ( ${t}_{\rm{gr}}$ ) is Greenwich hour angle, LHA ( ${t}_{\rm{local}}$ ) is local hour angle, R ( $R$ ) is the angular radius, Aries () is Greenwich hour angle of the first point of Aries, EqT ( $\eta $ ) is the equation of time, HP ( ${p}_{0}$ ) is parallax, m ( $m$ ) is the magnitude, Alt ( $h$ ) is the altitude of the celestial body above the horizon, Az ( $A$ ) is the azimuth of the star. When Alt and Az are highlighted in yellow if apparent horizontal coordinates were used and not highlighted if true horizontal coordinates. There is a possibility of calculating various events according to the algorithms described in Section 2.2.

3.2. StarFinder

The StarFinder program is designed to identify a star (stars and planets) by the measured parameters of altitude and azimuth. Calculations in the program are carried out according to the algorithms and methods from Section 2.3. The program user interface is shown in Figure 2.

Figure 2. Example of the user interface of StarFinder. Example in which the Epsilon Ursae Majoris (Alioth) star is defined. A sextant was used as a measuring device.

The following symbols are presented in the working area: $Z$ is the zenith distance, Dec ( $\delta $ ) is the declination, SHA ( ${\tau }^{\mathrm{*}}$ ) is sidereal hour angle. It is possible to choose used means of observation: a sextant or a telescope. The program uses atmospheric pressure and temperature in the calculations. When using a sextant, the values of altitude and azimuth are used; when using a telescope, the values of altitude and zenith distance are used, and the height of the eye is not considered.

As a result of the calculations, we have equatorial coordinates ${\tau }^{\mathrm{*}}$ and $\delta $ , which are required to identify the stars within the error ellipse.

In the presented example (Figure 2), the identification of the star is carried out based on the altitude and azimuth measured by the sextant on July 6, 2020 at 22:01:38. Location of the theoretical experiment: latitude $59^{\circ} 59.0^{\prime}\mathrm{}\, \rm{N}$ and longitude $29^{\circ}46.0^{\prime}\mathrm{}\, \rm{E}$ . The temperature value is $10\, {}^{\circ}\rm{C}$ . The value of atmospheric pressure is $760\, \mathrm{}\rm{mmHg}$ .

3.3. Position

Calculations for determining the location by the altitudes of the celestial bodies can be quickly and accurately performed on a computer using the Position program. The Position program is used for the determination of the ship position with the lines of position. Calculations in the program are carried out according to the algorithms and methods from Section 2.4. However, other methods are currently emerging that are of interest [see, for example, Lionheart et al. (2020)]. The program user interface is shown in Figure 3.

Figure 3. Example of the user interface of the Position program. Example of determining the ship’s position using four lines of position is shown.

3.4. Astroblank

Even though calculations for determining the location by the altitude can be quickly and accurately performed on a computer, manual calculations play a big role in the training of the officer on watch. The manual calculation requires special training. The Astroblank program is designed to check manual calculations. The program is a Russian form S-8B, in which you need to enter the initial data for the computation. Initially, the program was designed for teachers of Nautical Astronomy, but now it can be used by students for self-control.

The algorithms for calculating the equatorial coordinates of the celestial bodies are the same as in the Astronomical almanac program, but with an imitation of manual calculation using the Nautical Astronomical Yearbook or others. The calculation of the horizontal coordinates and the correction of the altitudes also reproduce the manual calculation, but with the same errors. Therefore, it is not recommended to use this program for a practical positioning, although there will be practically no difference in the final result from an accurate calculation.

For teachers of Nautical Astronomy and captain-mentors, this program comes with example tasks for collective use. Tasks for other combinations of celestial bodies or other years can be supplied separately, on request, including in text form. The program user interface is shown in Figure 4.

Figure 4. Example of the user interface of the Astroblank program.

3.5. Sunrise

The Sunrise program is designed to calculate the event of solar and lunar illumination for a given position and time zone. It is possible to print the results. Calculations in the program are carried out according to the algorithms and methods from Sections 2.1 and 2.2. The program user interface is shown in Figure 5. Calculated on March, 2021, for latitude $59^{\circ} 59.0^{\prime}\mathrm{}\, \rm{N}$ and longitude $29^{\circ} 46.0^{\prime}\mathrm{}\, \rm{E}$ .

Figure 5. Example of the user interface of the Sunrise program.

3.6. TBA

Classical tables for calculating altitudes and azimuths are intended to find the altitude and azimuth of the luminary when determining the position of the vessel at sea. The data of the table have a cumbersome appearance and when working with them, it is very easy to make mistakes.

The TBA (Tables of computed altitude and azimuth, TCA in Russian, sometimes transliterated as TBA) program is designed for self-control when learning to calculate the horizontal coordinates of celestials using Russian tables TBA-52 or TBA-57. The program simulates manual calculation using tables. To demonstrate the user interface, a screenshot of the program with the TBA-52 tab is presented in Figure 6. For example, the input data for program are latitude $33^{\circ} 33.0^{\prime}\mathrm{}\, \rm{N}$ , declination $26^{\circ} 35.0^{\prime}\mathrm{}\, \rm{N}$ and hour angle $00^{\circ} 0.0^{\prime}\mathrm{}\, \rm{E}$ .

Figure 6. Example of the user interface of the TBA program, TBA-52 tab.

3.7. Auxiliary celestial sphere

The auxiliary celestial sphere in Nautical Astronomy is a mathematical model, which allows us to consider the celestial bodies not in space, but on the surface of the sphere. Based on this mathematical model, the authors presented The Auxiliary celestial sphere program. The Auxiliary celestial sphere program is an interactive celestial sphere with the possibility of arbitrary rotation with the mouse cursor or automatic rotation in real time. The program user interface can be used to add various celestial bodies to the sphere and track their movement.

The Auxiliary celestial sphere is designed for approximate solutions of Nautical Astronomy problems and simulates the apparent motion of celestial bodies. The first part of the program (in the horizontal coordinate system) simulates the apparent diurnal motion of the Sun, the mean Sun, the Moon and Earth’s shadow on the Moon’s orbit, Venus, Mars, Jupiter and Saturn, as well as six other arbitrary celestial bodies (stars), whose coordinates are entered manually. The second part of the program (in the equatorial coordinate system) simulates the apparent annual motion of the Sun and the mean Sun, the apparent monthly motion of the Moon, and the apparent motion of the planets. Additionally, this program can be used to solve problems related to converting time from one system to another. Calculations use: Greenwich Mean Time (GMT); local time for any specified longitude; true solar time; standard time for a given longitude; time for any time zone; and true local sidereal time for a given longitude.

The program user interface is shown in Figure 7.

Figure 7. Example of the user interface of the Auxiliary celestial sphere second version.

4. Conclusion

Executed work is focused on various mathematical algorithms necessary for the functioning of the described software. Algorithms are formed and software is implemented. The algorithms include: method for calculating the equatorial and horizontal coordinates of navigational celestial bodies at any moment; calculation of the time of the apparent sunrise (sunset) and solar illumination, where for calculation of astronomical refraction is used the formula for calculating a correction for temperature and atmospheric pressure; a method for finding the local hour angle by the method of successive approximations in an inverse problem; a matrix method for implementing the least squares method for determining the coordinates of a place along the lines of position; an algorithm for identifying navigational planets, for which an error estimation method is proposed.

Thus, based on the formed complex of mathematical theories, the Astronomy Package for Nautical Astronomy has been developed. The software package has been developed in C++ for the Windows XP/7/8/8.1/10 operating system. The package includes the following programs: Astronomical almanac, StarFinder, Position, Astroblank, Sunrise, Tables of Altitudes and Azimuths, Calendarium, and Auxiliary celestial sphere.

The developed software is intended for full automation of the astronomical method of positioning at sea. In the future, the software can completely replace The Nautical Almanac and its analogues.