ZERMELO NAVIGATION PROBLEMS ON SURFACES OF REVOLUTION AND GEOMETRIC OPTIMAL CONTROL

. In this article, the historical study from Carath´eodory-Zermelo about computing the quickest nautical path is generalized to Zermelo navigation problems on surfaces of revolution, in the frame of geometric optimal control. Using the Maximum Principle, we present two methods dedicated to analyzing the geodesic ﬂow and to compute the conjugate and cut loci. We apply these calculations to investigate case studies related to applications in hydrodynamics, space mechanics and geometry.


Introduction
A Zermelo navigation problem on a surface of revolution M is defined by the pair (g, F 0 ) where g is the metric on M induced by the Euclidean metric of R 3 and F 0 is a vector field on M called the current. Using the control formalism, see [16], the Zermelo navigation problem can be set as the time minimal transfer between two points q 0 , q 1 for the control system: u i (t)F i (q(t)), u = (u 1 , u 2 ), u = u 2 1 + u 2 2 ≤ 1, where q = (r, θ) are the polar coordinates for the metric g which takes the form g = dr 2 + m 2 (r) dθ 2 , see [6], the current F 0 being invariant by rotation and F 1 , F 2 form an orthonormal frame. The surface M is split into rectangles r 0 < r < r 1 with weak current if F 0 g < 1 or strong current if F 0 g > 1. Such a problem is a generalization of the historical problem of the quickest nautical path analyzed by Carathéodory and Zermelo [18,34], which have provided a complete study in the case of a linear current.
The first contribution of this article is to set the problem in the frame of geometric control, starting with the Maximum Principle [28] to introduce two different methods to analyze the geodesic flow. First of all, using the heading angle α of the ship, the system can be extended to an affine single input system: dq(t) dt = X(q(t)) + v(t)Y (q(t)), (1.1) withq = (r, θ, α) and v is the time derivative of α. Such a transformation being called a Carathéodory-Zermelo-Goh (CZG) transformation in this article. Using this approach, geodesics correspond to the so-called singular trajectories associated to (1.1) with v ∈ R, see [10] for this concept, and the geodesics can be classified into normal and abnormal geodesics, the second being on the zero level of the induced Hamiltonian. This approach allows to compute conjugate points along normal geodesics, where optimality is lost for the C 1 -topology on the set of geodesics that is in a conic neighborhood defined by the heading angle, this using the results and technical approach in [14], based on the computation of semi-normal forms. Furthermore, using similar techniques a first result of this article is to define and compute conjugate points along abnormal geodesics. More precisely, as already detected in the historical study, they correspond to a cusp singularity of the abnormal geodesics, when meeting the transition set F 0 g = 1 between the strong and weak currents. The second main technical contribution of this article is (following the approach used by I. Kupka in SRgeometry [25]) to analyze the set of geodesics using a one-dimensional mechanical system, with an extended potential V (r, p θ ) where the r-dynamics takes the form:

2)
p θ being the adjoint variable of θ which is constant using the Clairaut relation in the Hamiltonian frame. This leads to the analysis of the geodesic flow using an extension of the Morse-Reeb classification for 2d-Hamiltonian dynamics, see [6]. p. 21, and in particular to provide a stratification of the set of geodesics into r-periodic or r-aperiodic curves (i.e. periodic or not with respect to the variable r). Also in this frame, complicated dynamics can occur, in particular related to the existence of Reeb components in the geodesic foliation of M , see [20,22] for such occurrences in the modern study of foliation and dynamical systems. Finally, the third contribution of this article is to use the above techniques to analyze in details different case studies which form the core of this article, motivated by geometry and control theory. In each case, our aim is to compute the time optimal synthesis in the sense of [30,31] in an adapted rectangular domain R of the initial point q 0 . This means to compute in each case the cut locus Σ(q 0 ), where optimality is lost along geodesics initiating from q 0 , when restricted to R. In particular, three cases are studied in details. The first case is to analyze the historical example in our frame. In this case, every geodesic is r-aperiodic and the cut locus contains a single branch of the abnormal geodesic, terminating with a cusp singularity, when meeting the set F 0 g = 1. The second case deals with the Riemannian metric on a two-sphere of revolution, which appears in space mechanics and was analyzed in full details in [7]. Introducing a small current corresponds to the Finsler case analyzed in [5] for which the properties of conjugate and cut loci are well known, thanks to the continuity of the value function, and it is similar to the Riemannian case [19], p. 267. Moreover, the techniques of Poincaré and Myers can be used to compute the cut locus [27]. But we extend the analysis to the case of a strong current, and we show that the cut locus splits into two branches. This phenomenon is related to the optimality status of the abnormal geodesics [15] and to the shape of the small-time balls [12]. The third case study concerns the extension of the evolution of a passive tracer near a vortex and it was analyzed in details in [13]. This study motivated by applications in hydrodynamics [1] but also in relation with the N-body dynamics [26] is generalized and indicate the complexity of the geodesic dynamics, in relation with many Reeb components.
This article is organized as follows. In Section 2, we introduce the general concepts and definitions, and a large collection of case studies. In Section 3, we introduce the geometric tools of this article in relation with the geodesic curves, solution of the Maximum Principle. We define two different parameterizations of those curves. The first one is the CZG-extension related to conjugate points computations in both normal and abnormal cases and integration of the geodesic curves, using quadratures, in relation with Clairaut condition on surfaces of revolution. The second parameterization is described in Section 4 introducing the generalized potential and the generalized Morse-Reeb classification of the geodesics. In Section 5, which is the core of this article, we investigate in details the case studies. The final Section 6 is the conclusion which summarizes our contributions and proposes further studies.

Definitions and notations
Let M be a smooth surface of revolution with g the induced Riemannian metric and let T * M be the cotangent bundle endowed with the Liouville canonical form α = p dq. We recall for Theorem 4.1 that a Lagrangian manifold is a 2d-submanifold where dα is zero. We denote by q = (r, θ) the normal (polar) coordinates on the covering Riemannian manifold M c defined as 0 < r <r and θ ∈ R for a certainr, cf. [9]. Considering these coordinates, the metric takes the form g = dr 2 + m 2 (r) dθ 2 setting m(r) > 0 and the vector fields F 1 = ∂ ∂r and F 2 = 1 m(r) ∂ ∂θ define the canonical orthonormal frame. The lines r = constant are called the parallels and the lines θ = constant are called the meridians. A Zermelo navigation problem of revolution is defined by a triplet (M, g, F 0 ) where the vector field F 0 defining the current is invariant by θ-rotation, and we shall assume to cover the case studies that F 0 is oriented along the parallels only so that on M c it can be written F 0 = µ(r) ∂ ∂θ . If µ(r) is constant (resp. linear) this is called the constant (resp. linear ) current case. We define an adapted neighborhood of a point q 0 ∈ M c as a rectangle R = [r 1 , From the control point of view, the Zermelo navigation problem can be written in q-coordinates as: minimize the transfer time between two points (q 0 , q 1 ) for the system with admissible controls in the set of measurable functions defined on [0 , +∞) and valued in {u | u ≤ 1}.
The heading angle α of the ship in the canonical frame is defined by u 1 = sin α, u 2 = cos α, whenever u = 1, where according to Clairaut interpretation, α is the angle with respect to the parallel. Due to the symmetry of revolution, we can decompose the covering space into rectangles where we have either (at least in their interiors) a weak current if F 0 g < 1 or a strong current if F 0 g > 1, the transition between the two cases, when F 0 g = 1, being called the case of a moderate current. Furthermore, we can cover M c by adapted neighborhoods on which we can restrict the dynamics. For such an adapted neighborhood denoted R and a fixed point q 0 ∈ R, the navigation problem is called geodesically complete on R from q 0 if for any q 1 in R, there exists a geodesic contained in R joining q 0 to q 1 . Besides, we denote by q 1 → T (q 0 , q 1 ) the time minimal value function representing the minimal transfer time from q 0 to q 1 .

The historical example
A founding problem in classical calculus of variations is the problem called quickest nautical path introduced by Carathéodory and Zermelo [18,34] for a ship navigating on a river and aiming to reach the opposite shore in minimum time. Hence, M is the 2d-Euclidean space with metric g = dx 2 + dy 2 in the coordinates q = (x, y), y being the distance to the shore. To make a complete analysis, they have considered a linear current of the form F 0 = y ∂ ∂x . We shall refer to this case all along this article as the historical example. Using our notation to fix parallels and meridians, we must set x = θ, y = r, so that the ambient manifold is the Euclidean space with metric g = dr 2 + dθ 2 and F 0 = r ∂ ∂θ .

The vortex case
This case was analyzed in [13] and will be generalized in our study. The ambient space is the punctured Euclidean space and the vortex is localized at the origin and the ship is a passive tracer in hydrodynamics whose motion is described by: where k > 0 is the circulation parameter. The problem is written in polar coordinates x = r cos θ, y = r sin θ so that the Euclidean metric takes the form g = dr 2 + r 2 dθ 2 and the current becomes F 0 = k r 2 ∂ ∂θ . The ambient manifold is defined by r ≥ 0, F 0 having a pole at the vortex identified to r = 0.

The averaged Kepler case
The Riemannian problem related to the averaged Kepler problem in space mechanics (see [7]) can be extended to a metric on a two-sphere of revolution defined in normal coordinates by m 2 (r) = cos 2 r 1−λ cos 2 r where λ is a homotopic parameter, deforming the round sphere (for λ = 0) to the singular metric called the Grushin case (for λ = 1) and λ = 4/5 corresponds to the averaged Kepler case. For this case, we will consider a constant current on the covering space.

Ellipsoid of revolution
This standard problem of geometry is analyzed in [24]. The ellipsoid of revolution is generated by the curve: y = sin ϕ, z = ε cos ϕ where 0 < ε < 1 corresponds to the oblate (flattened) case while ε > 1 corresponds to the prolate (elongated) case. The metric takes the form g = F 1 (ϕ) dϕ 2 + F 2 (ϕ) dθ 2 , with F 1 (ϕ) = cos 2 ϕ + ε 2 sin 2 ϕ, F 2 = sin 2 ϕ. The metric can be set in the polar form using a quadrature. This defines a family of metrics on a two-sphere of revolution, depending upon ε.

The Serret-Andoyer case
The Serret-Andoyer metric studied in [11] corresponds to a representation of a mechanical pendulum. It is given in the normal form by taking m 2 (r) = (A cn 2 (αr, k) + B sn 2 (αr, k)) −1 , where cn and sn are Jacobi elliptic functions so that m(r) is periodic and moreover m(r) = m(−r). we have Remark 2.1. For the two last metrics, we can also define a Zermelo navigation problem by adding for instance a constant or a linear current on the covering space. We refer to [23], for the analysis of the ellipsoid of revolution in the Finslerian case. Note also that on a two-sphere of revolution a constant current corresponds to a linear rotation with the axis 0z.
3. The geometric tools from control theory and the Hamiltonian analysis

Generalities and the Maximum Principle
If not mentioned, all the objects are in a smooth category. Recall that we consider a Zermelo navigation problem of revolution determined by a triplet (M, g, F 0 ). The vector field defining the current can be taken in the form F 0 = µ(r) ∂ ∂θ which permit to cover all the cases defined above. The Zermelo navigation problem on the covering space M c consists to a time minimal transfer between two points (q 0 , q 1 ) for the control system: with u = (u 1 , u 2 ), u ≤ 1 and the set of admissible controls U is the set of measurable mappings defined on [0 , +∞) and valued in the domain U = {u | u ≤ 1}. Given q 0 ∈ M c and u ∈ U we denote by q(·, q 0 , u) the solution of (3.1) with q(0) = q 0 , and defined on a maximal interval J. We introduce the pseudo-Hamiltonian associated to the cost (extended) system H(z, u) = H 0 (z) + u 1 H 1 (z) + u 2 H 2 (z) + p 0 with z = (q, p), p = (p r , p θ ) being called the adjoint vector, H i (z) = p · F i (q) being, for i = 0, 1, 2, the Hamiltonian lift of the vector field F i , where · denotes the standard inner product, and finally p 0 is a constant representing the dual variable of the cost. We define the maximized (or true) Hamiltonian by the maximization condition: H(z, u).
Let q 0 ∈ M c and (t f , u) an optimal solution of the Zermelo navigation problem with q(·) = q(·, q 0 , u) the associated optimal trajectory. According to the Pontryagin's Maximum Principle [28], there exists an absolutely continuous function p and a scalar p 0 ≤ 0 such that the 4-uplet (q, p, u, p 0 ) satisfies, for almost every t ∈ [0, t f ], the following conditions:q , the pair (p(·), p 0 ) never vanishes.
(3.2) Definition 3.1. An extremal is a solution z(·) = (q(·), p(·)) of (3.2) and a projection of an extremal is called a geodesic. It is called strict if p is unique up to a factor, normal if p 0 = 0 and abnormal (or exceptional ) if p 0 = 0. In the normal case it is called hyperbolic (resp. elliptic) if p 0 < 0 (resp. p 0 > 0).
Since F 1 and F 2 form a frame, then the following will hold in our case.
Proposition 3.2. The trajectory q(·) is the projection of the extremal z(·) = (q(·), p(·)), solution of the systeṁ where the Hamiltonian vector field #-H is given by #-H = (∇ p H, −∇ q H). Besides, the optimal control is given in feedback form by 4) and the maximized Hamiltonian reads H(z) = H 0 (z) + p g + p 0 , setting p g = H 2 1 (z) + H 2 2 (z). Introducing the following definition, we can use results from [10]. Chap. 3 to characterize normal and abnormal extremals.
Definition 3.3. The fixed time extremity mapping is the map E q0,t f : u → q(t f , q 0 , u) and the extremity mapping is the map E q0 : u → q(·, q 0 , u), the set of inputs u being defined on a subdomain of L ∞ , endowed with the L ∞ -norm topology. The accessibility set in time t f , denoted A(q 0 , t f ), is the image of E q0,t f and the accessibility set A(q 0 ) = t f ≥0 A(q 0 , t f ) is the image of the extremity mapping.
Proposition 3.4. Take a reference extremal z(·) = (q(·), p(·)) on [0, t f ] where the corresponding control is given by (3.4). If we endow the set of controls valued in u = 1 with the L ∞ -norm topology we have: 1. In the normal case, u(·) is a singularity of the fixed time extremity mapping. 2. In the abnormal case, u(·) is a singularity of the extremity mapping.
Definition 3.5. Let t → q(t) be a response of (3.1). It is called regular if it is a one-to-one immersion. From the Maximum Principle, the geodesics can be parameterized by the initial heading angle α 0 and fixing q(0) = q 0 , we can define the exponential mapping as the map exp q0 : (α 0 , t) → Π(exp(t #-H)(q 0 , p 0 (α 0 ))) where Π is the q-projection: (q, p) → q. The cut point along a given geodesic is the first point where it ceases to be optimal and they will form the cut locus Σ(q 0 ). The separating line L(q 0 ) is the set of points where two minimizing geodesics starting from q 0 are intersecting.

Carathéodory-Zermelo-Goh transformation
In the historical example, Carathéodory-Zermelo integrated the dynamics using the heading angle α to parameterize the geodesics, see [16], p. 77. This corresponds to the Goh transformation in optimal control, see [10], p. 98. Next, we use this crucial point to make computations in the Lie algebraic frame and to relate the analysis of the navigation problem to the general study of [14].
Definition 3.6. Let us consider the control system (3.1) with the restriction u = 1. Then, we can set u = (cos α, sin α), α being the heading angle of the ship. We denote byq = (q, α) the extended state and we introduce the vector fields X(q) = F 0 (q) + cos α F 1 (q) + sin α F 2 (q) and Y (q) = ∂ ∂α . This leads us to prolongate (3.1) into the single-input affine system: and the derivative of the heading angle v(t) = dα dt (t) is called the accessory control. Denotingz = (q,p) = ((q, α), (p, p α )), this leads to define the extended pseudo-Hamiltoniañ Parameterization of the geodesic curves in this extension. From [10], Chapter 6, in this representation, geodesic curves extend to singular trajectories of (3.5), where the accessory control v belongs to the whole set R.
Definition 3.7. The Lie bracket of two vector fields U , V is computed with the convention and is related to the Poisson bracket where H U and H V are the Hamiltonian lifts of U and V .
Lemma 3.8. One has the following: Assume that D(q) is not vanishing, then the singular control v(·) associated to the geodesics is given by the feedback and the geodesics extend to smooth solutions of For this dynamics, one has: hyperbolic geodesics are in DD > 0, elliptic geodesics are in DD < 0, abnormal (or exceptional) geodesics are located in D = 0.
Proof. From [10], Section 3.4, taking a singular control-trajectory pair (q(·), v(·)), denotingz = (q,p) one has and this leads to Hence, sincep ∈ R 3 \ {0},p can be eliminated using (3.9) and the singular control reads as the singular feedback Since D(q) is never vanishing, one has that Y (q), [Y, X](q) are linearly independent and the adjoint vectorp is unique up to a factor. Thus, every geodesic is strict. Moreover, the following condition called the strict generalized Legendre-Clebsch condition is satisfied: along any geodesic extension. This amounts to the strict Legendre-Clebsch condition: ∂ 2 H ∂α 2 = 0. The classification of the geodesics (in the strict case) into hyperbolic, elliptic and abnormal geodesics defines for the extension, the set DD" > 0, DD < 0 and D = 0.
Proposition 3.10. The dynamics of the geodesics in polar coordinates is given by: Hence, the Hamiltonian geodesic flow is Liouville integrable [2], p. 269 with two involutive first integrals H and p θ .
Proof. Computing with polar coordinatesq = (r, θ, α) one has: and Lie brackets computations give: This leads to: Then, (3.10) follows from (3.7). On the other hand, the pseudo-Hamiltonian in theq-representation takes the form H = p r cos α + p θ µ(r) + sin α m(r) + p 0 and from the maximization condition with v ∈ R, one has ∂H ∂α = 0. This gives the Clairaut relation p r sin α = p θ m(r) cos α. Assuming p θ = 0, then we have sin α = 0 and then plugging such p r into H, we obtain the relation relating r to α and this leads to the parameterization of the geodesics solution of (3.10). The case p θ = 0 is clear since in this case, sin α = 0 and soα = 0.
Application to the historical example.

Moderate current domain and regularity properties of the geodesic curves
By definition in polar coordinates, the moderate current domain F 0 (q) g = 1 is given by the relation: and so one can compute the decomposition of the domain into weak and strong current domains using the graph of the function f (r) = µ 2 (r) m 2 (r).
Definition 3.12. Weak and strong current domains are defined respectively by: f (r) < 1 and f (r) > 1. A moderate parallel r 0 solution of f (r) = 1 is called regular if the intersection of the line L(r) = 1 with the graph of f is transversal at r 0 , i.e. the root r 0 is simple, so that at r 0 the moderate parallel defines a non-empty transition between weak and strong current domains.
Lemma 3.13. Letσ = (q, α) be a reference geodesic solution of (3.10) on [0, t f ] and intersecting at time t f the moderate current domain. Then: 1.σ is either hyperbolic or abnormal.
Proof. Elliptic geodesics are contained in the strong current domain and only hyperbolic and abnormal geodesics can intersect the moderate current domain F 0 (q) g = f (r) = 1. Take a point q = (r, θ) in the moderate current domain. Then, there exist a heading angle α such that the following holds: From (3.10), this leads to the following relations for the dynamics: Hence we obtain: Sinceṙ =θ = 0, the geodesic inersecting the moderate current domain with the heading angle π 2 + kπ is abnormal. In hyperbolic case, one hasq(t f ) = 0. The lemma is proved.

Computation of conjugate points using the CZG-transformation
The aim of this section is to provide algorithms to compute conjugate points. The regular case is covered in [14] whereas the recent article [15] analyzes in a general framework conjugate points along non-immersed abnormal geodesics in Zermelo navigation problems. Note that in both cases, the integrability property of the geodesic flow is not required, but in fine the semi-normal form is integrable cf. [14].
First of all, we shall analyze the case whereσ(t) is a regular geodesic. It is based on [14]. [14] to determine conjugate points in the regular case Semi-normal forms. We assume that the reference geodesic curve t →σ(t) on [0, t f ] is a one-to-one immersion. In this case, we can choose coordinates and a feedback so thatσ : t → (t, 0, 0) and it can be taken as the response of the control v = 0. Further normalizations are obtained in the jet-spaces of (X, Y ) in the neighborhood ofσ.

A brief recap of
Normal case. We can choose coordinatesq = (q 1 , q 2 , q 3 ) such that the system takes the form: where ε 1 (q) can be neglected in the analysis and with a 33 < 0 (resp. a 33 > 0) in the hyperbolic (resp. elliptic) case.
Abnormal case. We can choose coordinatesq = (q 1 , q 2 , q 3 ) such that the system takes the form: where ε 2 (q) can be neglected in the analysis and where a is strictly positive. See [14] for details about the computations and the descriptions of the mappings q → ε 1 (q), ε 2 (q). In both cases, since q 1 (t) = t along the reference geodesic σ, one can replace q 1 by t in the semi-normal form (restricting our analysis to a conic neighborhood) and this allows to evaluate the accessibility set and its boundary. Hence, computing conjugate points to deduce the optimality status of the reference geodesic.
with X s = X + v s Y and v s given by (3.7). A Jacobi field J(t) is a solution of (3.16) which is said semi-vertical at t = 0 if J(0) ∈ RY (σ(0)).
1. In the normal case, if t c is a conjugate time, then, there exists a non-trivial Jacobi field J semi-vertical at t = 0 such that Let t 1c be the first conjugate time, then, if the geodesic is hyperbolic (resp. elliptic) it is time minimizing (resp. maximizing) with respect to all geodesics in a conic neighborhood of the reference geodesic σ, up to time t 1c . 2. In the abnormal case, the reference geodesic is both minimizing and maximizing in a conic neighborhood ofσ.
Remark 3.17. The result is clear in the abnormal case due to (3.15), since q 3 (·) is strictly positive, unless the geodesic curve is the reference geodesic. It was already observed by Carathéodory and Zermelo, see [18] where the abnormal geodesic are called "limit curves".
The concept of generalized curvature using CZG-transformation in the normal case. It was defined in [10], p. 163 and it can be used in our Zermelo problem. The Jacobi equation alongσ : t → (t, 0, 0) takes the form: where a(t), b(t), c(t) denote in short the coefficients of a 33 , a 23 + a 32 , a 22 in formula (3.14). The existence of conjugate time t 1c means that there exists a non-trivial solution of (3.17) satisfying δy(0) = δy(t 1c ) = 0. It can be written in the normal formẍ and K(t) is the generalized curvature defined by Remark 3.18. Note that the generalized curvature depends upon the reference geodesic parameterization contrary to the Gauss curvature in Riemannian geometry. In Section 5, the averaged Kepler case will exhibit some geodesic images such that each of them is parameterized by an hyperbolic trajectory and an elliptic trajectory (for instance the equator).

Cusp singularity associated to a conjugate point in the non-regular (abnormal) case
Next, we use [15] to describe the cusp singularity of the abnormal geodesics when meeting the transition between the strong and weak current domains. It is based on [32] in the algebraic case and [3] in the C ∞ -case. Let γ : t → (x(t), y(t)) be a planar smooth curve. We recall that if a point γ(t cusp ) is a cusp point of order (p, q), 2 ≤ p ≤ q, then γ (p) (t cusp ) and γ (q) (t cusp ) are independent. It is called an ordinary (or a semicubical ) cusp point if it is a cusp point of order (p, q) with p = 2 and q = 3. From [32], p. 56, an algebraic model in R[x, y] at an ordinary cusp γ(t cusp ) is given by the equation x 3 − y 2 = 0. Moreover it is the transition between a R-node solution of the equation x 3 − x 2 + y 2 = 0, where the origin is a double point with two distinct tangents at 0: x ± y = 0 and a C−node solution of x 3 + x 2 + y 2 = 0 with two complex tangents at 0 given by x ± iy = 0 and with two distinct components x = y = 0 and a smooth real branch (see also [15] for more details and Fig. 2 for an illustration).
A neat description from singularity theory suitable in our analysis is given by [3], p. 65 and is associated to a typical perestroika of a plane curve depending on a parameter and having a semicubical cusp point for some value of the parameters, where the curves sweep an umbrella while their inflectional tangents sweep another umbrella surface. Below, we give some results to describe the properties of semicubical cusp points in relation with abnormal geodesics as well as the optimality of abnormal and hyperbolic geodesics in the neighborhood of cusp points. Proposition 3.19. Consider a reference abnormal geodesicσ a = (σ a , α a ) defined on [t 1 , t 2 ], t 1 < 0 < t 2 so that at t = 0, σ a meets a regular moderate parallel r 0 at q 1 = σ a (0). Denote σ(t 1 ) = q 0 and assume thatα a (0) = 0. Then, there exists a neighborhood V of q 1 such that for every geodesic starting from q 0 one has: 1. The abnormal geodesic σ a meets the boundary at a semicubical cusp with a vertical tangent. 2. Hyperbolic geodesics are self-intersecting curves corresponding to a R-node. 3. elliptic geodesics exist only in the strong current domain and correspond to a C-node.
Proof. The proof follows from [15] and is illustrated by Figure 3. 1. The abnormal arc is optimal from q 0 up to the cusp point q 1 included which corresponds to a conjugate point in the abnormal case. 2. Self-intersecting geodesics starting from q 0 are optimal up to their intersection point q 2 with the abnormal arc, the point q 2 being excluded. Hence, q 2 corresponds to a cut point. 3. The time minimal value function T : q 2 → T (q 0 , q 2 ) is discontinuous for each q 2 = q 1 on the abnormal geodesic.
Proof. The proof follows from [15] and is illustrated by Figure 3. See also the cut loci of Figure 4 for an illustration.
Corollary 3.21. Under the previous assumptions, the abnormal arc σ a is optimal in the whole domain of strong current up to r 0 , restricting to a conic neighborhood.
Proof. Letσ a : t →σ a (t) = (σ a (t), α a (t)), t ∈ [t 1 , 0], t 1 < 0 (where −t 1 not necessarily small) be the reference abnormal geodesic. From the previous theorem, the abnormal arc starting from q 2 = σ a (t 2 ) with t 1 < t 2 < 0, t 2 small enough, is optimal from q 2 to q 1 . Moreover, the abnormal arc from q 0 = σ(t 1 ) corresponds to an immersed curve and thus from [14], it is optimal on the whole strong current domain, restricting to a conic neighborhood. Hence, concatenating the two arcs allows us to conclude.   . Time minimal synthesis in two different adapted rectangles containing the limit loop l loop (q 0 ). The second rectangle being delimited by the limit loop to emphasize the dependence of the cut locus to the considered adapted rectangle. In both cases, the cut locus (in black) contains the whole abnormal arc terminating at the cusp singularity. In gray is represented the strong current domain. The red curves are hyperbolic geodesics while the green represent abnormals.

Time minimal synthesis in the historical example in an adapted rectangle
In this section we use the global parameterization of the geodesic curves given by Proposition 3.10 to compute the time minimal synthesis in an adapted rectangle containing the limit loop starting from q 0 . The limit loop starting from q 0 , see Figure 4 for an illustration, is the geodesic starting from q 0 which returns to q 0 in minimum (positive) time. In the historical example, starting from q 0 in the strong current domain, there exists a limit loop denoted l loop (q 0 ) with return time denoted t 0 . The time minimal synthesis is represented on Figure 4 and the main properties are the following. Proposition 3.22. In the historical example, taking a point q 0 in the strong current domain, then we can choose an adapted rectangle R containing the limit loop l loop (q 0 ) such that: In the domain, the cut locus contains the whole branch of the abnormal geodesic arc q 0 q 1 and is the union of the separating line L(q 0 ) where abnormal and normal minimizing arcs intersect with unequal time and the terminating point q 1 which is a conjugate point of the non-immersed abnormal arc.

Mechanical system and generalized Morse-Reeb classification
4.1. Mechanical representation in the Riemannian averaged Kepler case [7,8] In this section, we give a brief description of the mecanical representation in the Riemannian case, following the example of the averaged Kepler metric. In the Riemannian setting, roughly speaking, the time and energy minimization problems are equivalent. Classically, the (maximized) Hamiltonian associated to the energy minimization problem is preferred. This Hamiltonian is given by where the ambient manifold M is the (compact) 2-dimensional sphere of revolution and the metric g on the covering space M c is given by Moreover, since θ is a cyclic variable, then G = p θ is an additional first integral given by the Clairaut relation. Trajectories of #-H split into three cases: the meridians defined by θ = constant with p θ = 0, the equator identified to r = 0 with p θ = m(r), and remaining geodesics with |p θ | ∈ (0 , m(r)) formed by r-periodic oscillating solutions of the mechanical system, defined by the so-called characteristic equation the term V (r, p θ ) = p 2 θ m 2 (r) representing the potential. A further integration is necessary in order to recover the θ-variable using the Hamiltonian dynamics. Parameterizing by r on each ascending branch of the characteristic equation, we have: This allows to compute the variation denoted ∆θ/2 of the angle θ starting from the equator and on the ascending branch, the total variation to return to the equator being ∆θ. Note that in the limit case of the equator solution, the rotation is stationary since r is constant. This gives the complete description of the Liouville torus T 2 defined in Theorem 4.1 hereinafter. Note that the geodesics split into periodic orbits if ∆θ/2π is rational and dense orbits in the Liouville Torus if ∆θ/2π is irrational. 1. T ξ is a smooth Lagrangian manifold invariant by the flow of #-H and #-G. 2. If T ξ is connected and compact, then T ξ is diffeomorphic to the 2-dimensional torus T 2 (called Liouville torus). 3. The Liouville foliation is trivial, that is, there exists a neighborhood U of T ξ such that U is the direct product of the torus T 2 and the disk D 2 . 4. In the neighborhood U there exists symplectic coordinates (s, ϕ) (also called action-angle variables) such that the dynamics #-H can be written in the form dsi dt = 0, dϕi dt = α i (s 1 , s 2 ), i = 1, 2, and for which the motion is quasi-periodic.
Next, we present a generalization of the previous mechanical representation in the Zermelo case.

Mechanical representation in the Zermelo case
In this case, since we have a current, there is no longer equivalence between the time and energy minimization problems. So, we consider the Hamiltonian associated to the time minimization problem given by The following theorem sets the framework for the classification of the geodesic flow. 1. The evolution of the system in the (r, p r )-space is described by the Hamiltonian dynamics dr dt = p r p g , dp r dt = −p θ µ (r) + p 2 θ m (r) m 3 (r) p g . (4.1) 2. Restricted dynamics (4.1) can be integrated using the mechanical system representation where the generalized potential is given by with ε = p 0 normalized to −1, 0, +1 respectively in the hyperbolic, abnormal and elliptic cases. 3. Moreover, the following relation holds Proof. Item 1 is a direct compution of equation (3.2) considering the Hamiltonian defined above. Moreover, since the Hamiltonian H = p g + p θ µ(r) + p 0 = 0 and p g = p 2 r + p 2 θ m 2 (r) 1/2 , we have: which gives item 3. Besides, from the restricted Hamiltonian system in (r, p r ) we have: which gives item 2. Orbits are classified according to the following definitions.
Definition 4.4. We consider the hyperbolic case ε = −1. A pair (r * , p * θ ) is called an equator if (r * , 0) is an equilibrium point of the restricted Hamiltonian dynamics, parameterized by p * θ . It is called L-elliptic if the linearized dynamics is with spectrum {±iα, α = 0}, L-hyperbolic if the spectrum is of the form {±λ, λ ∈ R \ 0} and L-parabolic if the spectrum is zero. The L-elliptic and L-hyperbolic situations correspond respectively to a stable case associated to a minimum of the potential and to an unstable case associated to a maximum. An equator defines a stationary rotation in the (r, θ)-space, it is called a positive (resp. negative) rotation if θ is rotating with a positive (resp. negative) frequency. A separatrix geodesic parameterized by p * θ is a geodesic (r(t), θ(t)) such that r(t) → r * as t → ∞ where (r * , p * θ ) is an equator. We give next some characterizations of equators and separatrices.  2. An equator (r * , p * θ ) is L-hyperbolic (resp. L-elliptic) if and only if:
These two equations then lead us to Proposition then follows.
Definition 4.6. Let U be an adapted rectangular neighborhood of q 0 . Geodesics starting from q 0 decompose into starting ascending or descending r-branches. At the limit tangential case, monotonicity is given by positive or negative acceleration d 2 r dt 2 (0). The first return to the equator (resp. the meridian) associated to a geodesic is the first point such that the geodesic intersects again the equator (resp. the meridian).
Proposition 4.7. Let U be a rectangular adapted neighborhood of q 0 = (r 0 , θ 0 ). Level sets of the Hamiltonian in the GMR-classification split into compact levels corresponding to r-periodic geodesics and non-compact level sets corresponding to r-aperiodic geodesics, when the dynamics is restricted to the neighborhood U . If (r * , p * θ ) is a L-elliptic equator, then locally the Liouville foliation by Liouville torus is preserved.
Remark 4.8. Let q 0 be a fixed initial condition, then using the GMR-classification for each adapted rectangular neighborhood of q 0 we can stratify the set of geodesics emanating from q 0 into micro-local conic sectors corresponding to compact and non-compact geodesics.
Remark 4.9. The decomposition depends upon the adapted rectangular neighborhood and can be described using the generalized potential restricted to the domain. we can easily have situations with two compact sectors separated by a singular level with a separatrix geodesic, or an equator for which when restricted to the domain, the singular level separates compact and non-compact orbits. This property can be illustrated for example in the Serret-Andoyer case, described in details in [11]. In this case the meridian can be identified to r = 0 and starting from the meridian, geodesics split into r-periodic curves and r-aperiodic curves with a limit curve corresponding to a separatrix. Only r-periodic curves have conjugate points. They are a representation of the standard pendulum equation where the oscillating solutions correspond to r-periodic solution while rotating solutions correspond to r-aperiodic solutions. But if they are interpreted on the cylinder, both types of solutions are periodic, oscillating trajectories are homotopic to a point, but not the rotating ones.

The averaged Kepler case
The aim of this section is to analyze the Zermelo navigation problem associated to the Riemannian case studied in [7,8]. We first recap the properties of the Zermelo problem associated to the Riemannian metric, when the current is zero and then extend the results to the Zermelo case with strong current.

Riemannian case
Recall that the metric takes the form g = dϕ 2 + m 2 λ (ϕ) dθ 2 , with where ϕ represents the angle on the two-sphere of revolution, where ϕ = 0 (resp. ϕ = π) corresponds to the north (resp. south) pole and λ ∈ [0, 1] is an homotopic parameter, λ = 0 being the round sphere, λ = 4/5 is associated to Kepler orbits transfer and λ = 1 is the so-called Grushin case with a singularity at the equator. The Gauss curvature is given by and is strictly negative in the limit case λ = 1. The only equator solution corresponds to ϕ = π/2 and we introduce r = π/2 − ϕ to normalize this equator to zero. The metric is then written g = dr 2 + m 2 (r) dθ 2 , where we set m(r) = m λ (π/2 − r) and it is symmetric with respect to the equator, that is m(r) = m(−r), which is crucial for the explicit computations of the conjugate and cut loci, following [9]. Using the Hamiltonian formalism we associate to the metric the Hamiltonian Fixing the parameterization to the arc-length amounts to set p 0 = −1. So that the characteristic equation takes the form: A geodesic is either a meridian if p θ = 0, the equator if p θ = m(r) and any solution is such that r is periodic and oscillates between −r + , r + and is entirely determined by a branch of the characteristic equation evaluated on the quarter of period T /4, r(t) being restricted to [0, r + ], where r + is the positive root of the equation V = 1, the period being given by the integral which depends upon p θ . By symmetry with respect to the equator it can be supposed non-negative and belonging to (0, m(r(0))). To make the analysis we introduce the application called the period mapping associated to the first return to the equator and defined by: The geodesic flow is integrated by quadrature using the characteristic equation and the transcendence of the solutions is basically determined by the transcendence of the period mapping. In this context, this case study is rather straightforward, since only elementary functions are necessary to parameterize the solutions. In particular, it is related to the historical example, replacing the ambient Euclidean space by a two-sphere. To integrate, we can assume that r(0) = 0 and θ(0) = 0, since every oscillating trajectory is such that r is intersecting the equator and we use: One gets that where n ∈ N counts the number of intersections with the equator and by symmetry, we can assume that the number of intersections is odd. The function ∆θ for p θ ∈ (0, m(0)) is the first return mapping to the equator introduced in Definition 4.6. The following proposition coming from [9] is crucial in the optimality analysis.
Proposition 5.1. Restricting the initial point to q 0 = (0, 0) and assuming that the first return mapping to the equator is tame, that is monotone non-increasing, then the first conjugate time is given by the equation where θ is parameterized by r according to and the first conjugate time being between T /2 and T /2 + T /4.
Integration of the geodesics. We take dr dt and we denote by Z + and Z − the roots of 1 + p 2 θ (λ − 1) = Z 2 (1 + λp 2 θ ), where Z = sin r and the period reads .
Normalizing the amplitude of the oscillation by Z = Z + Y one has . Proposition 5.2. The period is given by Moreover we have arcsin Y (t) = 1 + λp 2 θ t. This defines the re-normalized time s = √ 1 + λp θ t and the θ-variable is integrated using A quadrature gives the following. The θ-dynamics is given by This leads to a complete parameterization of the geodesics. Moreover, we can compute explicitly the conjugate and cut loci in the tame case according to the following propositions.
Proposition 5.3. In the tame case, the cut locus of a point on the equator is a sub-arc of the equator and the injectivity radius is the distance to the cusp extremity of the conjugate locus on the equator.
Proposition 5.4. Assume that the problem is tame and moreover, suppose that the first return mapping ∆θ is such that ∆θ < 0 < ∆θ on (0, m((r(0))), then: 1. The cut locus of a point not a pole is a segment of the antipodal parallel. 2. For a point not a pole, the conjugate locus has exactly four cusp points.
Computing, we have.
This can be applied to our case for λ ∈ (0 , 1). Note also that the conjugate locus of the equator is a standard astroid with four cusps. The limit Grushin case λ = 1 can be analyzed similarly, except that the equator is not a geodesic and the injectivity radius is zero. This gives a complete analysis of the Riemannian case and this leads to the following analysis.

Transition from the Riemannian case to the Zermelo case with a constant current
Recall that the problem with constant current is given on the covering space by the pair where v is a non-zero constant. Depending on the current at the initial point q 0 = (r 0 , θ 0 ), we say that we are in the weak (current) case if sin 2 r 0 < 1 v 2 +λ , strong case if sin 2 r 0 > 1 v 2 +λ and moderate case if sin 2 r 0 = 1 v 2 +λ . In the case where v 2 + λ < 1, the current will be weak on the whole domain. So we shall assume: v 2 + λ > 1. The following is a crucial geometric property.
Proposition 5.6. On the two-sphere of revolution embedded in R 3 , the vector field F 0 defines a linear vector field, tangent to the sphere, and it corresponds to a uniform rotation whose axis is the axis of revolution. For the metric the equator solution is also a stationary rotation since dθ dt is constant so that the effect of the current can be added to this rotation.
Integration of the geodesics. From the previous proposition, the integration follows from the Riemannian case. Introducing the generalized potential, recall that the r-dynamics is given by: where ε = p 0 < 0, = 0, > 0 correspond respectively to the hyperbolic, abnormal and elliptic cases. Taking the ascending branch starting from the equator r 0 = π/2, we have dr dt = p 2 θ (1 − λ sin 2 r) sin 2 r(ε + p θ v) 2 1/2 , Since H = 0, p g = −(p θ v + ε), then, using a time reparameterization, one gets: , which is like the r-dynamics in the Riemannian case, with the addition of v. Then, we can determine the first return mapping to the equator r 0 = π/2: where r + is the maximum of r(t).
Again the geodesic curves are symmetric with respect to the equator, the cone of admissible direction being symmetric with respect to the equator. This leads to the following stratification of the set of geodesics, using the variable p θ instead of the heading angle in the historical example.
Proposition 5.7. Assume that λ = 4/5 and v = 0.9, then starting from the equator and considering only the ascending branch, geodesics split into: Abnormal given by p a θ = −1/v; Hyperbolic geodesics parameterized by p θ ∈ (p a θ , m(r 0 )); Elliptic geodesics parameterized by p θ ∈ (−m(r 0 ), p a θ ). Moreover, in the hyperbolic case, the set of geodesics can be stratified in four different classes (see Fig. 5 for illustration): The equator which corresponds to r = π/2, p r = 0 and p θ = m(r 0 ). The two pseudo-meridians (ascending and descending ones) which correspond, on the covering space, to non-compact case where p θ = 0. Generic r-periodic orbits which split in two different families namely orbits without self-intersections, parameterized by p θ ∈ (0, m(r 0 )) and orbits with self-intersections, parameterized by p θ ∈ (p a θ , 0) and ±p r (0) corresponding to the symmetric orbits.
Remark 5.8. The other geodesics in the flow are obtained by a symmetry which respect to the equator.
The cut locus in this case will split into two branches. The first branch is associated as in the historical example to the cusp singularity of the abnormal directions, which are symmetric with respect to the equator (see Fig. 6 for an illustration). The second branch of the cut locus is the persistence of the segment formed by the equator and related to the tame behavior of the first return mapping corresponding to non self-intersecting geodesics. The conjugate points can be numerically evaluated. They exist for different types of geodesics but occur after the intersection of the geodesics with the equator. See Figures 7 to 9 for illustrations. Theorem 5.9. Assume that the equator r 0 = π/2 is in the strong current domain. Then the cut locus has two branches, the first branch being formed by the abnormal curves occurring in the neighborhood of the cusp point and associated to self-intersecting geodesics and the second branch being a segment of the equator, starting by a cusp point of the conjugate locus and associated to non self-intersecting geodesics.

Complexity of the Hamiltonian dynamics in the generalized vortex case
A preliminary study of a navigation problem with a vortex localized at the origin is studied in [13], but the aim of this section is to generalize this situation to get more intricate dynamics, in relation with the N-body problem. A first step is to generalize the existence theorem to get a geodesically complete framework.

Existence of optimal solutions
We consider the case where the Zermelo navigation problem is defined by: In our preliminary study, we have µ(r) = k r 2 but in the general case we shall assume that µ(r) has a pole of order β ∈ (1, +∞) at zero, so that near zero, we can take the approximation µ(r) ∼ 1 r β and moreover we assume that µ(r) → 0 as r → +∞. We shall generalize the argument of [13] and relate the proof to existence of minimizing solutions avoiding collisions in the N-body problem [21]. Also, we point the relation between extending the solutions beyond the vortex and the Levi-Civita regularization for double collision [29].
Theorem 5.10. Take q 0 , q 1 in the punctured plane R 2 \{0}, then there exists a time minimizing trajectory to transfer q 0 to q 1 . Moreover, q 0 = (r 0 , θ 0 ) can be transferred to the origin in minimum time t min = r 0 . Hence, we can extend the geodesic flow using a Levi-Civita type regularization beyond the collision with the pole by reversing the time parameterization when crossing the vortex.
Proof. The geodesic dynamics in polar coordinates reads As in [13], to prove the existence about minimizing trajectories it is sufficient to prove that the minimizing trajectories avoid the collision. Using the expansion near the pole and comparing the time to make a rotation around the pole on a circle of radius r denoted T θ (r) and the minimal time to reach a circle of radius ε denoted  (Fig. 7 top), Finslerian (Fig. 7 bottom) and Zermelian (top) cases. In red (resp. in blue) is the hyperbolic (resp. elliptic) geodesics. Conjugate locus for the hyperbolic (resp. elliptic) geodesics is represented in black (resp. in dashed blue). The gray sector represents the strong current domain. Note that the conjugate locus is a closed curve on the manifold, cf. the bottom figure. T r (ε) a direct computation gives Hence, the argument of [13] to replace a trajectory reaching a circle with small radius ε by the trajectory making a rotation around the pole (see Fig. 10) is still valid and the existence result follows. Clearly from the equations the time to reach the pole from q 0 is obtained for p θ = 0 and is given by r 0 . Following the Levi-Civita regularization we reverse the geodesics orientations when crossing the vortex. It amounts in particular to replace µ(r) by −µ(r) and p θ by −p θ in the geodesics equations.
Remark 5.11. To replace the pole of the vortex potential by β > 1 in the general vortex case is similar to the modification done by Poincaré in the Keplerian potential where he replaced the pole by β ≥ 2 (β ∈ N) in order to avoid collisions. In our case the bound of the pole is given by the conditions T θ (r) < T r (ε) i.e.

The single vortex case in hydrodynamics
On the punctured plane we consider the case of an Euclidean metric in polar coordinates with The generalized potential is given by The geodesic curves can be classified using the potential and the main features are described hereinafter. See [13] for more details and see also Figures 11 and 12 for illustrations. We first introduce the following novel definition.
Definition 5.12. A Reeb component is a foliation invariant by θ-rotation generated by a separatrix geodesic (r(t), θ(t)) such that r(t) converges when t → ±∞ to two different equators, with different orientations.
Then, in the single vortex case, whe have the following results.  . The geodesics below the separatrix (blue), i.e. parameterized by p θ < p * θ (p * θ parameterizing the separatrix), converge towards the vortex and those above the separatrix (red) go to infinity. However, in the neighborhood of the separatrix the geodesics come from the vortex. (Right) Illustration of the two pseudo-meridians represented in red and blue. Note that both orbits do not coincide because of the drift. The dotted lines represent the semi-orbits computed backward in time.
Proposition 5.13. Considering the single vortex case, for a given q 0 ∈ R * \ {0}, we have: The domain of strong current is near the vortex and limited by the circle of radius r = k of moderate current. The only equator solution is in the domain of weak current and is defined by the circle with radius r * = 2k. There exists a unique separatrix emanating from the vortex and converging to the equator with dθ dt = 0 on the circle of radius 2k/ √ 3. This separatrix forms a Reeb component in the interior of the disk delimited by the equator, whose foliation as a singularity at the vortex.
At the exterior of the circle of radius r * there exists a unique separatrix emanating from the equator and converging to the infinity. The separatrices split the geodesic flow into trajectories that converge towards the vortex and those that go to infinity. There exists two pseudo-meridians converging with maximal radial speed either towards the vortex or to the infinity.

Generalized single vortex case
In this section, we describe a generalization of the situation analyzed in Section 5.2.2 in which we have several equators and separatrices in relation with similar studied in dynamical systems, see for instance [22]. We assume m = 1 and µ(r) rational but complex dynamics can be obtained with this restriction. We consider a system with a current of the form: λ, β ∈ R \ {0}, so that the generalized potential is now given by Moreover we assume that λ 2 > −3β and β < 0. To integrate the flow, as in the Kepler case, we first use the characteristic equation to integrate the r-dynamics. This leads us to dr dt = 1 − V ε (r, p θ ) = (εr 3 + p θ λr + β) 2 − p θ r 4 εr 3 + p θ λr + β) 2 .
Then, we use an additional quadrature to integrate the θ-dynamics. Indeed, the θ dynamics given by dθ dt = ∂H ∂p θ can be re-parameterized by Finally, the geodesic flow can be classified using the potential and the Clairaut parameter p θ . At the end, we have the following.

2)
and respectively associated to the Clairaut constants: and .      From the potential plot, we deduce that the trajectories starting close to r 2 remain contained in an annulus (see the midlle figure for an example); those coming from the vortex return to it and those starting sufficiently far from the vortex go to infinity (cf. right figure).
Equator r 1 forms a Reeb component delimited by the vortex and r 1 , as in the simple vortex case (see Fig. 13). while equator r 3 is an homoclinic trajectory (see Fig. 15). When t → +∞, the micro-local classification gives: • Equators r 1 and r 3 have two sectors formed by geodesics converging either to the vortex when t → +∞ or to infinity. The two sectors being delimited by a pseudo-meridian (see Fig. 16). • Equator r 2 has three sectors, two of them being as in the previous item, the third one corresponds to an elliptic sector (see Fig. 17).

Algorithm
We can deduce from the previous studies the method of analysis to handle a general case, and we proceed as follows. In the normal coordinates (r, θ) on the covering manifold M c we have r ∈ (0, R) and we can decompose the domain into disks c i < r < c i+1 with alternatively weak and strong currents. We compute the equators solutions listed as 0 < r 1 < r 2 < · · · < r p < R and they can be classified according to their optimality status into L-hyperbolic or L-elliptic equators. Then, taking a point q 0 , we can parameterize the geodesics using the mechanical representation with the generalized potential using improper integrals. This allows to construct the time minimal synthesis fixing an adapted neighborhood, using the first return mappings to equators and meridians combined with con jugate point analysis, in both normal and abnormal cases. Note that in the strong current domain the size of the adapted neighborhood is defined by the limit loop of the self-intersecting geodesics related to the abnormal geodesics. This can be extended to a larger domain by gluing different adapted neighborhoods, see [30,31] for such a procedure for general control problems in the plane.

The gluing process
In the previous section we describe the method to obtain global syntheses, by gluing together the syntheses using different adapted neighborhoods but intricate situations can be constructed starting from case studies by gluing such cases using the normal coordinates (r, θ). Indeed, each case is described by a pair of covariant (µ i (r), m i (r)), parameterizing respectively the current and the metric. They can be glued together in the C ∞category using bump functions. For instance the vortex case with Euclidean metric can be glued to the averaged Kepler case, to represent the Zermelo navigation problem of a passive tracer, swallowed by the vortex to enter into a Kepler domain, to visit an equator solution, with non zero curvature and with different types of geodesics.

Conclusion
The main contribution of this article is starting from the historical navigation problem studied by Zermelo and Carathéodory to introduce general tools from geometric control to analyze navigation problems on surfaces of revolution, based on the inspection of the geodesic curves. Two complementary techniques have been introduced. First of all, the system can be extended to a single-input control system using the Goh-transformation and this extension is suitable to parameterize the geodesics by quadrature, using the heading angle of the ship and to compute conjugate points in normal and abnormal cases. Second, we observe that the dynamics can be integrated using a generalization of the methods used for 2d-integrable mechanical systems as a step to compute actionangle variables. This leads to introduce in our frame, an extended potential. The application is to analyze the geodesics using an extension of the Morse-Reeb classification [6].
These techniques are used to study three case studies which are the core of the article. The first one is to recover in a neat geometric context the seminal results in the historical example. The second study is related to celestial and space mechanics and is analyzed by homotopy starting from the weak current to strong current, and deforming metrics on two-spheres of revolution. The interest of this case is to prove that in general the cut locus splits into two branches. The first branch being formed by the standard Riemann-Finsler case, thanks to continuity of the value function. But a second branch occurs in the strong current domain and is associated to the semicubical singularity of the abnormal direction, with neighboring self-intersecting normal extremals. Cut point occurs as intersection of abnormal and normal minimizing arcs but with distinct minimal times. This phenomenon, already observed in the historical example, is shown to be related to non-continuity property of the value function and is interpreted in a general frame of singularity theory. The third case study concerns a generalization of the navigation of a passive tracer, to generate complex dynamics, in particular with several equators and separatrices. All these cases can be glued to provide a serie of case studies, using combination of mathematical and numerical computations, with nice 2d-representations and using adapted softwares.
This article paves the road to many further studies. First of all, the techniques extend to general Zermelo navigation problems on surfaces of revolution, since the assumption about parallel current can be relaxed, although the Morse-Reeb classification is more complex. Second, the results about the semicubical singularity of the abnormal geodesic when meeting the boundary of the strong domain can be generalized to generic Zermelo navigation problems, since the integrability of the geodesic flow is not an issue [15]. Finally, a new branch of the cut locus is detected in the problem in relation with non-continuity of the value function. A mathematical challenge is to make a complete description of the cut points, at least for 2d-Zermelo navigation problems and a first step towards is a generalization to planar time minimal control problems.
Finally, this article is a step towards the aim of an automatic computations of planar time minimal syntheses, in relation with classification of solutions of dynamical systems. Note in particular the relation between Reeb's graphs [6] and the generalized Morse-Reeb classification in our study.