Demo entry 6736309

cal

   

Submitted by anonymous on Apr 26, 2018 at 14:08
Language: Matlab. Code size: 4.1 kB.

RINEXN='F:\大三下\GNSS\作业&&实习\实习1:计算卫星轨道\数据\55191102.18N';  
cc=20;  
yy=18;  
mm=4;  
dd=20;  
hh=4;  
min=00;  
ss=00;    
function position=CalSatPos(RINEXN,cc,yy,mm,dd,hh,min,ss)  
%RINEXN为N文件地址  
%cc是公元年份的前两位数  
%yy是已知公元年份的后两位数?  
%mm是月数  
%dd是日数  
clear  
clc  
%原理:蔡勒公式计算星期几,0代表星期日(西方星期日为第一天)  
W=rem((floor(cc/4)-2*cc+yy+floor(yy/4)+floor(26*(mm+1)/10)+dd-1),7);  
gpstoc=(W*24+hh)*3600+min*60+ss;  
  
%跳过文件头读取数据,每行数据读取结果为一个字符串  
data=textread(RINEXN,'%s','delimiter','\t','headerlines',7);  
ls=size(data,1);  
%N代表该星历文件中涵盖的卫星数  
N=ls/8;  
%data1记录星历文件数据  
data1=zeros(8,4,N);  
Result=zeros(N,4);  
  
%通过正则化将数据存入数组  
for i=1:ls  
    %判断第几颗卫星  
    id=ceil(i/8);  
    r=rem(i,8);  
    %判断是否为每颗卫星的第一行数据  
    if(r==1)  
        %将字符串按照空格分割,并去除空白数  
        S = regexp(data(i), '\s', 'split');  
        xs=S{1};  
        tf=not(cellfun(@isempty,xs));  
        new_cell=xs(tf);  
        %存储卫星号  
        xa=new_cell(1);  
        xb=strtok(xa,'G');  
        data1(1,1,id)=str2num(xb{1});  
        Result(id,1)=str2num(xb{1});  
        %存储每颗卫星前三个参数  
        xc=new_cell(8);  
        xd=new_cell(9);  
        xe=new_cell(10);  
        data1(1,2,id)=str2num(xc{1});  
        data1(1,3,id)=str2num(xd{1});  
        data1(1,4,id)=str2num(xe{1});  
    else  
        for j=1:4  
            S=regexp(data(i),' ', 'split');  
            xs=S{1};  
            tf=not(cellfun(@isempty,xs));  
            new_cell=xs(tf);  
            xx=new_cell(j);  
            if(r==0)  
                r=8;  
            end  
            data1(r,j,id)=str2num(xx{1});  
        end  
    end  
end  
  
%根据数据计算卫星位置  
for k=1:N  
    a0=data1(1,2,k);    %卫星时钟偏差  
    a1=data1(1,3,k);    %卫星时钟漂移  
    a2=data1(1,4,k);    %卫星时钟漂移  
    at=data1(2,1,k);    %星历数据的有效龄  
    Crs=data1(2,2,k);   %轨道半径正弦改正项  
    dn=data1(2,3,k);    %平均运动修正项  
    M0=data1(2,4,k);    %t0e时的平近点角  
    Cuc=data1(3,1,k);   %纬度幅角余弦改正项  
    e=data1(3,2,k);     %卫星轨道偏心率  
    Cus=data1(3,3,k);   %纬度幅角正弦改正项  
    sqA=data1(3,4,k);   %轨道长半径平方根  
    t0e=data1(4,1,k);   %星历的基准时  
    Cic=data1(4,2,k);   %轨道倾角余弦调和项  
    Om=data1(4,3,k);    %升交点赤经  
    Cis=data1(4,4,k);   %轨道倾角正弦项  
    i0=data1(5,1,k);    %轨道倾角  
    Crc=data1(5,2,k);   %轨道半径余弦调和项  
    w=data1(5,3,k);     %近地点角  
    Omdot=data1(5,4,k); %升交点赤经变化率  
    Idot=data1(6,1,k);  %轨道倾角变化率  
    L2=data1(6,2,k);    %L2上的码  
    weeks=data1(6,3,k); %星期数(TOE)  
    L2p=data1(6,4,k);   %p码数据标识  
    p=data1(7,1,k);     %卫星精度  
    SH=data1(7,2,k);    %卫星健康(MSB  
    TGD=data1(7,3,k);   %单频机延迟改正数  
    IODC=data1(7,4,k);  %时钟数据有效期  
    st=data1(8,1,k);    %电文发送时间  
    %计算卫星钟差改正,对观测时刻t归化到GPS时系  
    dt=a0+a1*(gpstoc-t0e)+a2*(gpstoc-t0e)*(gpstoc-t0e);  
    t=gpstoc+dt;  
    tk=t-t0e;  
    if(tk>302400)  
        tk=tk-604800;  
    elseif(tk<-302400)  
        tk=tk+604800  
    end;  
    tk=a0+a1*tk+a2*tk*tk;  
    %计算卫星运动的平均角速度n  
    GM=3.986005e14;  
    n0=sqrt(GM)/((sqA)^3);  
    n=n0+dn;  
    %计算观测瞬间卫星的平近点角M  
    M=M0+n*tk;                
    %计算偏近点角,迭代法解方程  
    oldE=M;  
    E=M+e*sin(oldE);  
    de=1e-12;  
    while(abs(E-oldE)>de)  
        oldE=E;  
        E=M+e*sin(oldE);  
    end  
    %计算真近点角f  
    f=atan((sqrt(1-e^2)*sin(E))/(cos(E)-e));  
    %计算升交角距uu  
    uu=w+f;  
    %计算摄动改正项du,dr,di  
    du=Cuc*cos(2*uu)+Cus*sin(2*uu);  
    dr=Crc*cos(2*uu)+Crs*sin(2*uu);  
    di=Cic*cos(2*uu)+Cis*sin(2*uu);  
    %对uu,r',i0进行摄动改正  
    u=uu+du;  
    r=(sqA^2)*(1-e*cos(E))+dr;  
    i=i0+di+Idot*tk;  
    %计算卫星在轨道面坐标系中的位置  
    x=r*cos(u);  
    y=r*sin(u);  
    %计算观测瞬间升交点的经度L  
    we=7.292115*10^(-5);  
    L=Om+(Omdot-we)*t-Omdot*t0e;  
    %计算卫星在瞬时地球坐标系中的位置  
    x1=x*cos(L)-y*cos(i)*sin(L);  
    y1=x*sin(L)+y*cos(i)*cos(L);  
    z1=y*sin(i);  
    %计算卫星在协议地球坐标系中的位置  
    %[xcts ycts zcts]=[1,0,xp;0,1,-yp;-xp,yp,1].*[X Y Z];  
    Result(k,2)=x1;  
    Result(k,3)=y1;  
    Result(k,4)=z1;  
end  
position=Result;  
end  

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).