View source for Frozen orbit

{{Short description|Orbit in which natural drifting has been minimized}}

In [[orbital mechanics]], a '''frozen orbit''' is an [[orbit]] for an artificial [[satellite]] in which [[Perturbation (astronomy) |perturbations]] have been minimized by careful selection of the [[orbital parameters]]. Perturbations can result from natural drifting due to the [[central body]]'s shape, or other factors. Typically, the [[altitude]] of a satellite in a frozen orbit remains constant at the same point in each revolution over a long period of time.<ref name=cdeagle>{{cite web|last=Eagle |first=C. David |title=Frozen Orbit Design |url=http://www.cdeagle.com/omnum/pdf/frozen1.pdf |work=Orbital Mechanics with Numerit |accessdate=5 April 2012 |url-status=dead |archiveurl=https://web.archive.org/web/20111121072936/http://www.cdeagle.com/omnum/pdf/frozen1.pdf |archivedate=21 November 2011 }}</ref> Variations in the [[orbital inclination|inclination]], [[argument of perigee|position]] of the [[apsis]] of the orbit, and [[orbital eccentricity|eccentricity]] have been minimized by choosing [[initial value]]s so that their perturbations cancel out.<ref name=chobotov>{{cite book |last=Chobotov |first=Vladimir A. |title=Orbital Mechanics |edition=3rd |year=2002 |publisher=[[American Institute of Aeronautics and Astronautics]] |page=221 |url=http://www.knovel.com/web/portal/browse/display?_EXT_KNOVEL_DISPLAY_bookid=1586&VerticalID=0}}</ref> This results in a long-term stable orbit that minimizes the use of [[orbital station-keeping|station-keeping]] propellant.

==Background and motivation==
{{unsourced|section|date=March 2024}}
For spacecraft in orbit around the Earth, changes to orbital parameters are caused by the [[Geopotential model|oblateness of the Earth]], gravitational attraction from the Sun and Moon, [[solar radiation pressure]] and [[air drag]].<ref>{{cite book |last1=Ley |first1=Wilfried |last2=Wittmann |first2=Klaus |last3=Willi |first3=Hallmann |title=Handbuch der Raumfahrt |date=2019 |publisher=Carl Hanser Verlag München |isbn=978-3-446-45429-3 |page=109 |edition=5.}}</ref> These are called perturbing forces. They must be counteracted by maneuvers to keep the spacecraft in the desired orbit. For a [[Geostationary orbit|geostationary spacecraft]], correction maneuvers on the order of {{cvt|40–50|m/s|mph}} per year are required to counteract the gravitational forces from the Sun and Moon which move the orbital plane away from the equatorial plane of the Earth.{{cn|date=March 2024}}

For [[Sun-synchronous orbit|Sun-synchronous spacecraft]], intentional shifting of the orbit plane (called "precession") can be used for the benefit of the mission. For these missions, a near-circular orbit with an altitude of 600–900&nbsp;km is used. An appropriate inclination (97.8-99.0 degrees)<ref>{{cite book |last1=Ley |first1=Wilfried |last2=Wittmann |first2=Klaus |last3=Hallmann |first3=Willi |title=Handbuch der Raumfahrttechnik |publisher=Carl Hanser Verlag München |page=560 |edition=5}}</ref> is selected so that the precession of the orbital plane is equal to the rate of movement of the Earth around the Sun, about 1 degree per day.

As a result, the spacecraft will pass over points on the Earth that have the same time of day during every orbit. For instance, if the orbit is "square to the Sun", the vehicle will always pass over points at which it is 6 a.m. on the north-bound portion, and 6 p.m. on the south-bound portion (or vice versa). This is called a "Dawn-Dusk" orbit. Alternatively, if the Sun lies in the orbital plane, the vehicle will always pass over places where it is midday on the north-bound leg, and places where it is midnight on the south-bound leg (or vice versa). These are called "Noon-Midnight" orbits. Such orbits are desirable for many Earth observation missions such as weather, imagery, and mapping.

The perturbing force caused by the oblateness of the Earth will in general perturb not only the orbital plane but also the [[eccentricity vector]] of the orbit. There exists, however, an almost circular orbit for which there are no secular/long periodic perturbations of the eccentricity vector, only periodic perturbations with period equal to the orbital period. Such an orbit is then perfectly periodic (except for the orbital plane precession) and it is therefore called a "frozen orbit". Such an orbit is often the preferred choice for an Earth observation mission where repeated observations of the same area of the Earth should be made under as constant observation conditions as possible.

The [[Earth observation satellite]]s are often operated in Sun-synchronous frozen orbits due to the observational advantages they provide.

===Lunar frozen orbits===
; Low orbits
Through a study of many [[lunar orbit]]ing satellites, scientists have discovered that most [[Lunar orbit|low lunar orbit]]s (LLO) are unstable.<ref name=Elipe2003>[https://arc.aiaa.org/doi/abs/10.2514/2.5064?journalCode=jgcd Frozen Orbits About the Moon. 2003]</ref> Four ''frozen lunar orbits'' have been identified at 27°, 50°, 76°, and 86° inclination. NASA described this in 2006:

<blockquote>[[Lunar mascons]] make most low lunar orbits unstable ... As a satellite passes 50 or 60 miles overhead, the mascons pull it forward, back, left, right, or down, the exact direction and magnitude of the tugging depends on the satellite's trajectory. Absent any periodic boosts from onboard rockets to correct the orbit, most satellites released into low lunar orbits (under about 60 miles or 100 km) will eventually crash into the Moon. ... [There are] a number of 'frozen orbits' where a spacecraft can stay in a low lunar orbit indefinitely. They occur at four inclinations: 27°, 50°, 76°, and 86°"—the last one being nearly over the lunar poles. The orbit of the relatively long-lived Apollo 15 subsatellite [[PFS-1]] had an inclination of 28°, which turned out to be close to the inclination of one of the frozen orbits—but less fortunate [[PFS-2]] had an orbital inclination of only 11°.<ref name=nasa20061106>{{cite web |title=Bizarre Lunar Orbits |url=https://science.nasa.gov/science-news/science-at-nasa/2006/06nov_loworbit/ |last=Bell |first=Trudy E. |date=November 6, 2006 |editor-last=Phillips |editor-first=Tony |work=Science@NASA |publisher=[[NASA]] |accessdate=2017-09-08 |archive-date=2021-12-04 |archive-url=https://web.archive.org/web/20211204040014/https://science.nasa.gov/science-news/science-at-nasa/2006/06nov_loworbit/ |url-status=dead }}</ref></blockquote>

;Elliptical inclined orbits
For lunar orbits with altitudes in the {{cvt|500|to|20000|km|-1}} range, the gravity of Earth leads to orbit [[Perturbation (astronomy) |perturbations]]. Work published in 2005 showed a class of elliptical inclined lunar orbits resistant to this and are thus also frozen.<ref name=ely-2005>{{cite journal |first=Todd |last=Ely |title=Stable Constellations of Frozen Elliptical Inclined Lunar Orbits|date=July 2005 |journal=The Journal of the Astronautical Sciences |volume=53 |issue=3 |pages=301–316 |doi=10.1007/BF03546355 |bibcode=2005JAnSc..53..301E |url=https://www.researchgate.net/publication/264815005}}</ref>

==Classical theory==
{{cleanup
 |1=section
 |reason=[[Wikipedia:NOTTEXTBOOK|Wikipedia Is Not A Textbook]]. The provided derivation is lacking context for why it would be useful to readers.
|date=September 2013}}
{{cleanup|section|reason=variables and constants need to be defined for clarity|date=January 2025}}

The classical theory of frozen orbits is essentially based on the analytical [[Orbital perturbation analysis|perturbation analysis]] for artificial satellites of [[Dirk Brouwer]] made under contract with [[NASA]] and published in 1959.<ref>Dirk Brouwer: "Solution of the Problem of the Artificial Satellite Without Drag", ''Astronomical Journal'', 64 (1959)</ref>

This analysis can be carried out as follows:

In the article [[orbital perturbation analysis]], the secular perturbation of the orbital pole <math>\Delta \hat{z}\,</math> from the <math>J_2\,</math> term of the [[geopotential model]] is shown to be
{{NumBlk|::|<math>
\Delta \hat{z}\ =\ -2\pi\ \frac{J_2}{\mu\ p^2}\ \frac{3}{2}\ \cos i\ \sin i \quad \hat{g}
</math>|{{EquationRef|1}}}}

which can be expressed in terms of orbital elements thus:

{{NumBlk|::|<math>
\Delta i\ =\ 0
</math>|{{EquationRef|2}}}}
{{NumBlk|::|<math>
\Delta \Omega\ =\ -2\pi\ \frac{J_2}{\mu\ p^2}\ \frac{3}{2}\ \cos i
</math>|{{EquationRef|3}}}}

Making a similar analysis for the <math>J_3\,</math> term (corresponding to the fact that the earth is [[Figure of the Earth|slightly pear shaped]]), one gets
{{NumBlk|::|<math>
\Delta \hat{z}\ =\ 2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \cos i\ \left(\ e_h\ \left(1 - \frac{15}{4}\ \sin^2 i\right)\ \hat{g}\ -\ e_g\ \left(1 - \frac{5}{4}\ \sin^2 i\right)\ \hat{h}\right)
</math>|{{EquationRef|4}}}}

which can be expressed in terms of orbital elements as

{{NumBlk|::|<math>
\Delta i\ =\ -2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \cos i\ e_g\ \left(1 - \frac{5}{4}\ \sin^2 i\right)
</math>|{{EquationRef|5}}}}
{{NumBlk|::|<math>
\Delta \Omega\ =\ 2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \frac{\cos i}{\sin i}\ \ e_h\ \left(1 - \frac{15}{4}\ \sin^2 i\right)
</math>|{{EquationRef|6}}}}

In the same article the secular perturbation of the components of the [[Kepler orbit|eccentricity vector]] caused by the <math>J_2\,</math> is shown to be:

{{NumBlk|::|<math>
\begin{align}
  \left(\Delta e_g, \Delta e_h\right)\ =\ 
    & -2\pi\ \frac {J_2}{\mu\ p^2}\ \frac{3}{2} \left(\frac{3}{2}\ \sin^2 i\ -\ 1\right)\ (-e_h ,e_g)\ + \ 2\pi\ \frac {J_2}{\mu\ p^2}\ \frac{3}{2}\ \cos^2 i \left(-e_h, e_g\right)\ = \\
    & -2\pi\ \frac {J_2}{\mu\ p^2}\ 3 \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right)\ \left(-e_h, e_g\right)
\end{align}
</math>|{{EquationRef|7}}}}

where:

* The first term is the in-plane perturbation of the eccentricity vector caused by the in-plane component of the perturbing force
* The second term is the effect of the new position of the ascending node in the new orbital plane, the orbital plane being perturbed by the out-of-plane force component

Making the analysis for the <math>J_3\,</math> term one gets for the first term, i.e. for the perturbation of the eccentricity vector from the in-plane force component

{{NumBlk|::|<math>
2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \sin i\ \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right) 
\left(\left(1 - {e_g}^2 + 4\ {e_h}^2\right)\ \hat{g}\ -\ 5\ e_g\ e_h\ \hat{h}\right)
</math>|{{EquationRef|8}}}}

For inclinations in the range 97.8–99.0 deg, the <math>\Delta\Omega\,</math> value given by ({{EquationNote|6}}) is much smaller than the value given by ({{EquationNote|3}}) and can be ignored. Similarly the quadratic terms of the eccentricity vector components in ({{EquationNote|8}}) can be ignored for almost circular orbits, i.e. ({{EquationNote|8}}) can be approximated with

{{NumBlk|::|<math>
2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \sin i\ \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right) \ \hat{g}
</math>|{{EquationRef|9}}}}

Adding the <math>J_3\,</math> contribution 
<math>
2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \sin i\ \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right) \ (1,\ 0)
</math>

to ({{EquationNote|7}}) one gets
{{NumBlk|::|<math>
\left(\Delta e_g, \Delta e_h\right)\ =\ -2\pi\ \frac {J_2}{\mu\ p^2}\ 3 \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right)\ \left(-\left(e_h+\frac{J_3\ \sin i}{J_2\ 2\ p}\right),\ e_g \right) 
</math>|{{EquationRef|10}}}}

Now the difference equation shows that the eccentricity vector will describe a circle centered at the point <math>\left(\ 0,\ -\frac{J_3\ \sin i}{J_2\ 2\ p}\ \right)\,</math>; the polar argument of the eccentricity vector increases with <math> -2\pi\ \frac {J_2}{\mu\ p^2}\ 3 \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right)\,</math> radians between consecutive orbits.

As
:<math>\mu = 398600.440\text{ km}^3/s^2 \, </math>
:<math>J_2 = 1.7555\ 10^{10}\text{ km}^5/s^2 \, </math>
:<math>J_3 = -2.619\ 10^{11}\text{ km}^6/s^2 \, </math>

one gets for a polar orbit (<math> i = 90^\circ\,</math>) with <math> p = 7200\text{ km}\,</math> that the centre of the circle is at <math>(0,\ 0.001036)\,</math> and the change of polar argument is 0.00400 radians per orbit.

The latter figure means that the eccentricity vector will have described a full circle in 1569 orbits.
Selecting the initial mean eccentricity vector as <math> (0,\ 0.001036)\,</math> the mean eccentricity vector will stay constant for successive orbits, i.e. the orbit is frozen because the secular perturbations of the <math>J_2\,</math> term given by ({{EquationNote|7}}) and of the <math>J_3\,</math> term given by ({{EquationNote|9}}) cancel out.

In terms of classical orbital elements, this means that a frozen orbit should have the following mean elements:
:<math>e = -\frac{J_3\ \sin i}{J_2\ 2\ p}\,</math>
:<math>\omega =\ 90^\circ\,</math>

==Modern theory==
The modern theory of frozen orbits is based on the algorithm given in a 1989 article by Mats Rosengren.<ref>{{cite conference |bibcode=1989ommd.proc...49R |author=Mats Rosengren |title=Improved technique for Passive Eccentricity Control (AAS 89-155) |volume=69 |publisher=AAS/NASA |book-title=Advances in the Astronautical Sciences|year=1989}}</ref>

For this the analytical expression ({{EquationNote|7}}) is used to iteratively update the initial (mean) eccentricity vector to obtain that the (mean) eccentricity vector several orbits later computed by the precise numerical propagation takes precisely the same value. In this way the secular perturbation of the eccentricity vector caused by the <math>J_2\,</math> term is used to counteract all secular perturbations, not only those (dominating) caused by the <math>J_3\,</math> term. One such additional secular perturbation that in this way can be compensated for is the one caused by the [[solar radiation pressure]], this perturbation is discussed in the article "[[Orbital perturbation analysis (spacecraft)]]".

Applying this algorithm for the case discussed above, i.e. a polar orbit (<math> i=90^\circ\,</math>) with <math> p=7200\text{ km}\,</math> ignoring all perturbing forces other than the <math>J_2\,</math> and the <math>J_3\,</math> forces for the numerical propagation one gets exactly the same optimal average eccentricity vector as with the "classical theory", i.e. <math> (0,\ 0.001036)\,</math>.

When we also include the forces due to the higher zonal terms the optimal value changes to <math> (0,\ 0.001285)\,</math>.

Assuming in addition a reasonable solar pressure (a "cross-sectional-area" of {{nowrap|0.05 m<sup>2</sup>/kg}}, the direction to the sun in the direction towards the ascending node) the optimal value for the average eccentricity vector becomes <math> (0.000069,\ 0.001285)\,</math> which corresponds to :<math>\omega =\ 87^\circ\,</math>, i.e. the optimal value is '''not''' <math>\omega =\ 90^\circ</math> anymore.

This algorithm is implemented in the orbit control software used for the [[Earth observation satellite]]s [[European Remote-Sensing Satellite|ERS-1, ERS-2]] and [[Envisat]]

==Derivation of the closed form expressions for the ''J''<sub>3</sub> perturbation==
{{technical
 |reason=section is relevant as it explains why frozen orbits exist, but may be difficult for laymen to understand. The mathematics should be retained, but with additional information and explanation. 
|date=September 2013}}

The main perturbing force to be counteracted in order to have a frozen orbit is the "<math>J_3\,</math> force", i.e. the gravitational force caused by an imperfect symmetry north–south of the Earth, and the "classical theory" is based on the closed form expression for this "<math>J_3\,</math> perturbation".  With the "modern theory" this explicit closed form expression is not directly used but it is certainly still worthwhile{{For whom|date=March 2024}} to derive it.
The derivation of this expression can be done as follows:

The potential from a zonal term is rotational symmetric around the  polar axis of the Earth and corresponding force is entirely in a longitudinal plane with one component <math>F_r\ \hat{r}\,</math> in the radial direction and one component <math>F_\lambda\ \hat{\lambda}\,</math>  with the unit vector <math>\hat{\lambda}\,</math> orthogonal to the radial direction towards north. These directions <math>\hat{r}\,</math> and <math>\hat{\lambda}\,</math> are illustrated in Figure 1.

[[File:Spherical coordinates unit vectors.svg|thumb|right|Figure 1: The unit vectors <math>\hat{\phi},\ \hat{\lambda},\ \hat{r}</math>]]

In the article [[Geopotential model]] it is shown that these force components caused by the <math>J_3\,</math> term are

{{NumBlk|::|<math>
\begin{align}
&F_r = J_3\ \frac{1}{r^5}\ 2\ \sin\lambda\ \left(5\sin^2\lambda\ -\ 3\right) \\
&F_\lambda = -J_3\ \frac{1}{r^5}\ \frac{3}{2}\ \cos\lambda\ \left(5\ \sin^2\lambda\ -1\right)
\end{align}
</math>|{{EquationRef|11}}}}

To be able to apply relations derived in the article [[Orbital perturbation analysis (spacecraft)]] the force component <math>F_\lambda\ \hat{\lambda}\,</math> must be split into two orthogonal components <math>F_t\ \hat{t}</math> and <math>F_z\ \hat{z}</math> as illustrated in figure 2

[[File:Zonal term force components.svg|thumb|right|Figure 2: The unit vector <math>\hat{t}\,</math> orthogonal to <math>\hat{r}\,</math> in the direction of motion and the orbital pole <math>\hat{z}\,</math>. The force component <math>F_\lambda</math> is marked as "F"]]
Let <math>\hat{a},\ \hat{b},\ \hat{n}\,</math> make up a rectangular coordinate system with origin in the center of the Earth (in the center of the [[Reference ellipsoid]])  such that <math>\hat{n}\,</math> points in the direction north and such that <math>\hat{a},\ \hat{b}\,</math> are in the equatorial plane of the Earth with <math>\hat{a}\,</math> pointing towards the [[Orbital node|ascending node]], i.e. towards the blue point of Figure 2.

The components of the unit vectors
:<math>\hat{r},\ \hat{t},\ \hat{z}\,</math>

making up the local coordinate system (of which <math>\hat{t},\ \hat{z},</math> are illustrated in figure 2), and expressing their relation with <math>\hat{a},\ \hat{b},\ \hat{n}\,</math>, are as follows:

:<math>r_a= \cos u\,</math>
:<math>r_b= \cos i \ \sin u\,</math>
:<math>r_n= \sin i \ \sin u\,</math>
:<math>t_a=-\sin u\,</math>
:<math>t_b= \cos i \ \cos u\,</math>
:<math>t_n= \sin i \ \cos u\,</math>
:<math>z_a= 0\,</math>
:<math>z_b=-\sin i\,</math>
:<math>z_n= \cos i\,</math>

where <math>u\,</math> is the polar argument of <math>\hat{r}\,</math> relative the orthogonal unit vectors <math>\hat{g}=\hat{a}\,</math> and <math>\hat{h}=\cos i\ \hat{b}\ +\ \sin i\ \hat{n}\,</math> in the orbital plane

Firstly

:<math>\sin \lambda =\ r_n\ =\ \sin i \ \sin u\,</math>

where  <math>\lambda\,</math> is the angle between the equator plane and <math>\hat{r}\,</math> (between the green points of figure 2) and from equation (12) of the article [[Geopotential model]] one therefore obtains
{{NumBlk|::|<math>
F_r = J_3\ \frac{1}{r^5}\ 2\ \sin i \ \sin u\,\ \left(5\sin^2 i \ \sin^2 u\ -\ 3\right)
</math>|{{EquationRef|12}}}}

Secondly the projection of direction north, <math>\hat{n}\,</math>, on the plane spanned by <math>\hat{t},\ \hat{z},</math> is

:<math>\sin i \ \cos u \ \hat{t}\ +\ \cos i \ \hat{z}\,</math>

and this projection is 
:<math>\cos \lambda \ \hat{\lambda}\,</math>

where <math>\hat{\lambda}\,</math> is the unit vector <math>\hat{\lambda}</math> orthogonal to the radial direction towards north illustrated in figure 1.

From equation ({{EquationNote|11}}) we see that

:<math>F_\lambda \ \hat{\lambda}\ = -J_3\ \frac{1}{r^5}\ \frac{3}{2}\ \left(5\ \sin^2\lambda\ -1\right)\ \cos\lambda\ \hat{\lambda}\ =\ -J_3\ \frac{1}{r^5}\ \frac{3}{2}\ \left(5\ \sin^2\lambda\ -1\right)\ (\sin i \ \cos u \ \hat{t}\ +\ \cos i \ \hat{z})\,</math>

and therefore:
{{NumBlk|::|<math>
F_t  =\ -J_3\ \frac{1}{r^5}\ \frac{3}{2}\ \left(5\ \sin^2 i \ \sin^2 u\ -1\right)\ \sin i \ \cos u
</math>|{{EquationRef|13}}}}
{{NumBlk|::|<math>
F_z =\ -J_3\ \frac{1}{r^5}\ \frac{3}{2}\ \left(5\ \sin^2 i \ \sin^2 u\ -1\right)\ \cos i
</math>|{{EquationRef|14}}}}

In the article [[Orbital perturbation analysis (spacecraft)]] it is further shown that the secular perturbation of the orbital pole <math>\hat{z}\,</math> is

{{NumBlk|::|<math>
\Delta \hat{z}\ =\ \frac{1}{\mu p}\left[\hat{g}\int\limits_{0}^{2\pi}F_z r^3 \cos u \ du
+\ \hat{h}\int\limits_{0}^{2\pi}F_z r^3 \sin u \ du \right]\quad \times \ \hat{z}
</math>|{{EquationRef|15}}}}

Introducing the expression for <math>F_z\,</math> of ({{EquationNote|14}}) in ({{EquationNote|15}}) one gets
{{NumBlk|::|<math>
\begin{align}
&\Delta \hat{z}\ = -\frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \cos i\ \cdot \\
&\left[ \hat{g}\int\limits_{0}^{2\pi}{\left(\frac{p}{r}\right)}^2\left(5\ \sin^2 i \ \sin^2 u\ -1\right) \cos u \ du\ + 
\hat{h}\int\limits_{0}^{2\pi}{\left(\frac{p}{r}\right)}^2\left(5\ \sin^2 i \ \sin^2 u\ -1\right) \sin u \ du \right]\quad \times \ \hat{z} 
\end{align}
</math>|{{EquationRef|16}}}}

The fraction <math>\frac{p}{r}\,</math> is
:<math>\frac {p}{r}\ =\ 1 + e \cdot \cos \theta\ =\ 1 + e_g \cdot \cos u + e_h \cdot \sin u</math>

where

:<math>e_g =\ e\ \cos \omega</math>
:<math>e_h =\ e\ \sin \omega</math>

are the components of the eccentricity vector in the <math>\hat{g},\ \hat{h}\,</math> coordinate system.

As all integrals of type
:<math>\int\limits_{0}^{2\pi} \cos^m u \ \sin^n u\ du\,</math>
are zero if not both <math>n\,</math> and <math>m\,</math> are even, we see that

{{NumBlk|::|<math>
\begin{align}
\int\limits_{0}^{2\pi}{\left(\frac{p}{r}\right)}^2\left(5\ \sin^2 i \ \sin^2 u\ -1\right) \cos u \ du 
&=\ 2\ e_g\ \left(5\ \sin^2 i \ \int\limits_{0}^{2\pi} \sin^2 u \cos^2 u \ du \ -
\int\limits_{0}^{2\pi} \cos^2 u \ du \right) \\
&=\ 2\pi\ e_g \ (\frac{5}{4} \sin^2 i-1)
\end{align}
</math>|{{EquationRef|17}}}}

and

{{NumBlk|::|<math>
\begin{align}
\int\limits_{0}^{2\pi}{\left(\frac{p}{r}\right)}^2\left(5\ \sin^2 i \ \sin^2 u\ -1\right) \sin u \ du
&=\ 2\ e_h\ \left(5\ \sin^2 i \ \int\limits_{0}^{2\pi} \sin^4 u \ du \ -
\int\limits_{0}^{2\pi} \sin^2 u \ du \right) \\
&=\ 2\pi\ e_h \ (\frac{15}{4} \sin^2 i-1)
\end{align}
</math>|{{EquationRef|18}}}}

It follows that

{{NumBlk|::|<math>
\begin{align}
\Delta \hat{z}\ &=\  2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \cos i\ \left[e_g \ (1-\frac{5}{4} \sin^2 i)\ \hat{g}
+\ e_h \ (1-\frac{15}{4} \sin^2 i)\ \hat{h}\right]\quad \times \ \hat{z} \\
&=\ 2\pi\ \frac{J_3}{\mu\ p^3}\ \frac{3}{2}\ \cos i\ \left[
\ e_h \ (1-\frac{15}{4} \sin^2 i)\ \hat{g}\ - e_g \ (1-\frac{5}{4} \sin^2 i)\ \hat{h}\right]
\end{align}
</math>|{{EquationRef|19}}}}

where

:<math>\hat{g}\,</math> and <math>\hat{h}\,</math> are the base vectors of the rectangular coordinate system in the plane of the reference Kepler orbit with <math>\hat{g}\,</math> in the equatorial plane towards the ascending node and <math>u\,</math> is the polar argument relative this equatorial coordinate system

:<math>f_z\,</math> is the force component (per unit mass) in the direction of the orbit pole <math>\hat{z}\,</math>

In the article [[Orbital perturbation analysis (spacecraft)]] it is shown that the secular perturbation of the eccentricity vector is

{{NumBlk|::|<math>
\Delta \bar{e}\ =\frac {1}{\mu}\ \int\limits_{0}^{2\pi}\left(-\hat{t}\ f_r\ + \ \left(2\ \hat{r}-\frac{V_r}{V_t}\ \hat{t}\right)\ f_t\right)\ r^2\ du
</math>|{{EquationRef|20}}}}

where

*<math>\hat{r},\hat{t}\,</math> is the usual local coordinate system with unit vector <math>\hat{r}\,</math> directed away from the Earth
*<math>V_r = \sqrt{\frac {\mu}{p}} \cdot e \cdot \sin \theta</math> - the velocity component in direction <math>\hat{r}\,</math>
*<math>V_t = \sqrt{\frac {\mu}{p}} \cdot (1 + e \cdot \cos \theta)</math> - the velocity component in direction <math>\hat{t}\,</math>

Introducing the expression for <math>F_r,\ F_t\,</math> of ({{EquationNote|12}}) and ({{EquationNote|13}}) in ({{EquationNote|20}}) one gets
{{NumBlk|::|<math>
\begin{align}
&\Delta \bar{e}\ =\frac {J_3}{\mu\ p^3}\ \sin i\  \cdot \\ 
&\int\limits_{0}^{2\pi}\left(-\hat{t}\ {\left(\frac{p}{r}\right)}^3\ 2 \ \sin u\,\ \left(5\sin^2 i \ \sin^2 u\ -\ 3\right)\ - \ \left(2\ \hat{r}-\frac{V_r}{V_t}\ \hat{t}\right)\ {\left(\frac{p}{r}\right)}^3\ \frac{3}{2}\ \left(5\ \sin^2 i \ \sin^2 u\ -1\right) \ \cos u\right)du
\end{align}
</math>|{{EquationRef|21}}}}

Using that

:<math>\frac {V_r}{V_t} = \frac {e_g \cdot \sin u\ -\ e_h \cdot \cos u}{\frac {p}{r}}</math>

the integral above can be split in 8 terms:

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi}\left(-\hat{t}\ {\left(\frac{p}{r}\right)}^3\ 2 \ \sin u\,\ \left(5\sin^2 i \ \sin^2 u\ -\ 3\right)\ - \ \left(2\ \hat{r}-\frac{V_r}{V_t}\ \hat{t}\right)\ {\left(\frac{p}{r}\right)}^3\ \frac{3}{2}\ \left(5\ \sin^2 i \ \sin^2 u\ -1\right) \ \cos u\right)\ du\ = \\
&-10\sin^2 i \ \int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^3\ \sin^3 u\ du \\
&+6\ \int\limits_{0}^{2\pi}\hat{t}\ {\left(\frac{p}{r}\right)}^3\ \sin u\ du \\
&-15\ \sin^2 i \int\limits_{0}^{2\pi} \hat{r}\ {\left(\frac{p}{r}\right)}^3\ \sin^2 u  \ \cos u\ du \\
&+3\ \int\limits_{0}^{2\pi} \hat{r}\ {\left(\frac{p}{r}\right)}^3\ \cos u\ du  \\
&+\frac{15}{2}\sin^2 i\ e_g\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \ \ \sin^3 u \ \cos u\ du \\
&-\frac{15}{2}\sin^2 i\ e_h \ \int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \ \ \sin^2 u \ \cos^2 u\ du \\
&-\frac{3}{2}\ e_g \ \int\limits_{0}^{2\pi} \ \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \ \sin u \ \cos u\ du \\
&+\frac{3}{2}\ e_h \ \int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \ \cos^2 u\ du
\end{align}
</math>|{{EquationRef|22}}}}

Given that
:<math>\hat{r}=\cos u\ \hat{g}\ +\ \sin u\ \hat{h}</math>
:<math>\hat{t}=-\sin u\ \hat{g}\ +\ \cos u\ \hat{h}</math>

we obtain

:<math>\frac {p}{r}\ =\ 1 + e \cdot \cos \theta\ =\ 1 + e_g \cdot \cos u + e_h \cdot \sin u</math>

and that all integrals of type

:<math>\int\limits_{0}^{2\pi} \cos^m u \ \sin^n u\ du\,</math>
are zero if not both <math>n\,</math> and <math>m\,</math> are even:

Term 1

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^3\ \sin^3 u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin^4 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin^3 u\ \cos u \ du\  = \\
&-\hat{g}\ \left(\int\limits_{0}^{2\pi}\ \sin^4 u \ du\ +\ 
3\ {e_g}^2\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^4 u \ du\ \ +\ 
3\ {e_h}^2\ \int\limits_{0}^{2\pi}\ \sin^6 u \ du\ \right) \\
&+\hat{h}\ 6\ e_g\ e_h\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^4 u \ du = \\
&-\hat{g}\ \left(2\pi \left(\frac{3}{8}\ +\ \frac{3}{16}\ {e_g}^2\ +\ \frac{15}{16}\ {e_h}^2\right)\right)
+\hat{h}\ \left(2\pi \left(\frac{3}{8}\ e_g\ e_h\right)\right)
\end{align}
</math>|{{EquationRef|23}}}}

Term 2

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^3\ \sin u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin^2 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin u\ \cos u \ du\  = \\
&-\hat{g}\ \left(\int\limits_{0}^{2\pi}\ \sin^2 u \ du\ +\ 
3\ {e_g}^2\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^2 u \ du\ \ +\ 
3\ {e_h}^2\ \int\limits_{0}^{2\pi}\ \sin^6 u \ du\ \right) \\
&+\hat{h}\ 6\ e_g\ e_h\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^2 u \ du = \\
&-\hat{g}\ \left(2\pi \left(\frac{1}{2}\ +\ \frac{3}{8}\ {e_g}^2\ +\ \frac{9}{8}\ {e_h}^2\right)\right)
+\hat{h}\ \left(2\pi \left(\frac{3}{4}\ e_g\ e_h\right)\right)
\end{align}
</math>|{{EquationRef|24}}}}

Term 3

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{r}\ {\left(\frac{p}{r}\right)}^3\ \sin^2 u\ \cos u\ du\ = 
\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin^2 u \cos^2 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin^3 u\ \cos u \ du\  = \\
&\hat{g}\ \left(\int\limits_{0}^{2\pi}\ \sin^2 u \cos^2 u \ du\ +\ 
3\ {e_g}^2\ \int\limits_{0}^{2\pi}\ \cos^4 u\ \sin^2 u \ du\ \ +\ 
3\ {e_h}^2\ \int\limits_{0}^{2\pi}\ \sin^4 u\ \cos^2 u \ du\ \right) \\
&+\hat{h}\ 6\ e_g\ e_h\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^4 u \ du = \\
&\hat{g}\ \left(2\pi \left(\frac{1}{8}\ +\ \frac{3}{16}\ {e_g}^2\ +\ \frac{3}{16}\ {e_h}^2\right)\right)
+\hat{h}\ \left(2\pi \left(\frac{3}{8}\ e_g\ e_h\right)\right)
\end{align}
</math>|{{EquationRef|25}}}}

Term 4

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{r}\ {\left(\frac{p}{r}\right)}^3\ \cos u\ du\ = 
\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \cos^2 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^3\ \sin u\ \cos u \ du\  = \\
&\hat{g}\ \left(\int\limits_{0}^{2\pi} \cos^2 u \ du\ +\ 
3\ {e_g}^2\ \int\limits_{0}^{2\pi}\ \cos^4 u \ du\ \ +\ 
3\ {e_h}^2\ \int\limits_{0}^{2\pi}\ \sin^2 u\ \cos^2 u \ du\ \right) \\
&+\hat{h}\ 6\ e_g\ e_h\ \int\limits_{0}^{2\pi}\ \cos^2 u\ \sin^2 u \ du = \\
&\hat{g}\ \left(2\pi \left(\frac{1}{2}\ +\ \frac{9}{8}\ {e_g}^2\ +\ \frac{3}{8}\ {e_h}^2\right)\right)
+\hat{h}\ \left(2\pi \left(\frac{3}{4}\ e_g\ e_h\right)\right)
\end{align}
</math>|{{EquationRef|26}}}}

Term 5

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \sin^3 u\ \cos u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin^4 u\ \cos u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin^3 u\ \cos^2 u \ du\  = \\
&-\hat{g}\ 2\ e_g\ \int\limits_{0}^{2\pi}\ \sin^4 u\ \cos^2 u \ du 
+\hat{h}\ 2\ e_h\ \int\limits_{0}^{2\pi}\ \sin^4 u\ \cos^2 u \ du =
-\hat{g}\ \left(2\pi \frac{1}{8}\ e_g\right) + \hat{h}\ \left(2\pi \frac{1}{8}\ e_h\right)
\end{align}
</math>|{{EquationRef|27}}}}

Term 6

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \sin^2 u\ \cos^2 u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin^3 u\ \cos^2 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin^2 u\ \cos^3 u \ du\  = \\
&-\hat{g}\ 2\ e_h\ \int\limits_{0}^{2\pi}\ \sin^4 u\ \cos^2 u \ du 
+\hat{h}\ 2\ e_g\ \int\limits_{0}^{2\pi}\ \sin^2 u\ \cos^4 u \ du =
-\hat{g}\ \left(2\pi \frac{1}{8}\ e_h\right) + \hat{h}\ \left(2\pi \frac{1}{8}\ e_g\right)
\end{align}
</math>|{{EquationRef|28}}}}

Term 7

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \sin u\ \cos u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin^2 u\ \cos u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin u\ \cos^2 u \ du\  = \\
&-\hat{g}\ 2\ e_g\ \int\limits_{0}^{2\pi}\ \sin^2 u\ \cos^2 u \ du 
+\hat{h}\ 2\ e_h\ \int\limits_{0}^{2\pi}\ \sin^2 u\ \cos^2 u \ du =
-\hat{g}\ \left(2\pi \frac{1}{4}\ e_g\right) + \hat{h}\ \left(2\pi \frac{1}{4}\ e_h\right)
\end{align}
</math>|{{EquationRef|29}}}}

Term 8

{{NumBlk|::|<math>
\begin{align}
&\int\limits_{0}^{2\pi} \hat{t}\ {\left(\frac{p}{r}\right)}^2\ \cos^2 u\ du\ = 
-\hat{g}\ \int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \sin u\ \cos^2 u \ du\  
+\hat{h}\int\limits_{0}^{2\pi}\ {\left(\frac{p}{r}\right)}^2\ \cos^3 u \ du\  = \\
&-\hat{g}\ 2\ e_h\ \int\limits_{0}^{2\pi}\ \sin^2 u\ \cos^2 u \ du 
+\hat{h}\ 2\ e_g\ \int\limits_{0}^{2\pi}\ \cos^4 u \ du =
-\hat{g}\ \left(2\pi \frac{1}{4}\ e_h\right) + \hat{h}\ \left(2\pi \frac{3}{4}\ e_g\right)
\end{align}
</math>|{{EquationRef|30}}}}

As

{{NumBlk|::|<math>
\begin{align}
&-10\sin^2 i \ \left(-\hat{g}\ \left(\frac{3}{8}\ +\ \frac{3}{16}\ {e_g}^2\ +\ \frac{15}{16}\ {e_h}^2\right)
+\hat{h}\ \left(\frac{3}{8}\ e_g\ e_h\right)\right) \\
&+6\ \left(-\hat{g}\ \left(\frac{1}{2}\ +\ \frac{3}{8}\ {e_g}^2\ +\ \frac{9}{8}\ {e_h}^2\right)
+\hat{h}\ \left(\frac{3}{4}\ e_g\ e_h\right)\right) \\
&-15\sin^2 i \ \left(\hat{g}\ \left(\frac{1}{8}\ +\ \frac{3}{16}\ {e_g}^2\ +\ \frac{3}{16}\ {e_h}^2\right)
+\hat{h}\ \left(\frac{3}{8}\ e_g\ e_h\right)\right) \\
&+3\left(\hat{g}\ \left(\frac{1}{2}\ +\ \frac{9}{8}\ {e_g}^2\ +\ \frac{3}{8}\ {e_h}^2\right)
+\hat{h}\ \left(\frac{3}{4}\ e_g\ e_h\right)\right) \\
&+\frac{15}{2}\sin^2 i\ e_g \ \left(-\hat{g}\ \left(\frac{1}{8}\ e_g\right) + \hat{h}\ \left( \frac{1}{8}\ e_h\right)\right) \\
&-\frac{15}{2}\sin^2 i\ e_h \ \left(-\hat{g}\ \left(\frac{1}{8}\ e_h\right) + \hat{h}\ \left( \frac{1}{8}\ e_g\right)\right) \\
&-\frac {3}{2}\         e_g \ \left(-\hat{g}\ \left(\frac{1}{4}\ e_g\right) + \hat{h}\ \left(\frac{1}{4}\ e_h\right)\right) \\
&+\frac{3}{2}\ e_h \ \left(-\hat{g}\ \left(\frac{1}{4}\ e_h\right) + \hat{h}\ \left(\frac{3}{4}\ e_g\right)\right) = \\
&\frac{3}{2}\ \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right)\left((1-{e_g}^2\ +\ 4\ {e_h}^2)\hat{g}\ -\ 5\ e_g\ e_h\ \hat{h}\right)
\end{align}
</math>|{{EquationRef|31}}}}

It follows that

{{NumBlk|::|<math>
\begin{align}
&\Delta \bar{e}\ =\frac {J_3}{\mu\ p^3}\ \sin i\  \cdot \\ 
&\int\limits_{0}^{2\pi}\left(-\hat{t}\ {\left(\frac{p}{r}\right)}^3\ 2 \ \sin u\,\ \left(5\sin^2 i \ \sin^2 u\ -\ 3\right)\ - \ \left(2\ \hat{r}-\frac{V_r}{V_t}\ \hat{t}\right)\ {\left(\frac{p}{r}\right)}^3\ \frac{3}{2}\ \left(5\ \sin^2 i \ \sin^2 u\ -1\right) \ \cos u\right)du = \\
&2\pi\ \frac {J_3}{\mu\ p^3}\ \sin i\ \frac{3}{2}\ \left(\frac{5}{4}\ \sin^2 i\ -\ 1\right)\left((1-{e_g}^2\ +\ 4\ {e_h}^2)\hat{g}\ -\ 5\ e_g\ e_h\ \hat{h}\right)
\end{align}
</math>|{{EquationRef|32}}}}

==References==
<references />

==Further reading==
* [https://web.archive.org/web/20101129002858/http://www.aeroespacial.org.br/jaesa/editions/repository/v01/n02/3-PardalEtal.pdf STUDY OF ORBITAL ELEMENTS ON THE NEIGHBOURHOOD OF A FROZEN ORBIT.] including atmospheric drag and J<sub>5</sub> terms

{{Portal bar|Astronomy|Stars|Spaceflight|Outer space|Solar System}}

{{DEFAULTSORT:Frozen Orbit}}
[[Category:Astrodynamics]]
[[Category:Orbits]]
[[Category:Spaceflight concepts]]