keyongxing.asv

RAIM 完好性分析

clear all;
clc;
buffer_t_gps =[];
buffer_t_galileo =[];
buffer_t = [];
HPL_buffer = [];

gps_flag = 1;
galileo_flag = 2;
bd_flag = 3;

sigma = 33;
tao = 120;
P_fa = 0.002/(3600/tao);
ky = 0;
ky1 = 0;
js = 0;
for t=0:10:86400
    js = js+1;
    %--------GPS---------------------------------------
    [Xs_gps,Ys_gps,Zs_gps]=gps_sat(t); %解卫星在地心地固坐标系中的位置
%     %--------Galileo-----------------------------------
%     [Xs_galileo,Ys_galileo,Zs_galileo]=galileo_sat(t);
%     %--------Beidou-----------------------------------
%     [Xs_beidou,Ys_beidou,Zs_beidou]=bd_sat(t); %解卫星在地心地固坐标系中的位置
    %------用户轨迹------------------------------------     
    [long,lat,h]= usercontrail(t);
    [xu,yu,zu] = user_ecef(long,lat,h);
    [xu1,yu1,zu1]=ecef_g(long,lat,xu,yu,zu);             %转换成用户地点的地理坐标
    
    [Xs1_gps,Ys1_gps,Zs1_gps]=ecef_g(long,lat,Xs_gps,Ys_gps,Zs_gps);
%     [Xs1_galileo,Ys1_galileo,Zs1_galileo]=ecef_g(long,lat,Xs_galil                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                eo,Ys_galileo,Zs_galileo);
%     [Xs1_beidou,Ys1_beidou,Zs1_beidou]=ecef_g(long,lat,Xs_beidou,Ys_beidou,Zs_beidou);
        
    [num_gps,k_gps] = get_seesate(Xs1_gps,Ys1_gps,Zs1_gps,xu1,yu1,zu1,gps_flag); 
%     [num_galileo,k_galileo] = get_seesate(Xs1_galileo,Ys1_galileo,Zs1_galileo,xu1,yu1,zu1,galileo_flag); 
%     [num_beidou,k_beidou] = get_seesate(Xs1_beidou,Ys1_beidou,Zs1_beidou,xu1,yu1,zu1,bd_flag);  

    [X_gps,Y_gps,Z_gps]=get_k(num_gps,k_gps,Xs_gps,Ys_gps,Zs_gps,gps_flag);        %算出k颗可见星的坐标(地心地固坐标)
%     [X_galileo,Y_galileo,Z_galileo]=get_k(num_galileo,k_galileo,Xs_galileo,Ys_galileo,Zs_galileo,galileo_flag);        %算出k颗可见星的坐标(地心地固坐标)
%     [X_beidou,Y_beidou,Z_beidou] = get_k(num_beidou,k_beidou,Xs_beidou,Ys_beidou,Zs_beidou,bd_flag);
% 
    [XS_gps,YS_gps,ZS_gps] = ecef_g(long,lat,X_gps,Y_gps,Z_gps); % 转换为地理坐标
%     [XS_galileo,YS_galileo,ZS_galileo] = ecef_g(long,lat,X_galileo,Y_galileo,Z_galileo); % 转换为地理坐标
%     [XS_beidou,YS_beidou,ZS_beidou] = ecef_g(long,lat,X_beidou,Y_beidou,Z_beidou); % 转换为地理坐标     
    
    H_gps = Single_get_H(k_gps,XS_gps',YS_gps',ZS_gps',xu1,yu1,zu1);
    slopemax_gps = get_slopemax(H_gps,k_gps,1);
    freedom_gps = k_gps-4;
%     [T_gps,Td_gps] = kaifang_ceiling(sigma,freedom_gps,P_fa);
%     T_min = [233.58,172.40,144.97,128.50,117.20,108.81,102.27,96.97];
%     % sigma = 33
      T_min = [233.58,172.40,144.97,128.50,117.20,108.81,102.27,96.97];

    HPL_gps=slopemax_gps*T_min(freedom_gps);

    HPL_buffer = [HPL_buffer,HPL_gps];
    HAL_gps = 556;

    if HPL_gps<HAL_gps
        ky = ky+1;
    else
        buffer_t_gps = [buffer_t_gps,t];
        ky = ky;
    end
%     H_galileo=Single_get_H(k_galileo,XS_galileo',YS_galileo',ZS_galileo',xu1,yu1,zu1);
%     slopemax_galileo=get_slopemax(H_galileo,k_galileo,1);
%     freedom_galileo = k_galileo-4;
%     [T_galileo,Td_galileo] = kaifang_ceiling(sigma,freedom_galileo,P_fa);    
%     ARP_galileo=slopemax_galileo*Td_galileo;
    
%     k = k_gps+k_galileo+k_beidou;
%     H = Third_get_H(k_gps,k_galileo,k_beidou,XS_gps',YS_gps',ZS_gps',XS_galileo',YS_galileo',ZS_galileo',XS_beidou',YS_beidou',ZS_beidou',xu1,yu1,zu1);
%     slopemax=get_slopemax(H,k,3);
%     freedom_3 = k-6;
%     [T,Td] = kaifang_ceiling(sigma,freedom_3,P_fa);
%     ARP=slopemax*Td;


%       arp=[287,304,314,324,332,340,346,355];% ????????????? bo be continue............  
% %     arp = [331,345,355,365,372,379,385,391,397,402,407,411,415,420,423,428,432,435];
%     if ARP_gps<arp(1,freedom_gps)
%     else
%         buffer_t_gps=[buffer_t_gps,t];   
%     end

%     if ARP_galileo<arp(1,freedom_galileo)
%     else
%         buffer_t_galileo=[buffer_t_galileo,t];   
%     end
%     
%     if ARP<arp(1,freedom_3)
%     else
%         buffer_t=[buffer_t,t];   
%     end
  end
 
 rate = ky/js


% plot(biasbuffer,ratebuffer,'b-',biasbuffer,ratebuffer_s,'r-','LineWidth',2);
% legend('gps/galileo/beidou','gps');
% title('故障检测率');
% xlabel('偏差(m)');
% ylabel('故障检测率');
% grid on