GPS卫星运动及定位matlab仿真.毕业设计说明

.作者：专业班级：摘要指导老师：导航定位、定时和测速系统，现在在全球很多领域获得了应用。高，所以对卫星运动做了理想化处理，摄动力对卫星的影响忽略不计(所以为略分析，并使用matlab编程仿真实现了卫星的运功轨道平面、运动动态、可见卫星的分布及利用可见卫星计算出用户位置。.Abstract：Globalpositioningsystemisaglobal,versatility,all-weatheradvantageofnavigationandpositioning,timingandspeedsystem,nowtherehasmanyapplicationinmanyfields.GPSsatellitepositioningisacomplexsystem,whichincludesmanyparameters,suchastimeandspacecoordinatessystem.Thisdesignisisnothigh,sotodotheidealizedsatellitemovement,andignorethedisturbedmotion(socallitnon-disturbedmotion).UsingtheKeplerandleast-squaremethodforcalculatingtheparametersoforbitalmotion,forthecharacteristicsofmotiontomakeasimpleanalysis,andusethematlabsimulationtoprogramachievetheorbitalplaneofsatellite,thedynamicmotion,thedistributionofvisiblesatellitesandusingvisiblesatellitestocalculatehomeThroughthedesignhaveprimaryunderstandingfortheGPSsatellite,andunderstandingthestaticsingle-point,pseudorangeandsoon.b.目录第一章前言 1 1 .3GPS前景 2第二章GPS测量原理 42.1伪距测量的原理 42.1.1计算卫星位置 52.1.2用户位置的计算 52.1.3最小二乘法介绍 52载波相位测量原理 6第三章GPS的坐标、时间系统 10 3.1.2地球坐标系 123.2时间系统 133.2.1世界时系统 143.2.2原子时系统 15 3.2.5GPS时间系统 16 . 8 9 0 21 致.................................................................40 42.第一章前言1.1课题背景GPS军研制的一种子午仪卫星定位系统(Transit)，1958然海陆空三军及民用部门都感到迫切需要一种新的卫星导航系统[13]。PositioningSystem”的缩写词。其意为“卫星测时测距导航/全球定位系维位置和速度信息,并且还广播一种形式的世界协调时(UTC)。通过遍布全球的(21+3)GPS导航卫星，向全球围的用户全天候提供高精度的导航、跟踪定位1]。美国只向外国提供低精度的卫星信号。.8颗卫星才能确保覆盖俄罗斯全境;如要提供全球定位服务，则需要24颗卫星。斗”的指挥机和终端之间可以双向交流。震发生后,武警指挥中心和武警部队组网高峰1.2本课题研究的意义和方法的能力，训练仿真模拟的技巧和方法。理和实现方法！1.3GPS前景坐标、操作简便、全天候作业、功能多、应用广泛等特点著称。.球性的地球动态参数；用于建立以及海洋测绘；用于检测地球板实现仅有少量的地面控制或无地面控制的航测快速成图，导致地理信息系统、全球环境遥感监测的技术革命[4]。的空间资源环境，促使国际民间形成了一个共为单一导航手段的最高应用境界。国际民间的际环境。.第二章GPS测量原理GPS导航系统的基本原理是测量出已知位置的卫星到用户接收机之间的距离,然后综合多颗卫星的数据就可知道接收机的具体位置。要达到这一目的,卫星的位置可以根据星载时钟所记录的时间在卫星星历中查出。而用户到卫星的距离则通过纪录卫星信号传播到用户所经历的时间,再将其乘以光速得到由(于大气层电离层的干扰,这一距离并不是用户与卫星之间的真实距离,而是伪距 (PR):当GPS卫星正常工作时,会不断地用1和0二进制码元组成的伪随机码(简称伪码)发射导航电文。GPS系统使用的伪码一共有两种,分别是民用的C/A码和军用的P(Y码)。C/A码频率1.023MHz重,复周期一毫秒码,间距1微秒,相当于300m;P码频率10.23MHz,重复周期266.4天码,间距0.1微秒,相当于30m。而Y码是在P码的基础上形成的,性能更佳。受机使用的时钟与卫星星载时钟不可能总是同步所,以除了用户的三维坐标x、y、z外,还要引进一个Δt即卫星与接收机之间的时间差作为未知数,然后用4个方程将这4个未知数解出来。所以如果想知道接收机所处的位置,至少要能接收到4个卫星的信号。2.1伪距测量的原理时刻是由接收机钟确定，这就在测定的卫星至接收机的距离中，不可避免地包含.之间的实际距离，所以称之为伪距。星位置和伪距观测值，采用距离交会法求出接收机的三维坐标和时钟改正数。是测量定位精度低，但足以满足部分用户的需要。2.1.1计算卫星位置读入导航电文后首先根据需要调用广播轨道1至广播轨道5上的数据然后依次计算卫星的平均角速度归化时间平均近点角需要注意的是进行真近点角计算时要同时计算正弦和余弦以得到正确象限的角计算经校正的升交点精度时需要用到地球旋转速率在WGS-84中这一常数为[10]:2.1.2用户位置的计算首先利用近似的用户位置与伪距观测值计算出一个近似伪距利用该近似伪距可以计算出部分值然后计算出系数并生成一个Nx4的矩阵为参与运算的卫星数最后算出用户位移的坐标上述过程根据需要可以将计算出的用户坐标作为近似值反复迭代直至符合精度要求[10]。.介绍的最小二乘法原理带到上面的公式当中去计算。2.1.3最小二乘法介绍动时间的安排使得用户在地球的任意位置(两极个别地点除外)，都能够看到的利用了已知的数据信息，使得结果的偏差最小化。例如：对于下面的方程：(2-2)如果令xy(2-3)使用最小二乘法，用C,H表示X.令(2-4)?p?p(2-5)xy.xy(2-6)整理得到 (2-7)(2-8)写成矩阵的形式就是：哪么就能够得到 ==a1 x[HTH]1HTC==(2-9)==y2.2载波相位测量原理载波相位观测方法是GPS接收机用接收到的卫星载波 (L1:154f0,19.032cm;L2:120f0,24.42cm)与本地接收机产生的本振参考载波产相位差来计算的.(GPS所接收到得载波相位是不连续的，所以在进行相位测.过码相关等方法重新获取载波)。以Q(tk)表示k接收机在时刻tk所接收到的第j颗卫星接收到的载波相位的值；以kk表示k接收机在时刻tk本地载波信号的相位值，则接收机在接收机钟面时刻tk时观测j卫星所取得的相位观测量可写为[7]：kkkkkk(2-10)SSj(t0)QInt(φ)一初始取样时刻(t0)kSj(ti)0图2-1载波相位测量原理图如上图1所示，在初始时刻t0,测得小于一周的相位差为Q，其整周数为.N0j，此时包含整周数的相位观测值为：(2-11)(2-11)整波计数器记录从t0到ti时间的整周数变化量INT(Q),只要卫星j从to到tj时间信号没有中断，则整周模糊数N0就j为一个常数，任意时刻ti卫星到k接这样，观测量就包含了相位差的小数部分和累计的整数部分的整周数。假设在GPS系统时刻Ta(卫星a时刻)卫星Sj发射的载波信号相位为收机)。接收机产生的本地载波相位为，根据(2-10)得到：在Tb时刻，收机)。接收机产生的本地载波相位为，根据(2-10)得到：在Tb时刻，bba(2-13)66ab则有：则有：bbaa(2-14)由于卫星和接收机的频率都比较稳定，所以在一个小的时间间隔里面，我们可以近似的理解在时间[t,t+编t]，有.(2-15) (上式是考虑在频率前提下，所以没有在f前面乘以2冗)(2-16)由2-15和2-16代入到2-14得到：bbaaba(2(2-18)f(2-19)这样，上式即为接收机k对卫星k的载波相位的以米为单位的观测方程式。.第三章GPS的坐标、时间系统测量站的位置。下面就天球坐标系和地球坐标系做简要的说明。3.1.1天球坐标系中，我们会涉及到几个参考点，线，面[9]。做天极P,与地球北极相对应的是北天极Pn,与地球南极相对应的是南天极Ps.天影像的天极叫做平天极。球赤道。天球赤道平面的夹角叫做黄赤交角，为23.5度。黄北极和黄南极，分别用Kn和Ks表示。和天球赤道面是建立天球坐标系的基准点和基准面。图3-1春分点的一种顺时针圆周运动。章动指的是真北天极绕平北天极做得椭圆型运动。图3-2岁差和章动天球坐标系分为两种：真天球坐标系和平天球坐标系.点,Y轴垂直于XMZ平面。.点,Y轴垂直于XMZ平面。句话说就是要考虑到岁差旋转和章动旋转地影响。3.1.2地球坐标系地球坐标系也可以分为两种：即平地球坐标系和真地球坐标系。年测定的平均纬度所确定的平均地极位置)。X轴：指向格林威治起始子午面与地球平赤道的交点。真地球坐标系Z轴：指向地球瞬时极。X轴：指向格林威治起始子午面与地球瞬时真赤道的交点。图3-3世界地心坐标系瞬时真天球坐标到瞬时真地球坐标的转换为对应的平格林威治起始子午面的真春分点时角Ω。故仅仅需要绕Z轴旋转这个角度Ω，就能够做到二者的相互转换。相应的转换转动矩阵为：Rz(Ω)=3-1-SinΩSinΩ00在0CosΩ0面监控站点，他们之间的转换问题一般可以按照下面的步骤来分析：.真天球坐标3.2时间系统时间系统是卫星定位测量过程中的一个重要概念。现时的GPS测量的方法PS间的距离，然后通过一定的数学方法以确定接收机所在的具体位置，为得到接的精确时间，因而利用卫星技术进行精密的定位和导航，必须要获得高精度的时间信息，这需要一个精确的时间系统。现行的卫星定位测量中与之紧密相关的时间系统有三种:世界时，原子时和动力学时。3.2.1世界时系统时，太阳时，世界时三种。1.恒星时选定春分点(地球赤道平面与其绕太阳公转轨道的一个交点)作为参考点，由点连续两次经过本地子午圈的时间间隔为一恒星日，其1/24为一恒星时，由于其定义涉及到地方子午圈，因而恒星时具有地方性，又称地方恒星时。当从格林尼治子午线上观测时，所得的恒星时称为格林尼治恒星时。由于地球自转受岁差、章动的影响，春分点的空间位置并不唯一，有真春分点和平春分点之分，这导致恒星时可分为真恒星时和平恒星时，因而对格林尼治恒星时有格林尼治 (3-2)其中，△笋为黄经章动，‘为黄经交角。2.太阳时于真太阳周年运动的平均速度，这个假定的参考点，在天文学上被称为平太阳。常称地方平太阳时。3.世界时universalTimeuT准完全一致，仅仅是起算点有所不同。若有编表示平太阳相对格林尼治子午圈的时角，定义有世界时UTO可表示为::(3-3)由于地球自转的不均匀性，使地球自转轴产生了极移现象因而UTO并不均自转变化的季节性改正参数双由此可得世界时UTI和uTZ:((3-4)(3-5)其中观测瞬时地极相对国际协议地极原点(CIO的)极移修正△兄的表达式为((3-6);凡，汽分别为天文经度和纬度。地;凡，汽分别为天文经度和纬度。地度的季节性变化改正△兀有如下的经验公式:(3-7)t为自本年起始日起算的年小数部分即(为计算时年积日与该年全年积日的比例)。上述修正并不能完全消除地球自转速度变化率和地球自转季节3.2.2原子时系统与之又有微小误差，关系为:(3-8)分析，利用数据处理推算出了统一的原子时系统——国际原子时作为其高精度的时间基准。3.2.3动力学时系统.的严格关系为:3.2.4协调世界时渐慢的影响，逐渐比原子时慢，为避免两都之间误差的扩大，自1972年起，国际上念，当协调时与世界时的时刻相差超过士0.9(s时)，便于协调时中引入闰秒士l(s)，闰秒一般于12月31日或6月30日加入。协调时与TAI的关系如下:3.2.5GPS时间系统据协调时与TAI的关系可得:(3-11)(3-12).第四章卫星运动基本定律及其求解加以描述。在各种力引力，称为在称为无摄运动。卫星在地球的引力场当中所做的无摄运动，也称为开普勒运动，4.1开普勒第一定律质心重合。这一个定律表明了，在中心引力的作用下，卫星绕地球轨道运行的轨道面，勒椭圆，其形状和大点。轨道图形可以表示为如下图5：远地点远地点P’sM勒椭圆的偏心率；fs为真近点角，它描述了任意时刻，卫星在轨道上面相对于近地点的位置，是时间的函数，其定义见上图所示。律阐述了卫星运动轨道的基本形态及其与地心的关系。4.2开普勒第二定律地球质心与卫星质心间的距离向量，在相同的时间所扫过的空间面积是相等的。(如下图6所示)图4-2卫星地心向径在相同的时间间隔内扫描的面积 (图2).与任何其它的运动物体一样，在轨道上面运动的卫星，也具有两种的能量：的运动的速度为vs,则其动能为1mv2。根据能量守恒定律，卫星的势能与动2ss能的总量是不变的，即1mv2+GMms2ssr=常量(4-2)因此，当卫星运行到近地点的时候，其动能最大；在远地点的时候，其动能断变化的，在近地点处的速度最大，而在远地点的速度最小。三定律为：T2T42G(4-3)在这个式子当中，Ts为卫星的运行周期.如果我们假设卫星的平均角速度为N,则有:.于是，开普勒第三定律4-2就可以写成：As(4-4)(4-或者表示为常用的形式：N=5)GGAs3(4-6)且保持不变。参数卫星的无摄运动，一般的可以由下面的6个参数(图7)来描述：As-----------------------------卫星轨道的长半径Es----------------------------卫星轨道的偏心率Ω-------------------升交点的赤径i--------------------卫星轨道面的倾角ωs--------------近地点角距，即升交点与近地点的夹角fs-----------------卫星的真近点角，在轨道平面上为卫星与进地点的地心角距。图4-3开普勒轨道参数度，就被唯一的确定了！其求解在描述卫星无摄运动的6个参数当中，只有fs是关于时间的函数，其他的都是一般的参数。所以，计算卫星瞬时的位置的关键，计算出参数fs，并由此确定卫星的空间位置及其和时间的关系。的圆心角就是Es。近点角定义为：.Msntt(4-7)于一个确定的卫星来说，这个参数可以认为是常数。的关系成立：于是将上式带入到(4-1)中就得到：Cosfs=se1cosEss或者得到以下常用的形式：求解(4-10)(4-11)间的位置便可以随之确定。若以直角坐标的原点与地心M重合，轴指向近地点且垂直于轨道的ss.sGcosfs0 ns=sinfs(4-12).=ss)s0(4-13将坐标系(,n,G)依次做下面的变化旋转：sss (1)绕轴顺时针旋转角度ωs，使轴的指向由近地点变为升交点。 (2)绕轴顺时针旋转角度i，使轴与Z轴相同。ss (3)绕轴顺时针旋转角度Ω，使轴指向春分点。ss(4-14)(4-15)010-SinωsSinωs000Sini-SinΩ0100-Sini.(4-16)0SinΩ00天球坐标系下面的卫星位置坐标可以表示为：nns0(4-17)s1利用转换关系Rz(Θg)得到相应的地球坐标系的坐标了！CosΘgSinΘg0而Rz(Θg)=-SinΘg(4-18)0.第五章GPS的matlab仿真。一颗卫星信号能否被接收与下列因素是有关系的[7]：的影响。如果卫星处于图中的阴影部分，则对图中的飞机是不能够接收到该卫星的信号的!EEh水平面的角13.9度，这样就保证了一些飞行高度较高的航天器在高空可以更多.的接收到GPS卫星的信号。但是对于那些超出发射角的飞行器就收不到信号了。如果这个角度小于90度的话，就可以收到，反之不能够收到。特别的，当这个角度刚刚为90度的时候，我们一般认为是收不到的。程序主体见附录在进行仿真之前，有几个子程序段需要说明一下：functionplot3c(x,y,z,color)switch(color)case0plot3(x,y,z,'w--');plot3(x,y,z,'r--');case2plot3(x,y,z,'g--');case3plot3(x,y,z,'c--');case4plot3(x,y,z,'m--');case5plot3(x,y,z,'y--');.case6plot3(x,y,z,'b--');case7plot3(x,y,z,'k--');次代表的是红色,绿色,青绿色品.红色,黄色,蓝色,黑色.functionboxplot3(x,y,z,Lx,Ly,Lz,color)index=zeros(6,5);index(1,:)=[12341];index(2,:)=[56785];index(3,:)=[12651];index(4,:)=[43784];index(5,:)=[26732];index(6,:)=[15841];fork=1:6plot3c(x(index(k,:)),y(index(k,:)),z(index(k,:)),color).主要是用来表示用户的空间位置的。functiondrawearth(time)j1=[0:pi/10:2*pi];w1=[-pi/2:pi/10:pi/2];L2=length(j1);forn=1:L1z=ones(L2,1);z=z*r*sin(w1(n));x=temp*sin(j1);plot3(x,y,z);holdon;grid;.%figure(3);unit=ones(1,1);z0=ones(1,1);x0=ones(1,1);y0=ones(1,1);forn=1:L2z=r*sin(temp);z1=unit*z;x1=unit*x;y1=unit*y;z1];x1];y1];.z0(:,1)=[];x0(:,1)=[];y0(:,1)=[];plot3(x0,y0,z0);axisequal;axisoff;holdon;boxplot3(0,0,0,100,100,100,7);%标示出来地球球心的位置。方便观察.它的作用就是画出一个在上述地球坐标图当中的卫星。仿真程序一：绘制卫星的轨道平面程序见附录：地球地球质心图5-2卫星轨道平面仿真图仿真程序二：单颗卫星不同时刻的动态仿真程序见附录：图5-3单颗卫星动态仿真图由上图10可知，卫星的天空瞬时位置是随着时间的变动而发生变化的。对星是不相同的。这个仿真程序的功能实际上就是仿真了在一个轨道上面的卫星在不同的时候 (这里以一个小时为一个观测时元，进行动态的在屏幕上显示器位置)通过改变的相对位置的变化)。哪么，在其他轨道面上的卫星的运动也可以类似的模拟出仿真程序三：卫星在某个时刻的全轨道平面的分布和可见卫星程序见附录：图5-4全轨道平面的图形别如下：SatellitePosition=.011这个程序仿真了在某个时刻(在程序里面是在时刻timenow=0)的全部24可以在同一个时刻看到不同(4到11颗)的可见卫星。对于在这个程序当中的可见卫星的及时在轨坐标如下：用用户位置地心位置图5-5用户可见的卫星分布SatellitePosition=.2.28831.07591.33241.21270.44431.4000-0.8836-2.5255-1.8178-0.9774-0.9295-1.3073-2.0297-0.6021-1.5812-1.5756-1.4132-0.0000-2.1144-2.1757-1.3835-1.9058-0.7516-1.8463-0.7869-2.1091-1.3835-1.9217-2.1704-0.736500000000000000仿真程序四：用可见卫星计算用户位置程序见附录法逼近计算。对于前面的假设用户位置(6400，图5-6用户位置的计算由图13可见，根据最小二乘法的原理，计算出来了用户在一个时刻的位置标，并把它表示在了这个天球坐标系当中。calculaterecord=000根据程序当中的参数Error的设定而有所出入。图5-7用户位置的计算是很直观的反应出我们基于最小二乘法的伪距算法的主要思路。如下图所示：图5-8用户位置的计算用户的位置，第七箭头表示的用户的实际位置(用白色的柱体表示)。.结论所以本次设计难度很大，但在最后都还是通过查找资料或者老师同学帮助也都一一克服。间不断变化的，因此，在给出卫星运行位置的同时，必须给出相应的瞬时时刻。所以时间系统是GPS卫星计算里面很重要的一环。要计算卫星的运动轨迹及位素影响。通常把作用于卫星上的各种力，按其影响的大小分为两类：中心力和摄球引力场中无摄运动，也称为开普勒运动。所以本次设计中主要用开普勒定律计算其位置轨迹，在运算中还加入了最小二乘法计算参数数据。此次设计中卫星的幅角的计算。所以，根据卫星轨道的星座参数，找出在某时刻的近地点幅角(一般是零时刻)，然后在此基础上以时间为基础进行叠加，确定出每个时间点的瞬时近地点幅角值。程序从自定义的星座数据文件中建立起卫星星座。设计一个子函数，命名为calculateuserposition.m。其功能为其他程序得出一个全局变量.个模块：三维旋转图标，可以很清晰的从不同角度看到卫星和坐标原点；全部24颗卫星的轨道图形，以及对于用户来说在该时刻可见卫星位置；的位置坐标。也有了更直观的认识，对于以后遇到GPS的应用时有了一个简单的基础。.致同学表示万分的感！完成本次设计的青老师。在本次设计中，我遇到了许多都深深感染了我，将会让我受益终身。在此，我向老师表示深深的感，并致以崇高的敬意！华老师、周伟老师。大学四年期间，在大一大二无知的高中生脱变到具有独立思考能力的大学生，并提供给我许多锻炼自身能力的机会。在大三大四的时候周伟老师出任我们班主任，他作为一个专业课老师，要感理工大学核技术与自动化工程学院的领导们四年来给予的多方面的关心和帮助!感所有的朋友和同学，让我度过了一个快乐而又充实的大学时光！完成了学业。.[1]周忠谟,易杰军.GPS卫星测量原理与应用[M].:测绘,1997.[2]征航,黄劲松.GPS测量与数据处理[M].:大学,2005.[3]ElliottD.Kaplan.UnderstandingGPS:principlesandapplications[M]邱.致和,王万义,译.电:子工业,2002.[4]马中元.GPS卫星全球定位系统简介[J].科技经济市场2002(8)[5]郭秋英胡,振琪.GPS卫星坐标的计算[J].全球定位系统,2006(3).[6]夏振纯,金志华田,蔚风.基于MATLAB的GPS选星分析及实验研究[J].计算机仿真,2005,22(1):51-54.[7]俊，武奇生编著.GPA基本原理及其Matlab仿真[M].：电子科技大学[9]汤长存.GPS卫星定位仿真分析与系统程序设计[D].同济大学工学硕士学位论文LABGPSC中南分网第二十四次学术信息交流会论文集[11]毛敏.GPS卫星定位系统简介[Z]国.防工业职业技术学院学报第二十卷第一期。[12郭]秋英，胡振琪.GPS卫星坐标的计算[J].2006[13]朱明.全球卫星定位系统GPS及在现代交通中的应用[J].[14]楼益栋，施闯，茂荣，齐乐.GPS卫星实时精密定轨及初步结果分析[J].2008.8[15]永法.运用GPS卫星定位系统，提高车辆利用率[J].2008第.二十七期[16]齐乐.GPS导航星座及低轨卫星的精密定轨理论和软件研究[D].:大学,2004[17]ChenK.Real2TimePrecisePointPositioningandItsPotentialApplications[C].IONGNSS17thIn2ternationalTechnicalMeetingoftheSatelliteDivi2sion,LongBeach,CA,2004[18]LiuJingnan,GeMaorong.PANDASoftwareandItsPreliminaryResultofPositioningandOrbitDet2ermination[J].WuhanUniversityJournalofNaturalSciences,2003,8(B2):603-609[19]KoubaJ,HerouxP.GPSPrecisePointPositioningUsingIGSOrbitProjucts[J].GPSSolution,2001,5(2):12-28[20]..附录程序一(绘制卫星的轨道平面)代码functiondrawsatelliteorbitE=[0:0.1:2*pi];x=a*(cos(E)-e);y=a*sqrt((1-e^2))*sin(E);A1=[32.892.8152.8212.6272.8332.8];%卫星星座数据.fork=1:6sin(A)000-sin(A)0;cos(A)0;01;];0;cos(B)-sin(B);sin(B)cos(B);];.R3=[cos(C)-sin(C)0;sin(C)cos(C)0;001;];L1=length(E);Ans=R312*[x;y;z;];x1=Ans(1,:);y1=Ans(2,:);z1=Ans(3,:);plot3c(x1,y1,z1,k);boxplot3(0,0,0,200,200,200,7);holdon;axisequal;axisoff;程序二(单颗卫星不同时刻的动态仿真)代码clear;clc;closeall;j=0;.E=[0:0.1:2*pi];x=a*(cos(E)-e);y=a*sqrt((1-e^2))*sin(E);drawearth(0);holdon;axison;sin(A)000-sin(A)0;cos(A)0;01;];0;cos(B)-sin(B);sin(B)cos(B);];sin(C)0-sin(C)cos(C)00;0;1;];L1=length(E);.Ans=R312*[x;y;z;];x1=Ans(1,:);y1=Ans(2,:);z1=Ans(3,:);plot3c(x1,y1,z1,2);holdon;axisequal;axison;gridon;x=a*(cos(C)-e);y=a*sqrt((1-e^2))*sin(C);sin(A)00-sin(A)cosA)00cos(B)0;0;1;];0;-sin(B);.0sin(C)0sin(B)cos(B);];-sin(C)0;cos(C)0;01;];L1=length(E);Ans=R312*[x;y;z;];x1=Ans(1,:);y1=Ans(2,:);z1=Ans(3,:);drawsatellite(x1,y1,z1,6);axison;程序三(卫星在某个时刻的全轨道平面的分布和可见卫星)代码clear;clc;closeall;E=[0:0.1:2*pi];.x=a*(cos(E)-e);y=a*sqrt((1-e^2))*sin(E);globalSatellitePositionglobalttumSatellitePosition=ones(24,4);figure(1);drawearth(0);holdon;A1=[32.892.8152.8212.6272.8332.8];drawsatelliteorbit;fork=1:6.x=a*(cos(C)-e);y=a*sqrt((1-e^2))*sin(C);sin(A)000-sin(A)0;cos(A)0;01;];00;cos(B)-sin(B);sin(B)cos(B);];sin(C)0-sin(C)cos(C)00;0;1;];L1=length(E);Ans=R312*[x;y;z;];x1=Ans(1,:);y1=Ans(2,:);z1=Ans(3,:);.drawsatellite(x1,y1,z1,k);SatellitePosition(temp,:)=[x1y1z11];ttum(temp,:)=[x1y1z11];holdon;于90°,就认为改卫星是不可见的.earthcenterpos=[000];置fork=1:24temp=SatellitePosition(k,1:3)-userposition;temp=SatellitePosition(k,1:3)-earthcenterpos;jiajiao=acos((Dist2+Dist3-Dist1)/2/sqrt(Dist3)/sqrt(Dist2));if(jiajiao>=pi/2)SatellitePosition(k,4)=0;.%ttum(k,:)=SatellitePosition(k,1:3);figure(2)drawearth(0);holdon;drawsatelliteorbit;j=1;fork=1:6if(SatellitePosition(j,4)==1)tempx=SatellitePosition(j,1);tempy=SatellitePosition(j,2);tempz=SatellitePosition(j,3);dra