clear all;
close all;
clc;
gps_flag = 1;
galileo_flag = 2;
bd_flag = 3;
buffer_dop_gps = [];
buffer_dop_galileo = [];
buffer_dop = [];
buffer_t = [];
time_bin = 0;
time_sample = 10;
time_end = 24*3600;
tt = 0;
for t = time_bin:time_sample:time_end
    buffer_t = [buffer_t,t];
    tt = tt+1;
     %------用戶軌跡------------------------------------  
    [long,lat,h]= usercontrail(t);
    [xu,yu,zu] = user_ecef(long,lat,h);
    [xu1,yu1,zu1]=ecef_g(long,lat,xu,yu,zu);             %轉換成用戶地點的地理坐標
    
    %--------GPS---------------------------------------
    [Xs_gps,Ys_gps,Zs_gps]=gps_sat(t); %解衛星在地心地固坐標系中的位置
    [Xs1_gps,Ys1_gps,Zs1_gps]=ecef_g(long,lat,Xs_gps,Ys_gps,Zs_gps);
    [num_gps,k_gps] = get_seesate(Xs1_gps,Ys1_gps,Zs1_gps,xu1,yu1,zu1,gps_flag); 
    [X_gps,Y_gps,Z_gps]=get_k(num_gps,k_gps,Xs_gps,Ys_gps,Zs_gps,gps_flag);        %算出k顆可見星的坐標(地心地固坐標)
    [XS_gps,YS_gps,ZS_gps] = ecef_g(long,lat,X_gps,Y_gps,Z_gps); % 轉換為地理坐標
    xsg1_gps = XS_gps';ysg1_gps = YS_gps';zsg1_gps = ZS_gps';
    H_gps = Single_get_H(k_gps,XS_gps',YS_gps',ZS_gps',xu1,yu1,zu1); 
    [dop_gps] = get_dop(H_gps);
    buffer_dop_gps = [buffer_dop_gps,dop_gps];
    %--------Galileo-----------------------------------
    [Xs_galileo,Ys_galileo,Zs_galileo]=galileo_sat(t);
    [Xs1_galileo,Ys1_galileo,Zs1_galileo]=ecef_g(long,lat,Xs_galileo,Ys_galileo,Zs_galileo);
    [num_galileo,k_galileo] = get_seesate(Xs1_galileo,Ys1_galileo,Zs1_galileo,xu1,yu1,zu1,galileo_flag); 
    [X_galileo,Y_galileo,Z_galileo]=get_k(num_galileo,k_galileo,Xs_galileo,Ys_galileo,Zs_galileo,galileo_flag);        %算出k顆可見星的坐標(地心地固坐標)
    [XS_galileo,YS_galileo,ZS_galileo] = ecef_g(long,lat,X_galileo,Y_galileo,Z_galileo); % 轉換為地理坐標
    xsg1_galileo = XS_galileo';ysg1_galileo = YS_galileo';zsg1_galileo = ZS_galileo';
    H_galileo = Single_get_H(k_galileo,XS_galileo',YS_galileo',ZS_galileo',xu1,yu1,zu1); 
    [dop_galileo] = get_dop(H_galileo);
    buffer_dop_galileo = [buffer_dop_galileo,dop_galileo];
    
    %--------Beidou-----------------------------------
    [Xs_beidou,Ys_beidou,Zs_beidou]=bd_sat(t); %解衛星在地心地固坐標系中的位置
    [Xs1_beidou,Ys1_beidou,Zs1_beidou]=ecef_g(long,lat,Xs_beidou,Ys_beidou,Zs_beidou);
    [num_beidou,k_beidou] = get_seesate(Xs1_beidou,Ys1_beidou,Zs1_beidou,xu1,yu1,zu1,bd_flag);  
    [X_beidou,Y_beidou,Z_beidou] = get_k(num_beidou,k_beidou,Xs_beidou,Ys_beidou,Zs_beidou,bd_flag);
    [XS_beidou,YS_beidou,ZS_beidou] = ecef_g(long,lat,X_beidou,Y_beidou,Z_beidou); % 轉換為地理坐標
    xsg1_beidou = XS_beidou';ysg1_beidou = YS_beidou';zsg1_beidou = ZS_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);
    [gdop] = get_dop(H);
    buffer_dop = [buffer_dop,gdop];
end

plot(buffer_t,buffer_dop_gps,'b-',buffer_t,buffer_dop_galileo,'r-',buffer_t,buffer_dop,'k-');
legend('gps','galileo','gps/galileo/bd');
title('GDOP');
xlabel('time(s)');
ylabel('GDOP');
grid on

a1 = mean()
a2 = mean()
a3 = mean()
b1 = var()
b2 = var()
b3 = var()
c1 = min()
c2 = min()
c3 = min()
d1 = max()
d2 = max()
d3