www.gusucode.com > gps卫星位置预报，MATLAB编写，有动画界面源码程序 > code5/gps_final_new/TestAlmFromGPS.m

    clear;clc;close all;

PRN=zeros(32,1);
type=zeros(32,1);
Health=zeros(32,1);
es=zeros(32,1);
Toe=zeros(32,1);
delta_i=zeros(32,1);
dot_omega=zeros(32,1);
sqrt_as=zeros(32,1);
omega0=zeros(32,1);
w=zeros(32,1);
M0=zeros(32,1);
af0=zeros(32,1);
af1=zeros(32,1);
week=zeros(32,1);

SatPosX=zeros(32,1);
SatPosY=zeros(32,1);
SatPosZ=zeros(32,1);

UserPosLat=32.02*pi/180;
UserPosLon=118.48*pi/180;
UserPosHeight=20;
[UserPosX,UserPosY,UserPosZ]=llh2xyz(UserPosLat,UserPosLon,UserPosHeight);
UserPos=[UserPosX/1000.0 UserPosY/1000.0 UserPosZ/1000.0];

EarthCenPos=[0 0 0];

opentype=0;  % 0 from gps realtime almanac data; 1 from internet almanac data
if opentype==0
    fid=fopen('gps.log','r');
    BinData=fread(fid,85*32);
    fclose(fid);
    for i=0:31
        PRN(i+1)=bitand(BinData(5+i*85),63);
        type(i+1)=bitshift(bitand(BinData(5+i*85),192),-6);
        Toe(i+1)=BinData(6+i*85)+BinData(7+i*85)*256+BinData(8+i*85)*256*256+BinData(9+i*85)*256*256*256;
        af0(i+1)=ProcessDouble(BinData(10+i*85:17+i*85));
        af1(i+1)=ProcessDouble(BinData(18+i*85:25+i*85));
        M0(i+1)=ProcessDouble(BinData(26+i*85:33+i*85));
        w(i+1)=ProcessDouble(BinData(34+i*85:41+i*85));
        omega0(i+1)=ProcessDouble(BinData(42+i*85:49+i*85));
        sqrt_as(i+1)=ProcessDouble(BinData(50+i*85:57+i*85));
        dot_omega(i+1)=ProcessDouble(BinData(58+i*85:65+i*85));
        delta_i(i+1)=ProcessDouble(BinData(66+i*85:73+i*85));
        es(i+1)=ProcessDouble(BinData(74+i*85:81+i*85));
        week(i+1)=BinData(82+i*85)+BinData(83+i*85)*256;
        %[PRN(i),type(i),Toe(i),af0(i),af1(i),M0(i),w(i),omega0(i),sqrt_as(i),dot_omega(i),delta_i(i),es(i),week(i)]=ProcessGpsData(BinData((i-1)*85+1:i*85));
    end
else
    cnt=0;
    fid=fopen('almanac_yuma_week0552.txt','r');
    while 1
        tline=fgetl(fid);
        cnt=cnt+1;
        if ~ischar(tline),   break,   end
    end
    fclose(fid);

    recordnum=floor(cnt/15);

    fid=fopen('almanac_yuma_week0552.txt','r');
    for m=1:recordnum
        tline=fgetl(fid);
        tline=fgetl(fid);
        [temp,tempx]=strread(tline,'%s %d');
        i=tempx;
        PRN(i)=i;
        tline=fgetl(fid);
        [temp,Health(i)]=strread(tline,'%s %d');
        tline=fgetl(fid);
        [temp,es(i)]=strread(tline,'%s %f');
        tline=fgetl(fid);
        [temp,temp1,temp2,Toe(i)]=strread(tline,'%s %s %s %f');
        tline=fgetl(fid);
        [temp,temp1,delta_i(i)]=strread(tline,'%s %s %f');
        tline=fgetl(fid);
        [temp,temp1,temp2,temp3,dot_omega(i)]=strread(tline,'%s %s %s %s %f');
        tline=fgetl(fid);
        [temp,temp1,temp2,sqrt_as(i)]=strread(tline,'%s %s %s %f');
        tline=fgetl(fid);
        [temp,temp1,temp2,temp3,omega0(i)]=strread(tline,'%s %s %s %s %f');
        tline=fgetl(fid);
        [temp,temp1,temp2,w(i)]=strread(tline,'%s %s %s %f');
        tline=fgetl(fid);
        [temp,temp1,M0(i)]=strread(tline,'%s %s %f');
        tline=fgetl(fid);
        [temp,af0(i)]=strread(tline,'%s %f');
        tline=fgetl(fid);
        [temp,af1(i)]=strread(tline,'%s %f');
        tline=fgetl(fid);
        [temp,week(i)]=strread(tline,'%s %d');
        tline=fgetl(fid);
    end 
    fclose(fid);
end

for t=0:11
    
    figure(t+1)
    set (gcf,'Position',[0,30,1200,800], 'color','w')
    DrawEarth(0);%地球可视化
    DrawUser(UserPos);%用户可视化
    axis([-30000 30000 -30000 30000 -30000 30000]);
    hold on;
    
    T=Toe(1,1)+t*3600*2;
    
    for i=1:32
        %从历书计算卫星位置
        [SatPosX(i),SatPosY(i),SatPosZ(i)]=CompSatPosFromAlm(PRN(i),T,Toe(i),M0(i),sqrt_as(i),es(i),w(i),delta_i(i),omega0(i),dot_omega(i));        
        SatPos(i,1)=SatPosX(i);
        SatPos(i,2)=SatPosY(i);
        SatPos(i,3)=SatPosZ(i);
        SatPos(i,4)=1;
        
        %可见性计算
        Temp=SatPos(i,1:3)-UserPos;
        Dist1=Temp*Temp';
        Temp=UserPos-EarthCenPos;
        Dist2=Temp*Temp';
        Temp=SatPos(i,1:3)-EarthCenPos;
        Dist3=Temp*Temp';
        ViewAngle=acos((Dist1+Dist2-Dist3)/2/sqrt(Dist1)/sqrt(Dist2));
        if(ViewAngle<=pi/2)
            SatPos(i,4)=0;
        end
        
        %卫星可视化
        DrawSat(SatPos(i,:));
        DrawSatOrbit(PRN(i),T,Toe(i),M0(i),sqrt_as(i),es(i),w(i),delta_i(i),omega0(i),dot_omega(i));
    end

    drawnow;
    M(t+1)=getframe;
end

close all;
figure(111);
set (gcf,'Position',[0,30,1200,800], 'color','w')
axis equal;
axis off;
%动画
movie(M,12,2);