Generalized Force Approach to Point-to-Point Ionospheric Ray Tracing and Systematic Identification of High and Low Rays

A variant of the direct optimization method for point-to-point ionospheric ray tracing is presented. The method is well suited for applications where the launch direction of the radio wave ray is unknown, but the position of the receiver is specified instead. Iterative transformation of a candidate path to the sought-for ray is guided by a generalized force, where the definition of the force depends on the ray type. For high rays, the negative gradient of the optical path functional is used. For low rays, the transformation of the gradient is applied, converting the neighborhood of a saddle point to that of a local minimum. Knowledge about the character of the rays is used to establish a scheme for systematic identification of all relevant rays between the given points, without the need to provide an accurate initial estimate for each solution. Various applications of the method to isotropic ionosphere demonstrate its ability to resolve complex ray configurations including 3-D propagation and multi-path propagation where rays are close in the launch direction. Results of the application of the method to ray tracing between Khabarovsk and Tory show good quantitative agreement with the measured oblique ionograms.


I. INTRODUCTION
H IGH-FREQUENCY (HF) radio waves are widely used in various applications, including far-distance communication and probing of the ionosphere. For the ionospheric HF waves, the limit of geometrical optics is usually justified, which essentially turns the description of the wave propagation into the problem of identifying relevant ray paths -ray tracing. The theory of ionospheric radio waves based on the eikonal approximation was mostly developed during the 1950s and 60s [1]- [5] providing a tool for numerous ray tracing models including those for magnetoactive ionosphere [6]- [13]. Despite the challenges mainly arising from the complexity of the propagation medium, the prediction of ionospheric rays launched in a specified direction has largely become a welldefined, routine numerical problem, for which many software packages are readily available [6], [8], [11], [14]- [16]. However, in many important applications the launch direction is unknown, but the position of the receiver is specified instead. There, the problem of point-to-point ionospheric ray tracing arises, where relevant radio paths between the transmitter and receiver need to be found. A common approach to this boundary-value problem is the numerical integration of the dynamical equations for the rays combined with the shooting method, where the launch direction is systematically refined until the ray homes in on the specified landing point with the desired precision. The solution is often not unique. For example, when the operating frequency lies in a range from the critical frequency to the maximum usable frequency (MUF) of some of the ionospheric layers, multiple paths represented by high and low rays can occur [17], [18], all contributing to the radio communication between transmitter and receiver. In practice, multipaths are investigated by first scanning over elevation and azimuthal angles at the transmission site in order to isolate sought-for launch directions and provide information for homing-in methods. Care must be taken to send out enough rays in various directions so as not to skip relevant raypaths which satisfy the boundary condition.
Coleman [19] proposed an alternative approach to the pointto-point ionospheric ray tracing problem. The approach makes use of Fermat's variational principle directly and involves iterative transformation of some initially defined trajectory to the optimal configuration satisfying a stationary condition of the optical path functional. By considering only the paths with their endpoints fixed at the transmitter and receiver sites, the boundary conditions are satisfied automatically, which is clearly an advantage of the method. There is, however, a complication arising from different characters of extrema that need to be identified. While high ionospheric rays correspond to minima of the optical path functional and can therefore be found by a straightforward minimization of the objective function, low rays correspond to saddle points [19], [20], which are difficult to locate. Coleman [19] solves the problem of low rays by applying the Newton-Raphson method that requires a good initial guess for the solution. Such an initial estimate of a low ray is, however, not always available, especially for a heavily disturbed ionosphere where the ray paths are quite distorted. Multipaths cause an additional problem as some sort of sampling over initial estimates needs to be organized in order to enable the iterative procedure to converge on the various extrema.
In this article, the direct variational method of Coleman [19] is revisited in a sense that both high and low rays are found as a result of the same optimization procedure guided by a generalized force. The definition of the force depends on the character of the extremal ray paths. For high rays, a negative gradient of the optical path functional is used. For low rays, the gradient is modified so as to effectively convert the neighbourhood of a saddle point to that of a local minimum. In both cases, the component of the gradient parallel to the ray path is substituted by an artificial spring force which does not affect the position of the path in space, but ensures continuity of the solution. Knowledge about the character of the extrema of the optical path functional has also made it possible to develop an efficient computational scheme for global pointto-point ray tracing, where all relevant rays between the fixed points are found one after another in a systematic manner, without the need to provide an accurate initial estimate for each solution.
The article is organized as follows. In Section II, the method for point-to-point ray tracing based on the generalized force approach is described and illustrated using a simplified model of the ionosphere. Some practical advice for a successful implementation of the method is also given. Section III presents some applications of the method, including global search for radio paths between Khabarovsk and Tory in the ionosphere provided by the IRI2007 model and calculation of ionograms of oblique sounding, which are compared with the measured data. Comparison with the traditional shooting method is also given. Summary and outlook are presented in Section IV.

II. METHODS
We start with the equation for the optical path functional. For the isotropic ionosphere, it is given by the following line integral: Here, the integration is performed along the curve γ connecting the endpoints A and B which are associated with the transmitter and receiver of the radio ray, n( r) is the refractive index at point r, and dl is the length element along γ. Following [19], we discretize the functional S[γ] before applying the variational principle. While various methods for discretization are possible, we choose here a simple trapezoidal rule: where r i = (x i , y i , z i ) is the position of the ith vertex, n i = n( r i ), and r 0 = r A , r P +1 = r B . The curve γ is then approximated by a polygonal line connecting P + 2 vertices, where the end points are fixed according to the boundary conditions, but P intermediate vertices need to be adjusted to an optimal configuration representing the radio ray. The optical path functional is turned into a multidimensional function of vertex positions, S = S(r), with r = ( r 1 , r 2 , . . . , r P ) (everywhere in the article, unless stated otherwise, symbols with an arrow above denote three-dimensional vectors, whereas bold symbols denote 3P -dimensional vectors in the space of path configurations). According to Fermat's variational principle, radio paths correspond to extrema of S(r) so that point-topoint ionospheric ray tracing essentially becomes a problem of locating such stationary points.
It is important to realize that Fermat's principle does not specify the character of the stationary points and finding only the paths that make the phase distance minimum might be insufficient for the identification of all relevant rays. Indeed, only high rays and trans-ionospheric rays can be found by straightforward minimization of the optical path functional [19] using, for example, the negative gradient of S -the generalized force. Low rays correspond to saddle points, which are more difficult to locate [20]. The difficulty arises from the need to minimize the phase distance with respect to all but one variable for which a maximization should be carried out and it is not known a priori which variable should be treated differently. The problem of low rays is solved here by properly redefining the generalized force, as described below. The force is then integrated into the same standard optimization techniques used to find the high rays.
Before proceeding to precise definition of the forces that guide the optimization of the phase path, we discuss an important aspect of our approach -the force projection -which we adopt from the nudged elastic band method, originally developed to identify minimum energy paths and tunneling paths of chemical reactions, diffusion processes and magnetic transitions [21]- [25]. The force projection is motivated by the observation that only transverse displacements define position of the path in space. The transverse displacements are generated by the component of the gradient force perpendicular to the path. Therefore, only this component should be included in optimizing the position of the vertices of the polygonal line representation of the radio ray. The perpendicular component of the gradient is obtained by subtracting the parallel component: where ∇ i ≡ ∂/∂ r i . The unit tangent to the path at ith vertex, τ i , is estimated using the position of the neighboring vertices: A displacement of the vertices along the path does not affect the location of the path in space, but still needs to be taken care of so as to ensure good resolution of all parts of the radio ray. The distribution of the vertices along the path can be controlled by introducing artificial springs connecting adjacent vertices [21], [22]. The spring force acting on the ith vertex is evaluated as where κ is the spring constant. The value of κ is chosen here to be the same for all springs, which ensures uniform distribution of the vertices along the path. However, one is free to enforce any distribution as long as it gives an adequate representation of the ray. For example, it seems reasonable to have more vertices where the path is heavily distorted and less vertices where the path only slightly deviates from a straight line. For a particular number of vertices, such uneven distribution should give a better resolution of the radio ray. Uneven vertex distribution can be achieved by use of variable spring constants where the value of κ is chosen proportional to the local curvature of the ray. Because of the force projection, the spring force does not interfere with the gradient force. The former distributes the vertices along the ray, while the latter defines the position of the ray in space. Eqs. (3) -(5) provide necessary ingredients for the definition of the generalized force that iteratively brings the chain of vertices to the ionospheric rays, as described below.

A. Identification of high rays
High ionospheric rays are obtained by simply minimizing the phase distance S(r), see Eq. (2). In this case, the generalized force that acts on the vertices of the polygonalline representation of the path and guides the minimization is defined as the negative gradient of S(r), where the tangential component of the gradient is substituted by the spring force: where with ∇ i S(r) ⊥ and F s i defined in Eqs. (3) -(5). The spring force does not affect the position of the path in space, but ensures continuity of the solution. Initially generated position of the vertices is iteratively adjusted so as to zero the generalized force, F h , and the final, relaxed arrangement of the vertices gives a discrete representation of the high ray. The zeroing of the force can be achieved by various standard techniques including the method of steepest descent, conjugate gradient method, etc. Further discussion about the optimizer choice is presented in Section II-D.
When multiple rays are present between the fixed points, the optimization procedure involved in the ray tracing will most likely converge to the ray closest to the initially generated path. In order to find all relevant rays in such a situation, several initial estimates need to be produced so as to enable the optimization procedure to converge on the various solutions. Providing a good initial guess for each individual ray, particularly for three-dimensional, disturbed ionospheres, could be problematic. It is demonstrated in Section II-C that the issue of defining the initial estimates can largely be solved by adopting a strategy where the character of extremal paths is explicitly taken into account [26]. The rays are then systematically identified one after another and the initial guess for each next ray is determined by a previously found solution, as described below.
An application of the direct variational method to high ionospheric rays is illustrated in Fig. 1. The propagation medium is chosen to be an isotropic ionosphere characterized by two layers, E and F 2. The electron density N e as a function of altitude z is modeled by the following analytical formula [27]: Here, parameters N β , z β and ∆z β are the peak electron density, peak height and thickness of layer β, respectively, with β = 1 corresponding to the E layer and β = 2 corresponding to the F 2 layer. The parameter values are taken to be: km, ∆z 1 = 30 km and ∆z 2 = 75 km. The plasma frequency profile is shown in Fig. S1 of the Supplementary Material available at http://ieeexplore.ieee.org. Operating frequency was chosen to be 12 MHz. P = 13 movable vertices were used in the calculation. Initial estimates for the ray were obtained by placing all intermediate vertices at the same height of either 100 km or 200 km (see Fig. 1). Two local ray searches were then performed and in each search the vertices were iteratively brought from the initial arrangement to the nearest high ray using the velocity projection optimization (VPO) algorithm [21] (see also Appendix) equipped with the generalized force F h . Two high rays were revealed, see Fig. 1.
Depending on the initial vertex arrangement, the calculations converge either to the high ray of the E layer [see Fig. 1(a)] or to that of the F 2 layer [see Fig. 1(b)]. The convergence criterion was taken to be that the force on each vertex had dropped below 10 −9 .

B. Identification of low rays
Low ionospheric rays correspond to saddle points of the optical path functional, therefore neither minimization nor maximization of S(r) can be used to locate them. A special treatment is required and here we apply an adaptation of the method that has become a valuable tool in describing kinetics of thermally activated atomic rearrangements [28], [29]. An iterative optimization procedure involved in the method starts at a minimum of the objective function and follows two stages. First, the neighbourhood of the minimum is escaped by following, for example, the gradient of the objective function. Following the gradient eventually homes in on the maximum, which should be avoided. Therefore, the search direction is redefined at a certain point so as to carry out a maximization of the objective function along only one definite direction and minimization along all other directions. This ensures convergence on the saddle point. The direction along which the objective function needs to be maximized is determined using the eigenvalues and corresponding eigenvectors of the Hessian matrix. An important aspect of the adaptation of this method, referred to as the minimum mode following (MMF) [28], [29], to the ionospheric ray tracing is the evaluation of the Hessian in the subspace orthogonal to the ray. A detailed description of the method is presented below.
An iterative search for the low ray is initiated in the vicinity of a local minimum of S(r) represented, for example, by a previously found high ray or simply by a straight path between the transmitter and receiver located in the unionized surface layer. First, the path needs to be moved away from the convex region near the minimum, the region where all eigenvalues of the Hessian of S(r) are positive, and brought to the basin of attraction of a first order saddle point, where one of the eigenvalues of the Hessian is negative. The escape from the convex region can be accomplished in various ways, but we have noticed in our calculations that following either the gradient of S(r) or the eigenvector of the Hessian corresponding to the lowest eigenvalue gives particularly good results in terms of minimizing the number of iterations needed for convergence.
At a certain iteration of the search, the minimal eigenvalue of the Hessian turns negative, which is an indication that the region near a saddle point has been reached. When this happens, the trick is to define the generalized force, F l , in such a way that it corresponds to the neighbourhood of a local minimum of S(r) rather than a saddle point. This can be done by inverting the component of the gradient along the minimum mode, i.e. the eigenvector of the Hessian corresponding to the minimal eigenvalue, and using that in the definition of the force. Taking into account force projections introduced in Eqs. (3) -(5), the generalized force guiding the optimization of low rays is therefore defined as: Here, λ is the minimal eigenvalue of the Hessian and Q λ is the corresponding normalized eigenvector, the minimum mode. Similar to what is done in the high ray calculations, the vertices are relaxed using some standard optimization technique equipped with the generalized force F l . The final vertex configuration corresponding to zero of F l gives a discrete representation of the low ray. There is no special requirement on the algorithm for zeroing the force in the low ray calculations. The algorithm can be just the same as in the search for high rays. The difference in the identification of high and low rays is only in the definition of the generalized force. The method described above does not require an accurate initial estimate and the iterative search for the low ray can start far away from the solution. Typically, an initial guess for the low ray is chosen by making a small random displacement from the position of the high ray. By generating several such displacements, convergence on various low rays can be achieved. Thus, low ray searches rely on the identification of high rays. The opposite also applies. Indeed, a high ray can be found by displacing the path along the minimum mode from the known low ray solution, followed by a local minimization of the phase distance. The interconnection between high and low rays is the basis of the global ray tracing method presented in Section II-C. According to Eq. (10), identification of low rays relies on the calculation of eigenvalues and eigenvectors of the Hessian of S(r) at each iteration of the optimization procedure. Evaluation of the Hessian requires a careful description. Indeed, only transverse displacements of the path are relevant for the assessment of extrema corresponding to high and low rays. Therefore, the Hessian which describes the character of extrema must be defined in a subspace that is orthogonal to the path. A matrix representation of the Hessian in the orthogonal subspace can be obtained by a projection operator approach. This technique has been used before to define the Hessian in the curved configuration space of magnetic systems [30]. First, the matrix of second-order partial derivatives of S(r) is calculated: Equation (11) defines the Hessian in the 3P -dimensional embedding space of possible path configurations. This needs to be projected on the subspace of transverse displacements of the path in order to obtain a true Hessian, H ⊥ , which can be used to characterize the curvature of S(r) near extremal points: Here, V is a 3P × 2P transformation matrix whose columns represent the basis for the subspace. From the Hessian matrix, H ⊥ , the minimum mode can be evaluated using, for example, the Lanczos algorithm [31]. As the Hessian matrix given by Eq. (12) is defined in the 2P -dimensional subspace of transverse displacements, the evaluation of the minimum mode in the original 3P -dimensional configuration space requires the following transformation: whereQ λ is a 2P -dimensional eigenvector of H ⊥ corresponding to the lowest eigenvalue. The transformation matrix V can be defined for a threedimensional ionosphere as follows. For each vertex i along the path, a pair of orthonormal vectors, η i and ξ i , lying in the plane perpendicular to the path is chosen. For example, η i can be obtained by orthonormalization of a random vector with respect to the path tangent, τ i , and then ξ i can be generated using the cross-product: ξ i = [ η i × τ i ]. η i and ξ i define a basis for transverse displacements of a single vertex i. The corresponding basis transformation matrix, V i , is defined as: where η x i , η y i , η z i and ξ x i , ξ y i , ξ z i are x, y, z-components of vectors η i and ξ i , respectively. The transformation matrix for the whole path is simply a direct sum of V i : Low ray searches are illustrated in Fig. 2. The propagation medium and the radio wave frequency are taken to be the same as in the example presented in Section II-A. Several low ray searches were conducted. For each search, P = 13 movable vertices in the initial path were arranged along the previously found high ray of the E layer and small random displacements were added to start the calculation. The gradient of S(r) was followed to escape the convex region. Once λ, the minimum eigenvalue of the Hessian H ⊥ , turned negative (see corresponding vertex configuration in Fig. 2), the definition of the generalized force was switched, see Eq. (10), so as to enable convergence on the low ray. The VPO algorithm was used to both escape from the convex region and to ultimately zero the force F l . The calculation was considered converged when the force on each vertex had dropped below 10 −9 . The calculations revealed two low rays reflected from either the E layer or the F 2 layer. The latter can also be obtained by carrying out low ray searches initiated at the high ray of the F2 layer (see Fig. S2 of the Supplementary Material available at http://ieeexplore.ieee.org).

C. Global point-to-point ray tracing
Certainty about the character of extremal radio wave paths makes it possible to efficiently identify multiple ionospheric rays using the global optimization method that takes advantage of the observation that there should be at least one saddle point between two local minima of the objective function [29]. In particular, relevant solutions can be systematically explored by traversing from one high ray to another via low rays. For a given high ionospheric ray, several low ray searches are carried out starting from slightly perturbed positions of the vertices, as described in Section II-B. After a low ray has been identified, the path is displaced slightly along the minimum mode in the direction away from the previously found high ray. A new high ray is then found by a local minimization of the phase distance, as described in Section II-A. The algorithm for the global ionospheric ray tracing is presented in greater detail below.
1) For a given position of transmitter and receiver, create an initial path by constructing a chain of intermediate vertices. Initial configuration of the vertices is not critical for the successful outcome of the optimization procedure. For example, the vertices can simply be placed along a straight line connecting the fixed points. 2) Find the first high ray by carrying out a local minimization of the optical path functional (see Section II-A). 3) Generate a starting configuration for a low ray optimization by applying small random displacements to the vertex positions at the high ray minimum. The low ray search starts in the vicinity of the high ray. 4) Perform a low ray search using the MMF method (see Section II-B). 5) Repeat step 3 and step 4 until a predetermined number of unsuccessful attempts to locate new low rays has been exceeded. 6) For each low ray found, generate a starting point for a high ray optimization by slightly moving the vertices along the minimum mode in the direction away from the known high ray. The high ray search starts in the vicinity of the low ray. 7) Find an adjacent high ray by carrying out a local minimization of the optical path functional (see Section II-A). 8) For each new high ray, repeat steps 3 through 7 until all relevant rays have been found. Figure 3 illustrates the global point-to-point ionospheric ray tracing using the method described above. To demonstrate the ability of the method to deal with the three-dimensional propagation, a local inhomogeneity in the F 2 layer was added to the background model used in the examples from Sections II-A and II-B. The resulting electron density distribution is given by the following formula: Here, the profile N e (z) is defined in Eq. (9) with parameters identical to that used in Sections II-A and II-B and the inhomogeneity is described by the Gaussian function: where r d = (500, 500, 300) km is the position of the center of the inhomogeneity and σ = 100 km is the parameter characterizing the size of the inhomogeneity. The ray tracing was performed between the endpoints A and B fixed at r A = (0, 0, 0) km and r B = (1000, 1000, 0) km and the radio wave frequency was chosen to be 10 MHz. In all local searches for the rays, the generalized force was iteratively zeroed using the VPO algorithm. The calculations were considered converged when the force on each vertex had dropped below 10 −9 . The global ray tracing is started by placing intermediate vertices equidistantly along the straight line connecting points A and B. Technically, such straight path contained in the unionized surface layer is a ray corresponding to a local minimum of the phase distance (see ray 0 in Fig. 3). From the initial vertex arrangement, several small random displacements are produced followed by local saddle point searches using the MMF method. Only one low ray is found, ray 1 reflected from the E-layer of the ionosphere (see Fig. 3). After that, a small displacement along the minimum mode is made. Subsequent local minimization reveals ray 2, a high ray of the E-layer. Saddle point searches initiated in the vicinity of ray 2 converge onto three adjacent low rays, one of which, ray 1, has already been identified in the previous local optimizations originated from ray 0. Two newly found solutions, ray 3 and ray 4, are low rays reflected from the F 2-layer (see Fig. 3). Although the rays do not deviate from the vertical plane extending between the endpoints A and B, they are asymmetrical due to the presence of the ionospheric inhomogeneity, with their apices shifted towards either point A (ray 3) or point B (ray 4). When no inhomogeneity is present in the model ionosphere, rays 3 and 4 degenerate into a single, symmetrical low ray (see Fig. S3 of the Supplementary Material available at http://ieeexplore.ieee.org). From ray 3, a small displacement along the minimum mode is made and a new local minimization conducted to converge onto a high ray of the F 2-layer, ray 5. The same process is repeated for ray 4 and another high ray, ray 6, is discovered. Both newly found high rays demonstrate a three-dimensional character of radiowave propagation: They pass round the ionospheric inhomogeneity by deviating from the vertical plane containing the endpoints A and B (see Fig. 3). Saddle point searches launched from ray 5 reveal new low rays, rays 7 and 8, refracted on the ionospheric inhomogeneity. No new low rays are found in the searches originated from ray 6 and the ray tracing procedure is therefore terminated. As a result, eight rays are found by traversing from one minimum of the optical path S( r) to another via saddle points, including three high rays (rays 2, 5 and 6 in Fig. 3) and five low rays (rays 1, 3, 4, 7 and 8 in Fig. 3). The ray tracing is illustrated in the supplementary movie available at http://ieeexplore.ieee.org. The ionospheric inhomogeneity strongly affects the radio wave propagation. Compared to the unperturbed ionosphere, four additional rays appear (see Fig. S3 of the Supplementary Material available at http://ieeexplore.ieee.org for the results of the ray tracing for the unperturbed propagation medium). Two of the rays deviate from the vertical plane exhibiting three-dimensional structure. The example described above demonstrates an important aspect of the ray tracing method presented here. In particular, the knowledge about the character of extremal paths is explicitly used. This knowledge makes the method systematic in a sense that rays are identified one after another, and every local ray search starts from a previously found solution so that there is no need to provide an accurate initial estimate for each ray.

D. Practical advice
Below, we give some practical information concerning the convergence criterion for the ray searches, the choice of the algorithm for zeroing the generalized force, the strategy for escaping the convex region in the low ray searches, magnitude of the spring constant as well as calculations of the minimum mode. This information could be useful for a successful implementation of the ionospheric ray tracing method presented in this article.
Convergence criterion. For a given number of vertices included in the local optimization of the phase distance S(r), the calculation should be considered converged when the magnitude of the generalized force on each vertex drops below the set tolerance. However, even convergence with a tight force tolerance (such as the 10 −9 tolerance specified in the present calculations) may not guarantee a satisfactory resolution of the radio wave ray if the number of vertices is not large enough. On the other hand, including too many vertices in the calculation would result in an unnecessarily high computational effort. A good strategy is to start the ray search with a moderate number of vertices and first converge the path only to a rather high tolerance so as to bring the vertices relatively close to the ray at a reduced computational cost. After that, vertices should be progressively added to the path and optimization of the phase distance reiterated with a low force tolerance until S(r) stops changing. For example, all local ray searches in the test problem presented in Section II-C were started with P = 13 movable vertices, but then the number of vertices was systematically increased up to 39 depending on the ray so as to ensure a proper resolution (the process of adding the vertices to the path is demonstrated in the supplementary movie available at http://ieeexplore.ieee.org).
Zeroing of the force. The ray tracing method presented here does not imply a particular algorithm for zeroing of the generalized force. Any optimization method based only on the evaluation of the first derivatives of the objective function with respect to the configuration parameters can be used. The VPO method presented in the Appendix is stable and reliable, but probably not the most efficient. It is therefore beneficial to switch to some higher-order convergent optimization method when the vertices are relatively close to the ray, while the use of a more stable method, such as the VPO algorithm, is safer when starting from the arrangement of the vertices far away from the sought-for ray. For example, a combination of the VPO method and (limited-memory) Broyden-Fletcher-Goldfarb-Shanno algorithm (see, for example, [32]) should be a good choice. The former can be used to converge the ray only to a rather high tolerance in order to bring the vertices relatively close to the solution, but then a switch to the latter can be made so as to converge quickly on the ray with a low tolerance (such as the 10 −9 tolerance specified in the present calculations).
Escape from the convex region. If a low ray search starts in the vicinity of the high ray, the convex region around the minimum of S(r) needs to be escaped during the first stage of the optimization process. For doing so, two approaches were used [see Eq. (10)]: i) Following the minimum mode; ii) Following the gradient of S(r). It has been noticed that the first method usually brings the path to the neighborhood of the low ray that belongs to the same ionospheric layer as the high ray from which the search has started, while the second method tends to move the vertices to a different layer. In both cases, use of the VPO algorithm to guide the escape is a good strategy.
Minimum mode calculation. Calculation of the minimum mode is an essential part of the ray tracing method presented in this article. The minimum mode needs to be evaluated at each iteration of a low ray search. It is therefore important to perform corresponding matrix operations efficiently in order to minimize the overall computational effort. Fortunately, calculation of only the lowest eigenvalue of the Hessian matrix is required to find the minimum mode (it is actually a good practice to calculate the second lowest eigenvalue together with the lowest one so as to verify that the neighbourhood of the first order saddle point where one and only one eigenvalue of the Hessian is negative has been reached), and there is no need to solve the eigenvalue problem completely. A few of the lowest eigenvalues can be evaluated efficiently using the iterative Lanczos algorithm [31]. Use of the Lanczos algorithm is particulary recommended when the number of vertices included in the calculations becomes large.
Spring constant value. Since the spring force is decoupled from the perpendicular component of the gradient of S(r), the value of the spring constant κ is not critical. In fact, the spring constant can be varied by many orders of magnitude without affecting the results. However, assigning extreme values to κ may slow down convergence. In general, optimal performance of the method is achieved when the spring force is roughly of the same order of magnitude as the perpendicular component of the gradient force. It is recommended to use this criterion when deciding on the value of κ. For the calculations presented here, use of κ within the range from 10 −3 km −1 to 10 −1 km −1 gives particularly good results in terms of minimizing the number of iterations needed for convergence.

III. RESULTS
In this section, the method described above is applied to the radio wave ray tracing between Khabarovsk (47 • 36 , 134 • 42 ) and Tory (51 • 48 , 103 • 04 ). The properties of the propagation medium are specified using the International Reference Ionosphere 2007 (IRI2007) model [33]. For simplicity, no effects of Earth's magnetic field are included in the calculations. Figure 4 shows the simulated propagation at a frequency of f = 12 MHz. In this case, the electron density distribution corresponds to June 22, 2016 at 01:00 UT. Each local ray search in the global ray tracing procedure was started with P = 13 movable vertices, but then the number of vertices was systematically increased up to 41 depending on the ray so as to ensure the proper resolution of the ray. Application of the global ray tracing method reveals five radio wave rays, including three high rays (rays 1, 3, and 5 reflected from layer E, F1, and F2, respectively, as shown in Fig. 4) and two low rays (rays 2 and 4 reflected from layer F1, and F2, respectively, as shown in Fig. 4). The operating frequency is close to the lowest usable frequency (LUF) of the F2 ionospheric layer which results in that two of the rays -high ray of the F1 layer (ray 3) and low ray of the F2 layer (ray 4) -almost merge. Proximity of the rays does not pose any problem for the method.

A. Comparison with the shooting method
For comparison, the radio wave rays between Khabarovsk and Tory have also been calculated using a traditional ray tracing method based on the homing-in approach [10]. First, isolation of sought-for launch directions has been performed by sending out rays from Khabarovsk site and inspecting their landing points. In the calculations, the initial elevation angle was progressively increased from 0 • to 90 • with a predefined step ∆α. The desired launch elevation was considered isolated in the interval from α to α + ∆α if the ray sent at angle α overshot (undershot) the desired landing point (Tory), and the ray launched at α + ∆α undershot (overshot) it. Each isolated launch direction has then been systematically refined using the shooting method to produce rays landing at Tory with an accuracy of ±10 m. The resulting rays are almost identical to that obtained with the direct optimization approach (the rays produced by the two methods are indistinguishable in Fig. 4). However, isolation of all five solutions requires scanning over the launch elevations with a rather small increment, ∆α ≤ 0.1 • . Even if the increment is taken to be as small as 0.5 • , only rays 1, 2, 3 and 5 can be identified, while ray 4 remains unresolved. If ∆α ≥ 1 • , only two rays (1 and 5) are resolved with the traditional ray tracing method, while the other three solutions are lost.
The example above demonstrates how efficient the direct variational method can be in finding multipaths. Identification of the rays that become close in the launch direction (e.g. when the operating frequency approaches LUF or MUF of some of the ionospheric layers) could be problematic for the traditional shooting method as progressively more rays need to be calculated to isolate close solutions, leading to an increase in the computational effort. In contrast, the direct optimization method and global ray tracing technique presented in this article are rather insensitive to the configuration of the rays. In particular, proximity of the rays does not cause any apparent increase in the computational effort involved in the method and the rays can reliably be resolved without the need for fine adjustment of the calculation parameters.

B. Ionograms of oblique sounding
Simulation of oblique ionograms is an important application of point-to-point ionospheric ray tracing. By comparing calculated ionograms with the measured data, quality of ionospheric models can be verified. On the other hand, calculated ionograms can be used to monitor the state of the ionosphere over unpopulated regions.
Both measured and simulated ionograms of oblique ionospheric sounding between Khabarovsk and Tory are shown in Fig. 5. The measurements were performed on January 8, 2016 at 17:11 UT and January 9, 2016 at 05:31 UT using multifunctional ionosounders capable of producing continuous signals with complex phase modulation. The operating frequency range of the sounders is 4-30 MHz, the sweep rate is 500 kHz/s and the radiation power varies in a range from 15 W to 100 W. The ionosounders in Khabarovsk and Tory are a part of the chirp sounding network provided by the Institute of Solar-Terrestrial Physics of the Siberian Branch of the Russian Academy of Sciences. The network covers the East Siberian part of Russia and operates in a yearround monitoring schedule with a 5 minute time resolution enabling measurements of short-term ionospheric processes such as medium-scale travelling ionospheric disturbances [34] and responses of the ionosphere to X-ray solar flares [35]. The ionograms have also been simulated using the ray tracing method presented in this article. For each radio wave ray obtained with the direct variational method, the group delay has been calculated [17], [18] and plotted against the operating frequency (see Fig. 5). The results are in good quantitative agreement with the experimentally measured traces of onehop rays. Some deviations from the experimental results are partially attributed to the effect of Earth's magnetic field which was not included in the theoretical calculations. Particularly sensitive to the magnetic field are high rays [13], which explains generally better agreement between experiment and simulation for the low rays. Moreover, ionospheric weather variations may affect the radio wave propagation, but are not accounted for in the medium provided by IRI2007 model [36]. The neglect of ionospheric weather changes in the theoretical analysis could be another source of the deviations from the experimental results.
While the ionograms presented in this Section are probably too simple to truly demonstrate the performance of the ray tracing methods, they confirm applicability of the generalized force approach to real problems of ionospheric radio wave propagation. The ionograms also illustrate the interpretation of high and low rays in terms of minima and saddle points of the optical path functional. When the operating frequency is significantly different from the MUF, the minimum and the saddle point are far away from each other and the corresponding high and low rays as well as their traces on the ionogram are well separated [see Fig. 5(b) and the corresponding ray configuration for f = 15 MHz in Fig. 6(a)]. The stationary points move towards each other, merge and eventually disappear as the operating frequency approaches the MUF [see Fig. 6(b)]. This interpretation is consistent with the general knowledge about ionospheric propagation, but also provides a valuable insight into the nature of high and low rays. Certainty about the character of high and low rays is at the core of the global ray tracing method presented in this article.

IV. CONCLUSIONS AND DISCUSSION
In this article, we have presented a variant of the direct optimization method for solving the point-to-point ionospheric ray tracing problem. Within the method, the ray searches are guided by the generalized force. The generalized force approach makes it possible to apply a wide range of standard optimizers to identify both high and low ionospheric rays. This formulation represents a significant improvement, particularly for the low rays, which are impossible to locate by a straightforward minimization of the phase distance. Although essentially the same iterative optimization scheme is used to find both high and low rays, the rays are nevertheless discriminated because the definition of the generalized force depends on the ray type. The certainty about the character of the rays has actually made it possible to establish an efficient global ray tracing technique where all relevant rays are identified one after another and every local ray search starts at a previously found solution so that there is no need to provide an accurate initial estimate for each ray. The calculations presented here illustrate the ability of the direct optimization method and global ray tracing technique to resolve complex ray configurations such as three-dimensional propagation and rays which are close in the launch direction. We therefore believe that the method can become a valuable tool in various applications of radio science and possibly other fields of research where geometrical optics approximation can be used.
The developed ray tracing technique is applied here to isotropic propagation media where the phase distance is given by Eq. (1). Anisotropic radio wave propagation due to the presence of the geomagnetic field can be accessed by locating minima and saddle points of a more general functional describing optical path length in a magneto-ionic medium and allowing for the variation in both the radio path geometry and orientation of the wavefront (see e.g. Ref. [19]). Application of the ray tracing method developed in the present study to the radio wave propagation in magnetoactive ionosphere will be presented elsewhere.
Finally, we discuss further extensions of the direct optimization method. Single-hop high and low rays are in the focus of the present study, but multihops can make considerable contributions to the radio communication, too. Traces of multihop rays are evident from the experimentally measured ionograms shown in Fig. 5. Multi-hop propagation can be handled by including points of ground reflection in the optimization procedure. The reflection points should be treated differently compared to other movable vertices representing the ray. In particular, their displacements need to be constrained to lie strictly on the ground surface. A similar idea has already been realized by Coleman [19].
Low rays considered in the present study correspond to first order saddle points of the optical path functional. Such solutions are characterized by a single negative eigenvalue of the Hessian. Other extrema where two or more eigenvalues of the Hessian are negative also satisfy Fermat's principle and, therefore, represent radio wave rays, too. Such higherorder saddle points can be identified in a similar fashion as first order saddles by performing a maximization of the phase distance along directions defined by the Hessian's eigenvectors corresponding to the negative eigenvalues, and minimization along all other directions. Extension of the ray tracing method presented here to higher-order saddle-point-type solutions as well as multi-hop propagation is the subject of future research.

APPENDIX VELOCITY PROJECTION OPTIMIZATION
The velocity projection optimization (VPO) method essentially mimics simulation of damped dynamics of a point mass in the configuration space, where the effect of damping is modeled by a projection of the velocity on the direction of the force. The algorithm has turned out to be highly valuable in many areas of research where optimization problems arise. Here, the VPO optimizer is adapted to ionospheric ray tracing. Specifically, it is used to move the chain of vertices along the generalized force so as to converge on the radio ray.
In the VPO method, the coordinates of the vertices are updated according to the equation of motion: where m is a 'mass' of the point in the configuration space and the two dots above r denote the second derivative with respect to fictitious 'time' t. Equation (18) is iterated using some numerical scheme and in every iteration, only the component of the velocity parallel to the force is kept, unless it is pointing in a direction opposite to the force, at which point it is zeroed. The VPO algorithm coupled with the velocity Verlet integrator [37] is presented in greater detail below. 1) Define parameters of the optimizer including the time step ∆t and fictitious mass m. In our calculations, we used the following parameter values: ∆t = 0.1 and m = 1 km −1 .
2) Set up initial position, r (0) , and velocity, v (0) = 0, of the path. Calculate the generalized force, F (0) . Use Eq. (6) for high rays and Eq. (10) for low rays. Start iterative procedure. 3) If force on each vertex of the path is lower than some predefined threshold value, F tol , then stop. Otherwise continue. 4) Given a path configuration r (I) , velocity v (I) and force F (I) at a current iteration I, update position of the path using the following formula: 5) Calculate the force F (I+1) corresponding to position r (I+1) . Use Eq. (6) for high rays and Eq. (10) for low rays. 6) Update velocity of the path using the following formula: 7) Revise velocity v (I+1) using the following formula: 8) Go to step 3.