## 3 Trajectory Optimization Using an Indirect Method
In this context, the time variation in the optimal thrust angle is obtained, through Equation (12), by (numerically) solving the Euler-Lagrange equations [2; 5], which give the \(\tilde{t}\)-derivative of the four costates \(\{\lambda_{\tilde{r}},\,\lambda_{\theta},\,\lambda_{\tilde{u}},\,\lambda_{ \tilde{\phi}}\}\) as  
\[\lambda_{\tilde{r}}^{\prime}=-\frac{\partial\mathcal{H}}{\partial\tilde{r}}= \frac{\bar{v}(\lambda_{\theta}-\lambda_{\tilde{\phi}}\,\tilde{u}+\lambda_{ \tilde{u}}\,\bar{v})}{\bar{r}^{2}}-\frac{2\,\lambda_{\tilde{u}}}{\bar{r}^{3}} \tag{13}\]  
\[\lambda_{\theta}^{\prime}=-\frac{\partial\mathcal{H}}{\partial\theta}=0 \tag{14}\]  
\[\lambda_{\tilde{n}}^{\prime}=-\frac{\partial\mathcal{H}}{\partial\tilde{n}}= \frac{\lambda_{\phi}\,\bar{v}}{\bar{r}}-\lambda_{\tilde{r}} \tag{15}\]  
\[\lambda_{\theta}^{\prime}=-\frac{\partial\mathcal{H}}{\partial\bar{v}}= \frac{\lambda_{\theta}\,\bar{u}-\lambda_{\theta}-2\,\lambda_{\tilde{n}}\,\bar {v}}{\bar{r}} \tag{16}\]  
The expression of the Hamiltonian function is also used to evaluate the transversality condition [2], which in this case gives the following two scalar constraints at the final time  
\[\lambda_{\theta}(\bar{t}_{f})=0\quad,\quad\mathcal{H}(\bar{t}_{f})=1 \tag{17}\]  
In particular, taking Equation (14) into account, we obtain that  
\[\lambda_{\theta}=0 \tag{18}\]during the transfer, so Equations (13) and (16) are simplified as  
\[\mathcal{H}=\lambda_{F}\,\tilde{u}+\lambda_{\theta}\,\frac{\bar{\vartheta}}{ \bar{r}}+\lambda_{\theta}\left(\frac{\bar{\vartheta}^{2}}{\bar{r}}-\frac{1}{ \bar{r}^{2}}\right)-\lambda_{\theta}\,\frac{\tilde{u}\,\bar{\vartheta}}{\bar{r }}+\mathcal{H}_{c} \tag{10}\]0 \[\mathcal{H}=\lambda_{F}\,\tilde{u}+\lambda_{\theta}\,\frac{\bar{\vartheta}}{ \bar{r}}+\lambda_{\theta}\left(\frac{\bar{\vartheta}^{2}}{\bar{r}}-\frac{1}{ \bar{r}^{2}}\right)-\lambda_{\theta}\,\frac{\tilde{u}\,\bar{\vartheta}}{\bar{r }}+\mathcal{H}_{c} \tag{10}\]1  
On the other hand, observing that the dimensionless time \(\bar{t}\) does not compare explicitly in the expression of the Hamiltonian function given by Equation (10), one has that the value of \(\mathcal{H}\) is a constant of motion [2]. Accordingly, the second of Equation (17) can be rewritten as  
\[\mathcal{H}=\lambda_{F}\,\tilde{u}+\lambda_{\theta}\,\frac{\bar{\vartheta}}{ \bar{r}}+\lambda_{\theta}\left(\frac{\bar{\vartheta}^{2}}{\bar{r}}-\frac{1}{ \bar{r}^{2}}\right)-\lambda_{\theta}\,\frac{\tilde{u}\,\bar{\vartheta}}{\bar{r }}+\mathcal{H}_{c} \tag{10}\]2  
Bearing in mind Equations (10) and (11), the optimal control law given by Equation (12), the initial conditions reassumed in Equation (7), and the expression of \(\lambda_{\bar{o}}\) given by Equation (18), the previous equation gives the following constraint at the initial time instant  
\[\mathcal{H}=\lambda_{F}\,\tilde{u}+\lambda_{\theta}\,\frac{\bar{\vartheta}}{ \bar{r}}+\lambda_{\theta}\left(\frac{\bar{\vartheta}^{2}}{\bar{r}}-\frac{1}{ \bar{r}^{2}}\right)-\lambda_{\theta}\,\frac{\tilde{u}\,\bar{\vartheta}}{\bar{r }}+\mathcal{H}_{c} \tag{10}\]3  
where \(\lambda_{\bar{a}}(\bar{t}_{0})\) (or \(\lambda_{\bar{o}}(\bar{t}_{0})\)) is the initial value of the dimensionless costate \(\lambda_{\bar{a}}\) (or \(\lambda_{\bar{o}}\)), which is one of the unknowns of the associated TPBVP. The previous equation states that the point \(\{\lambda_{\bar{a}}(\bar{t}_{0}),\,\lambda_{\bar{o}}(\bar{t}_{0})\}\) lies on a circle of radius equal to \(1/\bar{a}_{m}\), as sketched in Figure 2, where the auxiliary angle \(\delta\in[-\pi/2,\,\pi/2]\) rad is introduced.