最近对星下点的轨迹比较感兴趣，星下点的轨迹直接反应了卫星在什么情况会经过哪个地区的上空，对在全球范围内测控站的选取以及卫星观测具有重要意义，于是想在地理信息平台中将星下点轨迹画出来。于是做了一些尝试：

首先应该对卫星轨道投影在地球表面的轨道方程做相关了解，通过相关推到我们直接可以得到卫星在某时刻的投影经纬度公式：

其中φ′= arcsin ( sin nΔt ·sini)

i为轨道倾角，n为卫星平均角速度，可以根据开普勒公式轻松得到：

t 时刻卫星经度λ为

λ = λ0 +Δλ

其中　tanΔλ = (- sinω′eΔt ·cos nΔt + cosω′eΔt ·sin nΔt ·cos i) / (cosω′eΔt ·cosnΔt +sinω′eΔt ·sin nΔt ·cos i)

具体相关公式推导和参数含义可以参照相关论文

其中我们需要在地图上反应轨道投影随时间t的变化情况，采用了C#平台上的Gmap.net通过调用这款控件，结合上述公式，可以得出卫星轨迹，核心代码如下：其中代码必须考虑轨道周期性变化以及及时将经度加减360度的问题

#region星下点轨迹计算

    //  现给定一组设计的初始轨道根数

　　//a = 7 716 343. 89 m

　　//i = 66. 01°

//λ0 = 107. 010 2°      

 double miu = 3.986005 * 100000000000000;

       doubleJ2 = 1.083 * 0.001;

       doubleae = 6378137;

       publicdouble[] getorbit(double[]canshu, double t)

       {

           doublea = canshu[0];//轨¨¬道쨤长¡è半ã?轴¨￠

           doublei = canshu[1];//轨¨¬道쨤倾?角?

           doublelamdda0 = canshu[2];//轨¨¬道쨤升|y交?点ì?经-度¨¨

           doublen = Math.Sqrt(miu) / (Math.Sqrt(a * a * a));

           doublefai = Math.Asin(Math.Sin(n* t) * Math.Sin(i));//此ä?处ä|仅?仅?是o?地ì?心?纬3度¨¨

           doublewomiga = -1.5 * J2 * Math.Sqrt(miu) * ae *ae * Math.Cos(i) / Math.Pow(a,3.5);

           doubleWepie = 7.292115 * Math.Pow(10, -5) -womiga;

           doubledetalamda = (Math.Atan((-Math.Sin(Wepie * t) * Math.Cos(n* t) + Math.Cos(Wepie * t) * Math.Sin(n * t) * Math.Cos(i))/

               (Math.Cos(Wepie* t) * Math.Cos(n * t) + Math.Sin(Wepie * t) * Math.Sin(n* t) * Math.Cos(i))));

           doubledetalamdajiao = detalamda * 180 / Math.PI;

           //doublelamda =  detalamdajiao;

           doublelamda = lamdda0 + detalamdajiao;

           if(lamda > 180)

           {

               lamda = lamda - 360;

           }

           double[]JingWei = newdouble[2];

           JingWei[0] = fai;

           JingWei[1] = lamda;

           returnJingWei;

       }

       #endregion

       privatevoid button7_Click(objectsender, EventArgs e)

       {

           this.MainMap.Zoom= 1;//设|¨¨置?控?件t当ì¡à前¡ã的ì?缩?放¤?比ਨ例¤y

           double[]canshu = newdouble[3];

           canshu[0] = 7716343.89;//卫¨¤星?长¡è半ã?轴¨￠

           canshu[1] = (66.01 / 180) * Math.PI; //轨¨¬道쨤倾?角?

           canshu[2] = 120;   //升|y交?点ì?经-度¨¨

           inttime = 70000;

           double[,]res1 = newdouble[time+ 1, 2];

           for(int t = 0; t <time; t = t + 1)

           {

               double[]res = getorbit(canshu, t);

               //double[]res = GetOb(t);

               doublewei = res[0] * 180 / Math.PI;

               doublejing = res[1];

               res1[t, 0] = wei;

               res1[t, 1] = jing;

               //AddOrbit(wei,jing);

           }

           for(int t = 0; t < time; t++)

           {

               if(res1[t + 1, 1] - res1[t, 1] > 100)

               {

                   for (int j = t+1; j < time; j++)

                   {

                       res1[j, 1] = res1[j, 1]- 180;

                   }

               }

               for(int j = t; j < time; j++)

               {

                    if(res1[t, 1] < -180)

                    {

                         res1[t, 1] += 360;

                     }

               }

               AddOrbit(res1[t,0], res1[t,1]);

           }

           }

   }

如果加上timer控件，还可实现轨迹的动画效果，程序效果图如下：

（时间取了70000秒）

此外，也可用GMT生成xingxiadian星下点轨迹图，但是很难zu到动画效果，所以采用C#的地理信息平台，如果需要生成卫星轨迹的三维模型OPENGL和谷歌地球的二次开发是个不错的选择。