Analytical Solution to the Zonal Harmonics Problem Using Koopman Operator Theory
aa r X i v : . [ a s t r o - ph . E P ] D ec Analytical Solution to the Zonal HarmonicsProblem Using Koopman Operator Theory
David Arnas ∗ , Richard Linares † Abstract
This work introduces the use of the Koopman Operator Theoryto generate analytical solutions for the zonal harmonics problem of asatellite orbiting a non spherical celestial body. Particularly, the solu-tion proposed directly provides the osculating evolution of the systemand can be automated to generate any approximation order to thesolution. Moreover, this manuscript defines a modified set of orbitalelements that can be applied to any kind of orbit and that allows theKoopman Operator to generate stable and accurate solutions. In thatregard, several examples of application are included, showing that theproposed methodology can be used in any kind of orbit, including cir-cular, elliptic, parabolic and hyperbolic orbits.
The study of the effects in the motion of satellites subjected to the zonalharmonic components of the gravitational potential of an oblate planet isone of the most elemental problems in celestial mechanics and is of greatimportance both for the design and control of space missions. However,despite being this problem one of the most simple models to study in astro-dynamics (apart from the Keplerian motion), the problem is know to haveno analytical solution [1, 2], and thus, a wide variety of approaches have ap-peared over the years to provide approximate solutions to this problem [3,4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. In this regard, the analytic and semi-analytic solutions are of special interest since they allow us to have a deeperunderstanding of the problem, to obtain a faster long term propagation oforbiting objects, to analyze the stability of orbits [16, 17, 18], or to defineand assess satellite constellations [19, 20, 21] and formation flying [22, 23,24, 25] among other topics.Among the zonal terms of the gravitational potential, the the main satel-lite problem (the study of the motion of a satellite subjected to the J term ∗ Massachusetts Institute of Technology, MA, USA. Email: [email protected] † Massachusetts Institute of Technology, MA, USA. Email: [email protected]
Page 1onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresof the Earth gravitational potential) has always been of special importanceas it is the most important gravitational perturbation for satellites orbitingthe Earth. Examples of that include the solutions proposed by Brouwer [3],Kozai [4], Deprit [6], or Liu [7], which are still extensively used both in the-oretical and applied problems. Particularly, Brouwer [3], based on the vonZeipel perturbation method, proposed a first order solution which was latercomplemented by Lyddane [26] and Cohen [27] to study orbits with eithersmall eccentricity or inclination, and by Coffey [28] to be able to effectivelyassess orbits close to the critical inclination. Kozai [4] proposed, on theother hand, a different approach based on the decomposition of the motionin first-order secular, second-order secular, short-periodic, and long-periodicterms. Afterwards, he extended Brouwer’s formulation by introducing asecond order solution [5] to the main satellite problem. Later, Deprit [6]approached the problem differently by making use of Lie series to definea set of canonical mappings that allows to generate an approximate solu-tion following an iterative set of transformations [29]. This methodologyled to the Lie-Deprit methods, a set of perturbation techniques that hasbeen extensively used in the literature to deal with the main satellite prob-lem [29, 30, 31, 32, 33, 34, 35, 36, 37]. Important examples of application ofthis approach include the method of elimination of the parallax [8], or theelimination of the perigee [11].In this work, we approach the problem from a different perspective. In-stead of performing near unitary transformations using a small parameter,we make use of operator theory to generate solutions to this problem di-rectly. Operator theory is based on functional analysis [38] and focus on thestudy of linear operators in functional spaces [39]. Operator theory is widelyused in theoretical Physics [40] (specially in quantum mechanics) and otherresearch fields [41, 42] such as Fluid-dynamics [43] or control [44, 45, 46], butis not widely used in astrodynamics. Particularly, this manuscript focuseson the use of the Koopman operator, a linear operator introduced by Koop-man [47] and later developed by Newmann [48] that is able to transform anon-linear system in a finite number of dimensions into a linear system inan infinite number of dimensions, on the zonal harmonics problem aroundan oblate celestial body. This allows to linearize the system (benefiting fromthe interesting properties of linear systems) while obtaining a good approx-imation to the solution. The Koopman operator has already been proposedto study problems in celestial mechanics, particularly the motion aroundlibration points in the 3rd body problem [49], and attitude dynamics [50],but it has never been successfully applied to the zonal harmonics problemaround a celestial body.To that end, this work presents the methodology to apply the Koopmanoperator to the zonal harmonics problem. To be more precise, the contribu-tions of this work are as follows; 1) a new set of orbital elements is developedbased on Arnas [51], 2) the development of a closed-form computation of thePage 2onal Harmonics Motion Using Koopman Operator D. Arnas, R. LinaresKoopman matrix using Legendre polynomials, 3) numerical error analysis ofthe analytical Koopman-based solution applied to the motion of a satelliteorbiting a non-spherical celestial body. Particularly, this manuscript makesuse of a modified set of orbital elements from Arnas [51] especially devisedfor their application alongside the Koopman operator. This set of orbitalelements is able to represent any kind of orbit, including circular, elliptic,parabolic and hyperbolic orbits. In that regard, several examples of applica-tion are included to show the performance of the methodology for differentorbits. Additionally, this work proposes the use of Legendre polynomials tocompute the Koopman matrix of the system. This provides several advan-tages, including lower computation complexity, easier implementation, anda better distribution of the error in the variable range defined.This manuscript is organized as follows. First, we include a summaryof the Koopman operator theory and the methodology followed to applyit to the zonal harmonics problem. Second, we propose two sets of orbitalelements to represent the dynamical system. The particularity of thesesets of orbital elements is that they allow us to represent the differentialequation governing the system in completely polynomial form, being, inaddition, their expressions cuasi-linear for the application of the Koopmanoperator methodology. These set of orbital elements are a set of variablesespecially devised for their use with the Koopman operator due to theirparticular properties. Third, we apply this methodology to a collectionof four common orbits used in astrodynamics, namely, a sun-synchronousfrozen orbit, a molniya orbit, a hyperbolic orbit, and a near equatorial orbit.This allows us to show the performance of the methodology presented in thiswork and its possibilities.
We depart from the classical initial value problem, represented by an au-tonomous system of differential equations expressed as: ( ddt x ( t ) = f ( x ) x ( t ) = x ; (1)where x ∈ R d is the set of d variables that evolves with time t and that definethe state of the system. In here, f : R d → R d represents the mathematicalmodel of a given dynamical system, d is the number of dimensions of theproblem, and x are the initial conditions.Let g ( x ) be an observable of the space, or in other words, g ( x ) is afunction in the set of variables x . Then, if we study the evolution of g ( x ) inthe dynamic defined by Eq. (1) we obtain: ddt g ( x ) = K ( g ( x )) = ( ∇ x g ( x )) ddt x ( t ) = ( ∇ x g ( x )) f ( x ) , (2)Page 3onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linareswhere ∇ x g = [ ∂g/∂x , ∂g/∂x , . . . , ∂g/∂x d ], and K ( g ( x )) is defined as theKoopman operator applied to an observable. This means that the evolutionof any observable under the effects of the dynamical system is provided bythe Koopman operator: K ( · ) = ( ∇ x · ) f ( x ) . (3)where ( ∇ x · ) = [ ∂ ( · ) /∂x , ∂ ( · ) /∂x , . . . , ∂ ( · ) /∂x d ]. At this point, it is im-portant to note that since the Koopman operator was generated applyingthe chain rule in the derivative of the observable, the operator is in factlinear [47], and thus: K ( C g ( x ) + C g ( x )) = C K ( g ( x )) + C K ( g ( x )) , (4)where C and C are two constants.The main idea of the Koopman operator methodology is to find an ap-proximated linear representation of the full dynamical system, and then usethe property of linearity to obtain the exact solution to this approximation.To that end, the Koopman operator makes use of an extended configurationspace of dimension m > d that is able to approximate the non-linearitiesof the original system. In fact, if we could be able to work with the com-plete infinite Hilbert space [44], that is, a configuration space in an infinitenumber of dimensions, the nonlinearities of the system would be completelyrepresented by this linear operator. However, since we are constrained bythe use of a limited set of dimensions in the extended space, the lineariza-tion only provides an approximation of the original system that is better thecloser is the system in the extended space to a linear one.In order to perform this linearization, this manuscript makes use of theGarlekin method to obtain an approximated linearized system of equations.The Garlekin method is based on the idea of performing this approxima-tion by finding the projection of the original system into a set of orthogonalfunctions that define the extended space of configuration for the Koopmanoperator. This projection is performed using inner products between func-tions. Particularly, let f ( x ) be the function to represent in the dynamic,and let g ( x ) be a basis function. Then, the inner product between bothfunctions is represented by: h f, g i = Z Ξ f ( x ) g ( x ) w ( x ) d x , (5)where w ( x ) is a weighting associated with the basis functions, and Ξ is thedomain in which the basis functions are defined. That way, the error resul-tant in the approximation is orthogonal to the configuration space definedby the set of basis functions. Therefore, a set of orthogonal basis functionshas to be defined in the extended space of dimension m . To that end, thiswork makes use of Legendre polynomials as basis functions. The reason forthat is the multiple advantages that they provide when computing the innerPage 4onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresproducts. First, they are polynomials, which simplifies the computation.Second, Legendre polynomials are defined in a bounded domain ([ − , R . Third, their weighting function is a constant, a fact that simplifiesthe expression of the inner product. This allows to ease the automatizedcomputation of all the inner products required in the methodology.Let L ( x ) be the set of all the basis functions of the problem in theoriginal variables written in vector form L ( x ) = [ L ( x ) , ..., L m ( x )] T , where L i ( x ) represents the basis function i from the set. Then, the derivative of abasis function L i ( x ) can be obtained using the Koopman operator: K ( L i ) = ddt L i = ( ∇ x L i ) f , (6)where it is important to note that ∇ x L i ( x ) can be represented by a combi-nation of basis functions from L ( x ) due to the orthogonality and polynomialnature of the basis functions. Therefore, our objective now is to approximatethis derivative with the basis functions defined in L ( x ). This is done usingthe inner product operation. Let K be the Koopman matrix, a matrix ofsize m × m whose components K ij contain the value of the projection of thederivative of the basis function i into the basis function j . In other words: K ij = h ( ∇ x L i ) f , L j i = Z Ξ ( ∇ x L i ) f L j wd x . (7)For this work, these inner products are always performed between polyno-mials, and thus, it is possible to directly compute the integral analyticallywhich significantly improves the computational speed of the methodologypresented. Then, by following this procedure for all the combinations ofbasis functions, a linear system in the form: ddt L = K L . (8)is generated, which represents the approximated linear problem of the dy-namic. Note that by using higher orders of basis functions, the Koopmanmatrix increases in size, but it is better able to better represent the non-linearities of the original system.Once this linearization is performed, it is only required to solve the linearsystem represented by Eq. (8). This can be done, for instance, using thespectral theorem when the matrix K is diagonalizable. In particular, let E i and v i be the eigenvalues and corresponding left eigenvectors of the matrix K . In addition, let V be the matrix containing all the left eigenvectors inits rows. Then, by pre-multiplying V in Eq. (8) we obtain: V ddt L = EK L = EV L , (9)Page 5onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linareswhere E is a diagonal matrix containing all the eigenvalues of the system.We know that V is a constant matrix since it is composed by the eigenvectorsof the system, therefore: V ddt L = ddt ( V L ) = EV L , (10)where we can perform the algebraic variable transformation Φ = V L with φ i = v i L , called the Koopman Canonical Transform [45], to obtain: ddt Φ = E Φ , (11)whose solution is: Φ ( t ) = exp( Et ) Φ ( t ) . (12)In this regard, it is important to note that the functions Φ = V L are in factthe eigenfunctions of the space defined by the dynamical system or, in otherwords, they are the set of functions defining the orthogonal directions of theHilbert space of dimension m subjected to the dynamic of the system.The objective is now to recuperate the original variables from the system.From Eq. (12) we can obtain the evolution of the basis functions L : L ( t ) = V − Φ = V − exp( Et ) V L ( t ) . (13)In addition, the original variables x can be represented using the basis func-tions L though a projection on this basis. Let T be a transformation matrixof size d × m defined as: T ij = h x i , L j i = Z Ξ x i L j wd x , (14)that is, the component T ij represents the projection of the original variable x i with i ∈ { , . . . , d } into the basis function j with j ∈ { , . . . , m } . T isalso called the Koopman modes matrix since it relates the original variablesfrom the problem with the eigendecomposition of the dynamic, and thus,the frequencies of the motion. Therefore, the original set of variables x ( t )can be expressed as: x ( t ) = T L ( t ) = T V − exp( Et ) V L ( x ) , (15)where x are the initial conditions of the original set of variables. In order to efficiently apply the Koopman operator to a perturbed system,we require that the system behaves like a linear system with a small pertur-bation. This is due to the fact that for a real application of the KoopmanPage 6onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresoperator we can only use a limited number of basis functions, and thus, thecloser we are to the linear system, the faster the perturbation method willconverge to the solution of the problem. Therefore, we require a formula-tion that is able to fully represent the zonal harmonics perturbation in thesystem and that is close to linear for the solution to converge quickly. Inthat regard, Ref [51] introduced a new set of orbital elements that fulfillsthe conditions we are looking for. However, due to the particularities ofthe Koopman operator methodology, some modifications to that formula-tion have to be performed. For that reason, we include a summary of theformulation presented in Ref [51] as a reference for this manuscript.The Hamiltonian in spherical coordinates of a satellite orbiting the Earthand subjected to the zonal harmonics terms of the gravitational perturbationis: H = 12 p r + p ϕ r + p λ r cos ( ϕ ) ! − µr − m X n =2 J n P n (sin( ϕ )) R n ⊕ r n ! , (16)where µ is the gravitational constant of the Earth, R ⊕ is the Earth Equa-torial radius, r is the radial distance from the satellite to the center of theEarth, ϕ is the latitude, λ is the longitude of the satellite in the inertialframe of reference for a given instant, J n are the zonal terms of the Earthgravitational potential, and P n (sin( ϕ )) are the Legendre polynomials of or-der n applied to the variable sin( ϕ ). In addition, the conjugate momenta ofthese variables are: p r = ˙ r ; p ϕ = r ˙ ϕ ; p λ = r cos ( ϕ ) ˙ λ ; (17)where we defined ˙ x as the derivative with respect to time of variable x . Theassociated Hamilton equations of this hamiltonian are: drdt = p r ; dp r dt = − µr + p ϕ r + p λ r cos ( ϕ ) + m X n =2 ( n + 1) µJ n P n (sin( ϕ )) R n ⊕ r n +2 ; dϕdt = p ϕ r ; dp ϕ dt = − p λ r sin( ϕ )cos ( ϕ ) − m X n =2 µJ n ∂P n (sin( ϕ )) ∂ϕ R n ⊕ r n +1 ; dλdt = p λ r cos ( ϕ ) ; dp λ dt = 0 . (18)Page 7onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linareswhich as it can be seen, are highly non-linear even for the unperturbed prob-lem and thus, they are not appropriate for their direct use of the Koopmanoperator. To solve this, Ref. [51] proposed to perform a time regularizationdefined by: dθdt = p θ r . (19)where, p θ = s p ϕ + p λ cos ( ϕ ) , (20)is the angular momentum of the orbit, followed by a set of variable trans-formations defined by: α = p θ r − µp θ = r q ˙ ϕ + cos ( ϕ ) ˙ λ − µr q ˙ ϕ + cos ( ϕ ) ˙ λ ; p r = ˙ r ; s = sin( ϕ ); γ = p ϕ p θ cos( ϕ ) = ˙ ϕ cos( ϕ ) q ˙ ϕ + cos ( ϕ ) ˙ λ ; I θ = s p ϕ + p λ cos ( ϕ ) = 1 r q ˙ ϕ + cos ( ϕ ) ˙ λ ; β = λ − arcsin tan( ϕ ) s p λ p θ − p λ ! = λ − arcsin tan( ϕ ) ˙ λ s cos ( ϕ )˙ ϕ + cos ( ϕ ) ˙ λ − cos ( ϕ ) ˙ λ ! ; ξ = p λ p θ − p λ = cos ( ϕ ) ˙ λ ˙ ϕ + cos ( ϕ ) ˙ λ − cos ( ϕ ) ˙ λ ; p λ = r cos ( ϕ ) ˙ λ. (21)Using this expanded set of orbital elements, the system of differential equa-Page 8onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linarestions from Eq. (27) transforms into [51]: dαdθ = − p r − m X n =2 µJ n R n ⊕ ∂P n ( s ) ∂s γI n +1 θ ( α + µI θ ) n − ( α + 2 µI θ ) ; dp r dθ = α + m X n =2 ( n + 1) µJ n R n ⊕ P n ( s ) I n +1 θ ( α + µI θ ) n ; dsdθ = γ ; dγdθ = − s − m X n =2 µJ n R n ⊕ ∂P n ( s ) ∂s p λ I n +3 θ ( α + µI θ ) n − ; dI θ dθ = m X n =2 µJ n R n ⊕ ∂P n ( s ) ∂s γI n +2 θ ( α + µI θ ) n − ; dβdθ = − m X n =2 µJ n R n ⊕ ∂P n ( s ) ∂s p λ sξI nθ ( α + µI θ ) n − ; dξdθ = m X n =2 µJ n R n ⊕ ∂P n ( s ) ∂s p λ γξ I n − θ ( α + µI θ ) n − ; dp λ dθ = 0 . (22)which, as it can be seen, it is completely polynomial for the variables con-sidered, since p λ is a constant of the motion. This property will allow usto perform the direct analytical integration of the inner products definedbefore following a simple and fast process.Unfortunately, this set of equations cannot be used directly with theKoopman operator for two reasons. First, we are looking for a Koopmanoperator matrix that is independent on the initial conditions since this willallow us to apply the Koopman matrix obtained by the Garlekin methodto any orbit. This means that we have to include p λ as a variable in theequations, making the system from Eq. (22) not linear any more under thisconsideration. Second, when using the Garlekin method we are performing aset of inner products of functions in a given range. Therefore, we require thatthe variables of the system are normalized in order to bound the maximumvalues as much as possible. In that regard, the variable ξ from the previousformulation has to be transformed since it is not bounded depending on theinitial conditions of the orbit. In this section we introduce a modified set of elements that allows to applythe Koopman operator while obtaining a good performance in the solution.Page 9onal Harmonics Motion Using Koopman Operator D. Arnas, R. LinaresTo that end, the final goal is to find a Koopman matrix that allows torepresent the widest range of orbits possible. In that regard, and since thevariable β is not well defined for equatorial orbits, the space of solutions willbe separated in two regions. The first one covers the general formulation fororbits that are far from the equatorial orbit. In particular, this region coversany orbit whose inclination is larger than 15 degrees. On the other hand,the second region covers the close to equatorial orbits up to inclinations of20 degrees. That way, no matter the initial conditions of the satellite, itsorbit will be able to be studied, having, in addition, an overlap between bothregions in case a mission is flying near the boundary between the definedregions. Moreover it is important to note that these regions can be extendedeven further at a cost of a gradual lost in precision in the solution. In order to overcome the limitations of the formulation provided in Eq. (22)when using the Koopman operator, we propose the following modified setof variables for the application in the zonal harmonics problem: Λ = s R ⊕ µ (cid:18) p θ r − µp θ (cid:19) = α s R ⊕ µ ; η = p r s R ⊕ µ ; s = sin( ϕ ); γ = p ϕ p θ cos( ϕ ); κ = s µR ⊕ p ϕ + µR ⊕ p λ cos ( ϕ ) = I θ p µR ⊕ ; β = λ − arcsin tan( ϕ ) s p λ p θ − p λ ! ; χ = p λ p θ (cid:16) µR ⊕ p ϕ + µR ⊕ p λ cos ( ϕ ) (cid:17) / sin ( ϕ ) + p λ p θ cos ( ϕ ) = p λ I θ ( µR ⊕ ) / s + γ ; ρ = p λ p θ = p λ I θ . (23)Note that in this set of equations, we have normalized several of the originalorbital elements making them non-dimensional, and substituted the variable ξ for the element χ which behaves much better than the original variabledue to the fact that its value is more constrained for a wide variety of orbits.It is also important to remark that χ is a first integral of the unperturbedproblem as in the case of the elements κ , β , and ρ . Page 10onal Harmonics Motion Using Koopman Operator D. Arnas, R. LinaresWhen dealing with the general zonal harmonics problem around anoblate celestial body, the derivatives of this new set of orbital elements are: dΛdθ = − η − m X n =2 J n ∂P n ( s ) ∂s γκ n +1 ( Λ + κ ) n − ( Λ + 2 κ ) ; dηdθ = Λ + m X n =2 ( n + 1) J n P n ( s ) κ n +1 ( Λ + κ ) n ; dsdθ = γ ; dγdθ = − s − m X n =2 J n ∂P n ( s ) ∂s ρ κ n +1 ( Λ + κ ) n − ; dκdθ = m X n =2 J n ∂P n ( s ) ∂s γκ n +2 ( Λ + κ ) n − ; dβdθ = − m X n =2 J n ∂P n ( s ) ∂s sχκ n − ( Λ + κ ) n − ; dχdθ = m X n =2 J n ∂P n ( s ) ∂s γχκ n +1 ( Λ + κ ) n − + m X n =2 J n ∂P n ( s ) ∂s ρχ κ n − ( Λ + κ ) n − = m X n =2 J n ∂P n ( s ) ∂s (cid:0) γκ + ρχ (cid:1) χκ n − ( Λ + κ ) n − ; dρdθ = m X n =2 J n ∂P n ( s ) ∂s γρκ n +1 ( Λ + κ ) n − . (24)As it can be seen the former system of differential equations is completelypolynomial (note also that the derivatives of the Legendre polynomials P n are also polynomials in s ), and, in addition, all the coefficients of the polyno-mial are either 1 (corresponding to the unperturbed problem), or a multipleof J n by a rational number (corresponding to each one of the terms fromthe zonal harmonics). In general, Eq. (23) can be used to transform spherical coordinates to theset of orbital elements proposed. However, the element β requires a specialattention in its definition due to the non uniqueness of its definition due thePage 11onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresarcsine function. This problem can be solved easily by defining β as: β = λ − arcsin tan( ϕ ) s p λ p θ − p λ ! if p ϕ ≥ λ + arcsin tan( ϕ ) s p λ p θ − p λ ! + π if p ϕ < . (25)which in that case has the same value of the right ascension of the ascendingnode of the orbit.For the inverse transformation, i.e., to transform these set of orbitalelements to spherical coordinates, we use the orbital element definitionsprovided by Eqs. (21) and (23), and the geometrical relations that existbetween these variables. Particularly, the spherical coordinates can be ob-tained following this sequence of operations: p θ = √ µR ⊕ κ ; r = p θ Λ r µR ⊕ + µp θ ; p r = r µR ⊕ η ; ϕ = arcsin( s ); p ϕ = γ cos( ϕ ) p θ ; si = p s + γ ; u = arcsin (cid:16) ssi (cid:17) if p ϕ ≥ π − arcsin (cid:16) ssi (cid:17) if p ϕ < λ = β ± arccos (cid:18) cos( u )cos( ϕ ) (cid:19) if s ≥ β ∓ arccos (cid:18) cos( u )cos( ϕ ) (cid:19) if s < p λ = ρp θ ; (26)where si is the sine of the osculating inclination of the orbit, u is the argu-ment of latitude of the orbit at a given time, and ± and ∓ represent thedifferent solutions when considering prograde (upper symbol) or retrograde(lower symbol) orbits. Page 12onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares In the cases where only the J H = 12 p r + p ϕ r + p λ r cos ( ϕ ) ! − µr + 12 µR ⊕ J r (cid:0) ( ϕ ) − (cid:1) , (27)whose associated Hamilton equations are: drdt = p r ; dp r dt = − µr + p ϕ r + p λ r cos ( ϕ ) + 32 µJ R ⊕ r (cid:0) ( ϕ ) − (cid:1) ; dϕdt = p ϕ r ; dp ϕ dt = − p λ r sin( ϕ )cos ( ϕ ) − µJ R ⊕ r sin( ϕ ) cos( ϕ ); dλdt = p λ r cos ( ϕ ) ; dp λ dt = 0 . (28)Then, if we perform the time regularization, the expansion of the configura-tion space and the variable transformation described for the zonal problemwe obtain the following system of equations: dΛdθ = − η − J sγκ ( Λ + κ ) ( Λ + 2 κ ) ; dηdθ = Λ + 32 J κ ( Λ + κ ) (3 s − dsdθ = γ ; dγdθ = − s − J sρ κ ( Λ + κ ) ; dκdθ = 3 J sγκ ( Λ + κ ) ; dβdθ = − J s χ ( Λ + κ ) ; dχdθ = 12 J sγχκ ( Λ + κ ) + 6 J sρχ ( Λ + κ ) ; dρdθ = 3 J sγρκ ( Λ + κ ) . (29)which is again completely polynomial. Page 13onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares Orbits that are close to the Equator present two difficulties when using theformer formulation. First, the variable β is not well defined in a similar wayin which the right ascension of the ascending node in keplerian variablesis not well defined for equatorial orbits. Second, the value of the variable χ can increase too much for the Garlekin method to properly approximatethe variable in the ranges of application. For this reason, we use a differentformulation for orbits that are close to equatorial.We start with the Hamilton equations expressed in spherical coordinatesprovided by Eq. (18). From this system of differential equations we perform atime regularization, but instead of using θ as the time regularization variableas in the previous case, we use a new independent variable τ defined by: dτdt = p θ r cos ( ϕ ) , (30)and introduce it in the system of differential equations from Eq. (18) toobtain: drdτ = p r p θ r cos ( ϕ ); dp r dτ = − µp θ cos ( ϕ ) + 1 r p ϕ p θ cos ( ϕ ) + 1 r p λ p θ + m X n =2 ( n + 1) µJ n P n (sin( ϕ )) R n ⊕ r n p θ cos ( ϕ ); dϕdτ = p ϕ p θ cos ( ϕ ); dp ϕ dτ = − p λ p θ tan( ϕ ) − m X n =2 µJ n ∂P n (sin( ϕ )) ∂ϕ p θ R n ⊕ r n − cos ( ϕ ); dλdτ = p λ p θ ; dp λ dτ = 0 . (31)Similarly to the previous section, we perform the variable transformationfor Λ , η , s , γ , and ρ , and expand the configuration space of the problem withvariable κ . In this regard, note that variable χ is not required to representthis region of the space. Following these transformations, the system ofPage 14onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresdifferential equations becomes: dΛdτ = − η (cid:0) − s (cid:1) − m X n =2 J n ∂P n ( s ) ∂s γκ n +1 ( Λ + κ ) n − ( Λ + 2 κ ) (cid:0) − s (cid:1) ; dηdτ = Λ (cid:0) − s (cid:1) + m X n =2 ( n + 1) J n P n ( s ) κ n +1 ( Λ + κ ) n (cid:0) − s (cid:1) ; dsdτ = γ (cid:0) − s (cid:1) ; dγdτ = − s (cid:0) − s (cid:1) − m X n =2 J n ∂P n ( s ) ∂s ρ κ n +1 ( Λ + κ ) n − (cid:0) − s (cid:1) ; dκdτ = m X n =2 J n ∂P n ( s ) ∂s γκ n +2 ( Λ + κ ) n − (cid:0) − s (cid:1) ; dλdτ = ρ ; dρdτ = m X n =2 J n ∂P n ( s ) ∂s γρκ n +1 ( Λ + κ ) n − (cid:0) − s (cid:1) . (32)Note that this system of differential equations is still polynomial, but con-trary to what happens in the previous formulation, the non-perturbed prob-lem is non-linear. This may seem like a limiting factor for the applicationof the Koopman operator, however, the non-linearities of the unperturbedproblem are in fact small when compared with the linear part, making thesystem nearly linear. This property allows us to directly apply the Koopmanoperator to this system of differential equations while still obtaining a goodperformance in the results. Unfortunately, even with these transformations,the solution for the variables s and γ provided by the Koopman operatorperform poorly as we increase the inclination of the orbit. Therefore, theseorbital elements must be modified to improve the performance. This is donevia a normalization of these variables.Let ψ be a constant value such that ψ > sin( i ) where i is the maximuminclination of the orbits we are interested in. In our case, and since wewant to cover all the space of solutions including some overlap, we define ψ = sin(20 ◦ ). Note that this value can be reduced to improve performanceif we are interested in orbits with lower inclinations. In particular, the closerthe value of ψ is to sin( i ) of the orbit, the better is the performance obtainedas long as ψ > sin( i ). On the other hand, it is not convenient to increasethe value of ψ over 1 since that will make the differential equations morenon-linear, reducing in consequence the performance of the solution.Page 15onal Harmonics Motion Using Koopman Operator D. Arnas, R. LinaresBased on this constant ψ we define two transformations: σ = sψ ;Γ = γψ ; (33)and introduce them into the set of differential equations to obtain: dΛdτ = − η (cid:0) − ψ σ (cid:1) − m X n =2 ψJ n ∂P n ( ψσ ) ∂ ( ψσ ) Γ κ n +1 ( Λ + κ ) n − ( Λ + 2 κ ) (cid:0) − ψ σ (cid:1) ; dηdτ = Λ (cid:0) − ψ σ (cid:1) + m X n =2 ( n + 1) J n P n ( ψσ ) κ n +1 ( Λ + κ ) n (cid:0) − ψ σ (cid:1) ; dσdτ = Γ (cid:0) − ψ σ (cid:1) ; d Γ dτ = − σ (cid:0) − ψ σ (cid:1) − m X n =2 ψ J n ∂P n ( ψσ ) ∂ ( ψσ ) ρ κ n +1 ( Λ + κ ) n − (cid:0) − ψ σ (cid:1) ; dκdτ = m X n =2 ψJ n ∂P n ( ψσ ) ∂ ( ψσ ) Γ κ n +2 ( Λ + κ ) n − (cid:0) − ψ σ (cid:1) ; dλdτ = ρ ; dρdτ = m X n =2 ψJ n ∂P n ( ψσ ) ∂ ( ψσ ) Γ ρκ n +1 ( Λ + κ ) n − (cid:0) − ψ σ (cid:1) ; (34)which can finally be applied directly to the Koopman operator formulationfor any orbit whose inclination is smaller than arcsin( ψ ). Note also that nowthe coefficients of the polynomials of the unperturbed problem that are non-linear ( ψ ) are small compared to the linear term. In fact, the Koopmanoperator treats these terms as another perturbation of the linear system,which is why it is possible to apply this non-linear formulation. As in the zonal case, we include the particular case of the J term of thegravitational potential for completeness and due to the fact that it is themost important term for orbiting satellites around the Earth. In that re-gard the system of differential equations when taking into account just thePage 16onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresoblateness of the Earth ( J term of the Earth gravitational potential) is: dΛdτ = − η (cid:0) − ψ σ (cid:1) − ψ J σ Γ κ ( Λ + κ ) ( Λ + 2 κ ) (cid:0) − ψ σ (cid:1) ; dηdτ = Λ (cid:0) − ψ σ (cid:1) + 32 J κ ( Λ + κ ) (cid:0) ψ σ − (cid:1) (cid:0) − ψ σ (cid:1) ; dσdτ = Γ (cid:0) − ψ σ (cid:1) ; d Γ dτ = − σ (cid:0) − ψ σ (cid:1) − J σρ κ ( Λ + κ ) (cid:0) − ψ σ (cid:1) ; dκdτ = 3 ψ J σ Γ κ ( Λ + κ ) (cid:0) − ψ σ (cid:1) ; dλdτ = ρ ; dρdτ = 3 ψ J n σ Γ ρκ ( Λ + κ ) (cid:0) − ψ σ (cid:1) . (35) In this section we show the result of applying the methodology described inthis work to different orbits. In particular, we show three cases of applica-tion: a near circular sun-synchronous orbit, a Molniya orbit, a hyperbolicorbit, and a near-equatorial orbit. This will allow us to show the perfor-mance of the formulation for some characteristic orbits. In addition, wealso include the results for different orders of basis functions to show howthe Koopman operator methodology performs and converges to the solutionof the system of differential equations. In that regard, it is important tonote that the set of differential equations provided by Eqs. (29) and (35)are based on polynomials in odd powers. This means that only the basisfunctions containing terms in an odd power in the polynomial contributeto the representation of the differential equation. For that reason, in theexamples provided we focus on the use of basis function with a maximumorder that is odd.The examples provided focus on LEO orbits due to the fact that thezonal harmonics perturbation is larger in this region. Therefore, a much bet-ter performance should be expected as the altitude of the orbits increases.Moreover, all of the examples included are compared to a numerical inte-gration of the system using a Runge-Kutta scheme of order 4-5 to show theperformance of the methodology for different scenarios.
For this example, a sun-synchronous frozen Earth observation orbit is se-lected with the initial orbital conditions: semi-major axis a = 7077 .
722 km,eccentricity e = 0 . i = 98 .
186 deg, argument of perigeePage 17onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares ω = 90 . . ν = 0 . r ), the latitude ( ϕ ), and the longitude ( λ ) of the orbit as a functionof the independent variable θ (defined as the time regularization variable inEq. (19)). Moreover, Figure 1 (right) shows the comparison of the Koop-man operator solution using basis functions of order 7 with the numericalintegration of the equations. From these results, it can be observed thatthe error of this solution is in the order of magnitude of the meters in theposition of the satellite (both taking into account both the radial distanceand the angles latitude and longitude of the orbit). -3 -4 -3 Figure 1: 7th order solution and error in { r, ϕ, λ } for a frozen sun-synchronous orbit.The result from Figure 1 can be improved by increasing the order of thebasis functions used in the Koopman operator approach. In that regard,Figure 2 shows the errors of the solution for the 9th (left) and 11th (right)order basis function solutions when compared to the numerical integration.As it can be seen, the maximum error is reduced to a couple of meters andless than 0.32 meters respectively. This shows that the solution is gettingcloser to the real evolution of the system by increasing the order of the basisfunctions used to represent the system of differential equations.Moreover, we can also study the long term evolution of these solutions.To that end, Figure 3 shows the error for a propagation of 100 orbital revo-lutions when using a set of basis functions of order 7th (left) and 9th (right).As it can be observed, the error of the solution quickly increases with time.Nevertheless, the error reduces with the increase on the order of the basisfunctions. Thus, it is possible to improve the long term performance of themethodology by further increasing the order of the basis functions used.Additionally, it is also possible to study the spectral structure of thePage 18onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares -3 -5 -4 -4 -6 -5 Figure 2: 9th and 11th order solution error in { r, ϕ, λ } for a frozen sun-synchronous orbit. Number of orbital periods -505 0 20 40 60 80 100-2-100 20 40 60 80 100-0.500.50 20 40 60 80 100
Number of orbital periods -0.500.511.5
Figure 3: 7th and 9th order solution error for a long term propagation in { r, ϕ, λ } for a frozen sun-synchronous orbit.solution. To that end, Figure 4 shows the distribution of eigenvalues ofthe system when using a set of basis functions of order 7 (left) and 9 (right)respectively. In that regard, it is important to note that, since the Koopmanmatrix is independent of the initial conditions of the orbit, these eigenvaluesare valid for any orbit that the dynamical system provided by Eq. (29) isable to represent. As it can be seen in the figure, the maximum values ofthe frequencies are 7 and 9 respectively. In that respect, we know that theeigenvalues of the unperturbed problem using the system from Eq. (29) are ± ı and 0. However, since the basis functions are built as a combination ofPage 19onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresthe initial set of variables, the maximum frequency of the system increaseswith the order of the basis functions used. Particularly, let φ and φ betwo eigenfunctions of the Koopman matrix with associated eigenvalues E and E respectively, therefore: dφ dt = E φ ; dφ dt = E φ . (36)Then, we can define a third function g as the multiplication of the first two,that is: g = φ φ . If we introduce it in the dynamical system we obtain: dgdt = dφ dt φ + φ dφ dt = E φ φ + E φ φ = ( E + E ) φ φ = ( E + E ) g, (37)which means that g is also an eigenfunction of the system with associatedeigenvalue E + E . This is the reason why the maximum frequency that weobserve in Figure 4 matches the maximum order of the basis functions used.On the other hand, Figure 4 also shows that there are some eigenvalues withreal part different to zero. These eigenvalues are generated by the multiplecombinations of the basis functions related with the secular variation of theright ascension of the ascending node (represented by variable β ) with therest of variables involved in the problem. -0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 Real part of the eigenvalues -8-6-4-202468 I m ag i na r y pa r t o f t he e i gen v a l ue s -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 Real part of the eigenvalues -10-50510 I m ag i na r y pa r t o f t he e i gen v a l ue s Figure 4: 7th and 9th order eigenvalues for the general formulation case.
We select a typical molniya orbit as an example of application to eccentricorbits. Particularly, we set the initial conditions for the propagation to: a = 26600 . e = 0 . i = 63 .
435 deg, ω = 270 . . ν = 0 . { r, ϕ, λ } andtheir error when using a set of basis functions of order 7. As it can be seen,the maximum error observed (less than 400 m) is now bigger than in thecase of the sun-synchronous orbit. Nevertheless, this is consistent with theresults from the previous example having a relative error in the order ofmagnitude of 10 − . Page 20onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares -4 -4 Figure 5: 7th order solution and error in { r, ϕ, λ } for a molniya orbit. -5 -5 -3 -6 -4 Figure 6: 9th and 11th order solution error in { r, ϕ, λ } for a molniya orbit.In addition, Figure 6 shows the error on the position when using a set ofbasis functions of order 9 (left) and 11 (right) respectively. As it can be seen,the error in this case is of 10s of meters. This shows that the formulationand methodology presented can be applied to very eccentric orbits whileobtaining a good accuracy. For the example of hyperbolic orbit, we select the following initial conditions: a = − . e = 1 . i = 50 . ω = 0 . . ν = 0 . θ = 120 ◦ in order to prevent thevariables to become extremely large. Nevertheless, it can be seen from thefigures that the precision of the solution is also maintained in this case,having a maximum error on the order of magnitude of meters even for thehyperbolic case. -3 -5 -5 Figure 7: 7th order solution and error in { r, ϕ, λ } for a hyperbolic orbit. -6 -4 -3 -7 -5 Figure 8: 9th and 11th order solution error in { r, ϕ, λ } for a hyperbolic orbit.Note that sometimes applying a set of higher order basis functions pro-duces a slightly decrease in accuracy in the methodology, especially in thelongitude of the orbit. This can be seen, for instance, in the errors in lati-Page 22onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linarestude of Figure 2 for the cases of 9th and 11th orders, or in error in latitudefor this hyperbolic example for orders 7 and 9 (Figures 7 and 8. However, ifthe order of the basis functions is further increased, this anomalous behaviordisappears. An example of that can be seen in Figure 8 in the latitude ofthe solution for basis functions of orders 9 and 11. This happens due tonumerical errors appearing during the eigendecomposition of the Koopmanmatrix and subsequent inversion of the matrix of eigenvectors. This final example is included to show the performance of the formulationfor close to equatorial orbits. Particularly, the initial conditions selectedfor this examples are: a = 7192 .
15 km, e = 0, i = 5 deg, ω = −
90 deg,Ω = 0 . ν = 180 . -4 -5 Figure 9: 7th order solution and error in { r, ϕ, λ } for a 5 deg inclinationorbit.Additionally, and since the Koopman matrix is different from the oneused in the previous examples (due to the fact that we are using a slightlydifferent formulation), it is interesting to study the difference in the distri-bution of eigenvalues of this system. In that regard, Figure 11 shows a mapcontaining all the eigenvalues when using eigenfunctions of order 7 and 9respectively. As it can be seen, the distribution is different from Figure 4,particularly the real parts of the eigenvalues. This is due to the differentPage 23onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares -3 -4 -4 -4 -5 -5 Figure 10: 9th and 11th order solution error in { r, ϕ, λ } for a 5 deg inclina-tion orbit.nature of both formulations. In the general case, we have as an orbitalelement a variable whose secular value matches the one of the right ascen-sion of the ascending node. However, in the close to equatorial case, thatvariable was substituted by another whose derivative was linear with theratio p λ /p θ . Additionally, the close to equatorial formulation does not re-quire the additional variable χ which also generates a large number of theseeigenvalues with non-zero real part. This shows the importance of selectingan adequate formulation for the use alongside with the Koopman operatorsince it greatly depends on the eigendecomposition of the system. -8 -6 -4 -2 0 2 4 6 8 Real part of the eigenvalues -3 -8-6-4-202468 I m ag i na r y pa r t o f t he e i gen v a l ue s -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 Real part of the eigenvalues -10-50510 I m ag i na r y pa r t o f t he e i gen v a l ue s Figure 11: 7th and 9th order eigenvalues for the close to equatorial formu-lation.
This manuscript focuses on the application of the Koopman operator the-ory to the zonal harmonics problem of a satellite orbiting an oblate celestialPage 24onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linaresbody. This represents the first attempt of using operator theory in this im-portant problem in astrodynamics. To that end, this work presents a sum-mary of the Koopman operator theory and how to apply it effectively intothe zonal harmonics problem using a set of modified orbital elements thatallow to represent any kind of orbit, including, circular, elliptic, parabolicand hyperbolic orbits.The Koopman operator methodology allows us to obtain in an auto-mated process any order of the solution by an extension of the set of basisfunctions used to represent the solution. Therefore, this methodology allowsthe error to be adjusted to any precision required. In that regard, four exam-ples of application are shown in this manuscript, showing the performanceresults of this methodology to very characteristic orbits in celestial mechan-ics, including a sun-synchronous frozen orbit, a molniya orbit, a hyperbolicorbit and a close to equatorial orbit. These examples show that it is possibleto obtain relative errors of 10 − when using this methodology for the differ-ent cases of application. Similar performance have been observed for otherorbits in the same range of altitudes. In that regard, it is important to notethat as the altitude of the orbit increases, the perturbation produced bythe zonal harmonic terms of the gravitational potential is also reduced, andthus, the performance of the methodology also improves accordingly. Forinstance, for a geo-synchronous orbit, less than 1 cm errors were observedfor the same order of basis functions considered in this work.Compared to other methodologies that focus on the zonal harmonicsproblem, the Koopman operator directly provides the osculating variationof all the variables involved in the problem, as opposed to other approacheswhere it is required to define some intermediaries or the transformation oforbital elements. Additionally, the Koopman operator allows us to handleautomatically any order of the zonal harmonics of the gravitational poten-tial, not requiring to modify the methodology in any way to obtain thesolution. Another interesting result of the Koopman operator is that themethodology provides the spectral behavior of the system, which can beused to study some properties of the problem. Finally, since the Koopmanmatrix is obtained as part of the methodology, this matrix represents thebest linear approximation of the problem with the set of basis functions andvariable ranges defined. This can be used for other applications, includingcontrol schemes. Acknowledgements
The authors wish to thank Daniel Jang at Massachusetts Institute of Tech-nology for his help in reviewing and improving the use of English in thismanuscript. Page 25onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares
References [1] M. Irigoyen and C. Simo. “Non integrability of theJ 2 problem”. In:
Celestial Mechanics and Dynamical Astronomy doi : .[2] A. Celletti and P. Negrini. “Non-integrability of the problem of mo-tion around an oblate planet”. In: Celestial Mechanics and DynamicalAstronomy doi : .[3] D. Brouwer. “Solution of the problem of artificial satellite theory with-out drag”. In: The Astronomical Journal
64 (1959), pp. 378–397.[4] Y. Kozai. “The motion of a close earth satellite”. In:
The AstronomicalJournal
64 (1959), p. 367.[5] Y. Kozai. “Second-order solution of artificial satellite theory withoutair drag”. In:
The Astronomical Journal
67 (1962), p. 446.[6] A. Deprit. “Canonical transformations depending on a small parame-ter”. In:
Celestial Mechanics doi : .[7] J. Liu. “Satellite motion about an oblate Earth”. In: AIAA Journal doi : .[8] A. Deprit. “The elimination of the parallax in satellite theory”. In: Ce-lestial Mechanics doi : .[9] A. Deprit. “The main problem in the theory of artificial satellites toorder four”. In: Journal of Guidance and Control doi : .[10] S. Coffey and A. Deprit. “Third-order solution to the main problemin satellite theory”. In: Journal of Guidance, Control, and Dynamics doi : .[11] K. T. Alfriend and S. L. Coffey. “Elimination of the perigee in thesatellite problem”. In: Celestial Mechanics
32 (1984), 163–172. doi : .[12] A. Kamel and T. Duhamel. “A second-order solution of the main prob-lem of artificial satellites using multiple scales”. In: Journal of Guid-ance, Control, and Dynamics doi : .[13] Y. Hashida and P. L. Palmer. “Epicyclic motion of satellites about anoblate planet”. In: Journal of Guidance, Control, and Dynamics doi : .[14] I. M. Ross. “Linearized dynamic equations for spacecraft subject toj perturbations”. In: Journal of guidance, control, and dynamics doi : . Page 26onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares[15] M. Humi. “J2 Effect in Cylindrical Coordinates”. In: Journal of guid-ance, control, and dynamics doi : .[16] R. Broucke. “Stability of periodic orbits in the elliptic, restricted three-body problem”. In: AIAA journal doi : .[17] C. R. McInnes. “Dynamics, stability, and control of displaced non-Keplerian orbits”. In: Journal of Guidance, Control, and Dynamics doi : .[18] D. Scheeres, M. Guman, and B. Villac. “Stability analysis of plane-tary satellite orbiters: application to the Europa orbiter”. In: Journalof Guidance, Control, and Dynamics doi : .[19] D. Arnas, D. Casanova, and E. Tresaco. “Relative and Absolute Station-Keeping for Two-Dimensional–Lattice Flower Constellations”. In: Jour-nal of Guidance, Control, and Dynamics doi : .[20] C. N. McGrath and M. Macdonald. “General perturbation method forsatellite constellation reconfiguration using low-thrust maneuvers”. In: Journal of Guidance, Control, and Dynamics doi : .[21] D. Arnas and D. Casanova. “Nominal definition of satellite constella-tions under the Earth gravitational potential”. In: Celestial Mechan-ics and Dynamical Astronomy
132 (2020), pp. 1–20. doi : .[22] J.-F. Hamel and J. de Lafontaine. “Linearized dynamics of formationflying spacecraft on a J2-perturbed elliptical orbit”. In: Journal ofGuidance, Control, and Dynamics doi : .[23] S. R. Vadali. “Model for linearized satellite relative motion about aJ2-perturbed mean circular orbit”. In: Journal of guidance, control,and dynamics doi : .[24] D. Arnas, P. Jurado, I. Barat, B. Duesmann, and R. Bock. “FLEX:A Parametric Study of Its Tandem Formation With Sentinel-3”. In: IEEE Journal of Selected Topics in Applied Earth Observations andRemote Sensing doi :
10 . 1109 / JSTARS.2019.2896196 .[25] D. Arnas. “Linearized model for satellite station-keeping and tan-dem formations under the effects of atmospheric drag”. In:
Acta As-tronautica
178 (2021), pp. 835–845. doi :
10 . 1016 / j . actaastro.2020.10.035 . Page 27onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares[26] R. Lyddane. “Small eccentricities or inclinations in the Brouwer theoryof the artificial satellite”. In:
The Astronomical Journal
68 (1963),p. 555. doi : .[27] C. Cohen and R. Lyddane. “Radius of convergence of Lie series forsome elliptic elements”. In: Celestial Mechanics doi : .[28] S. L. Coffey, A. Deprit, and B. R. Miller. “The critical inclination in ar-tificial satellite theory”. In: Celestial Mechanics doi : .[29] A. A. Kamel. “Expansion formulae in canonical transformations de-pending on a small parameter”. In: Celestial Mechanics doi : .[30] A. A. Kamel. “Perturbation method in the theory of nonlinear oscilla-tions”. In: Celestial Mechanics doi : .[31] A. Deprit and A. Rom. “The main problem of artificial satellite theoryfor small and moderate eccentricities”. In: Celestial Mechanics doi : .[32] A. Deprit. “Delaunay normalisations”. In: Celestial Mechanics doi : .[33] R. Cid, S. Ferrer, and M. Sein-Echaluce. “On the radial intermedi-aries and the time transformation in satellite theory”. In: CelestialMechanics doi : .[34] A Abad, J. San-Juan, and A Gav´ın. “Short term evolution of artificialsatellites”. In: Celestial Mechanics and Dynamical Astronomy doi : .[35] A. Abad, M. Calvo, and A. Elipe. “Integration of Deprit’s radial inter-mediary”. In: Acta Astronautica
173 (2020), pp. 19–21. doi : .[36] B. Mahajan, S. R. Vadali, and K. T. Alfriend. “Exact Delaunay nor-malization of the perturbed Keplerian Hamiltonian with tesseral har-monics”. In: Celestial Mechanics and Dynamical Astronomy doi : .[37] M. Lara. “Solution to the main problem of the artificial satellite by re-verse normalization”. In: Nonlinear Dynamics doi : .[38] J. B. Conway. A course in functional analysis . Vol. 96. Springer, 2019. isbn : 978-1-4757-3830-8. doi : .[39] K. Kowalski and W.-h. Steeb. Nonlinear dynamical systems and Car-leman linearization . World Scientific, 1991. isbn : 978-9810205874.Page 28onal Harmonics Motion Using Koopman Operator D. Arnas, R. Linares[40] I. Prigogine.
Non–equilibrium statistical mechanics . Dover Publica-tions, 1962. isbn : 978-0486815558.[41] A. Naylor and G. Sell.
Linear Operator Theory in Engineering and Sci-ence . Applied Mathematical Sciences. Springer, 2000. isbn : 9780387950013.[42] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley. “A data–drivenapproximation of the koopman operator: Extending dynamic mode de-composition”. In:
Journal of Nonlinear Science doi : .[43] I. Mezi´c. “Analysis of fluid flows via spectral properties of the Koop-man operator”. In: Annual Review of Fluid Mechanics
45 (2013),pp. 357–378. doi : .[44] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz. “Koop-man invariant subspaces and finite linear representations of nonlineardynamical systems for control”. In: PloS one doi : .[45] A. Surana and A. Banaszuk. “Linear observer synthesis for nonlinearsystems using Koopman operator framework”. In: IFAC-PapersOnLine doi : .[46] A. Surana. “Koopman operator based observer synthesis for control-affine nonlinear systems”. In: Decision and Control (CDC), 2016 IEEE55th Conference on . IEEE. 2016, pp. 6492–6499. doi : .[47] B. O. Koopman. “Hamiltonian systems and transformation in Hilbertspace”. In: Proceedings of the National Academy of Sciences doi : .[48] J. v. Neumann. “Zur Operatorenmethode in der klassischen Mechanik”.In: Annals of Mathematics (1932), pp. 587–642. doi : .[49] R. Linares. “Koopman Operator Theory Applied to the Motion ofSatellites”. In: Advances in the Astronautical Sciences . Vol. 171. AAS/AIAA.2019, AAS 19–821. isbn : 978-0-87703-665-4.[50] T. Chen and J. Shan. “Koopman-Operator-Based Attitude Dynam-ics and Control on SO (3)”. In:
Journal of Guidance, Control, andDynamics doi : .[51] D. Arnas and R. Linares. “A set of orbital elements to fully repre-sent the zonal harmonics around an oblate celestial body”. In: arXivpreprint arXiv:2010.12722arXivpreprint arXiv:2010.12722