## Appendix A  
The following forms of the rate of six orbit elements are derived from [PERSON]'s planetary equations for the averaged geopotential ([PERSON], 1979). They can be used for solutions of the orbital parameters for repeat ground track (RGT) orbits.  
\[\dot{a}=0\] (A1) \[\dot{e}=n\frac{1-e^{2}}{e}\sum\limits_{l=2}^{\infty}\sum\limits_{m= l_{s}}^{l-2}\sum\limits_{k=m}^{l-2}\left(\frac{R_{e}}{p}\right)^{l}C_{l,0}F_{l,m}( cosi)G_{l-1,k}(e)m\overline{l}_{mk}simn\left(\omega-\frac{\pi}{2}\right)\] (A2) \[\dot{i}=-n\frac{cosi}{sini}\sum\limits_{l=2}^{\infty}\sum\limits _{m=l_{s}}^{l-2}\sum\limits_{k=m}^{l-2}\left(\frac{R_{e}}{p}\right)^{l}C_{l,0}F _{l,m}(cosi)G_{l-1,k}(e)m\overline{l}_{mk}simn\left(\omega-\frac{\pi}{2}\right)\] (A3) \[\dot{\Omega}=-n\sum\limits_{l=2}^{\infty}\sum\limits_{m=l_{s}}^{ l-2}\sum\limits_{k=m}^{l-2}\left(\frac{R_{e}}{p}\right)^{l}C_{l,0}F^{l}_{l,m}( cosi)G_{l-1,k}(e)\overline{l}_{mk}cosm\left(\omega-\frac{\pi}{2}\right)\] (A4) \[\dot{\omega}=n\sum\limits_{l=2}^{\infty}\sum\limits_{m=l_{s}}^{ l-2}\sum\limits_{k=m}^{l-2}\left(\frac{R_{e}}{p}\right)^{l}C_{l,0}G_{l-1,k}(e) \overline{l}_{mk}cosm\left(\omega-\frac{\pi}{2}\right)\left\{\left[\frac{k}{ e^{2}}+(2l-k-1)\right]F_{l,m}(cosi)\right.\] (A5) \[\left.+cosiF^{\prime}{}_{l,m}(cosi)\right\}\] \[\dot{M}=n+n\sqrt{1-e^{2}}\sum\limits_{l=2}^{\infty}\sum\limits_{m =l_{s}}^{l-2}\sum\limits_{k=m}^{l-2}\left(\frac{R_{e}}{p}\right)^{l}C_{l,0}F_{ l,m}(cosi)G_{l-1,k}(e)\overline{l}_{mk}\left[(k+3)-\frac{k}{e^{2}} \right]cosm\left(\omega-\frac{\pi}{2}\right)\] (A6) \[F_{l,m}(cosi)=(2-\delta_{0,m})\frac{(l-m)!}{(l+m)!}P_{l,m}(0)P_{l, m}(cosi)\] (A7)\[G_{l-1,k}(e)=\left(\begin{array}{c}l-1\\ k\end{array}\right)e^{k} \tag{10}\]  
\[\overline{I}_{mk}=\int_{0}^{2\pi}\left(cos^{k}fcosmf\right)df \tag{11}\]  
\(n\) is the mean motion; \(R_{e}\) is Earth's equatorial radius; \(p\) is semi-latus axis; \(C_{l,0}\) geopotential coefficient of degree \(l\) and order 0; \(P_{l,m}\) is associated Legendre polynomial degree \(l\) and order m.  
In addition, a frozen orbit limits the variations of the orbital elements in the average sense. So, a frozen orbit is also acquired for a RGT orbit. For a frozen orbit, the rate of the semi-major axis, eccentricity, inclination and argument of perigee should be design to keep zero. When \(\omega\) can be set to be \(\frac{\pi}{2}\), the rate of the four orbit elements except \(\dot{\omega}\) are all equal to 0. In order to let \(\frac{d\omega}{dt}\) equal to 0, we should adjust the value of e and let \(\frac{d\omega}{dt}\) equal to 0. The method is described in details in the [39].  
Equation (1) can be rewritten as:  
\[\dot{M}+\dot{\omega}+\frac{N}{D}\dot{\Omega}=\frac{N}{D}\omega_{E} \tag{12}\]  
From Equations (11)-(12), we can know that \(\dot{M},\ \dot{\omega}\) and \(\dot{\Omega}\) have the same factor \(n\) (\(n=\sqrt{\mu/a^{3}}\)), so the following equation can be obtained.  
\[n\left(\frac{\dot{M}}{n}+\frac{\dot{\omega}}{n}+\frac{N}{D}\frac{\dot{\Omega} }{n}\right)=\frac{N}{D}\omega_{E} \tag{13}\]  
\(a\) can be expressed as:  
\[a=\left\{\sqrt{\mu}\frac{D}{N}\frac{\dot{M}/n+\dot{\omega}/n+\frac{N}{D}\dot{ \Omega}/n}{\omega_{E}}\right\}^{2/3} \tag{14}\]  
When \(\dot{\omega}=0\), Equation (14) can be simplified as:  
\[a=\left\{\sqrt{\mu}\frac{D}{N}\frac{\dot{M}/n+\frac{N}{D}\dot{\Omega}/n}{ \omega_{E}}\right\}^{2/3} \tag{15}\]  
Both \(\dot{M}\) and \(\dot{\Omega}\) are the function of \(a\), so we can start an iterative solution of Equation (15). The algorithm starts with the following initial guess for the required mean semi-major axis:  
\[a_{0}=\mu^{1/3}\left\{\frac{N}{D}\omega_{E}\right\}^{-2/3} \tag{16}\]