xraimda.m

GPS TOOLBOX包含以下内容： 1、GPS相关常量和转换因子； 2、角度变换； 3、坐标系转换： &#61656 点变换； &#61656 矩阵变换； &#61656 向量变换

else
   t_sample = input('Specify number of time samples -->  ');
end
nr_sp_time = npt * t_sample;           % total number of space time points

disp('  ');
disp('Enter time increment value, in seconds - the default value is 1800');
answer6 = input('Do you want to use the default value? (y/n)[y] ','s');
if isempty(answer6)
   answer6 = yes;
end   
if (strcmp(answer6,yes) == 1)
   t_incr = 1800.;           %  30 minutes - time increment
else
   t_incr = input('Specify time increment, in seconds -->  ');
end

%  Specify the output files (optional)

disp('  ');
answer7 = input('Do you want to save the statistics? (y/n)[n] ','s');
disp('  ');
if isempty(answer7)
   answer7 = 'n';
   pf1 = 1;               %  on sceen
end   
if (strcmp(answer7,yes) == 1)
   pf1 = input('Specify the statistics output filename (with extension) --> ','s');
   disp('  ');
end
answer8 = input('Do you want to save the unavailable cases? (y/n)[n] ','s');
disp('  ');
if isempty(answer8)
   answer8 = 'n';
   pf2 = 0;
end   
if (strcmp(answer8,yes) == 1)
   pf2 = input('Specify the RAIM unavailable output filename (with extension) --> ','s');
   disp('  ');
end

%  Initialization of arrays and counters 

count = zeros(1,12);
inad = zeros(1,12);
psat = zeros(3,nr_sat);
ulos = zeros(3,nr_sat);
eangle = zeros(1,nr_sat);
ind_vis = zeros(1,nr_sat);
hpl = zeros(1,nr_sp_time);

count_a = 0;            %  availability counter       
count_una = 0;          %  unavailability counter 
kt = 1;

%  Start main computation loop

for kkk = 1:t_sample    %  cycle for all time samples

   fprintf('*****   Computation in progress for time sample = %4.0f\n',kkk);

%  Determine satellite position

   for k = 1:nr_sat,       %  cycle for all satellites  

%     Compute satellite position in ECEF

      temp = almanac(k,2:9);
      psat_temp = (svpalm(tsim,temp))';
      psat(1:3,k) = psat_temp;

   end        %  end of the satellite position loop   (index  k)

   clear ind_vis
   for kk = 1:npt          %  cycle for all locations

      user = user_ecef(:,kk);

      for k = 1:nr_sat,      %  cycle for all satellites  
   
%     Compute elevation angle

         temp2 = psat(:,k);
         [temp3,temp1] = eleva(user,temp2);
         eangle(k) = temp3;
         ulos(:,k) = vecefenu(temp1,lat(kk),lon(kk));

      end                    %  end of the satellite loop (index  k)
            
%  Determine the number of visible satellites and save the index and elevation
%  angle for all visible satellites for a specific location

      clear  htransp
      nr_vis = 0;
      for k = 1:nr_sat,
         if  eangle(k) >= ele_lim
            nr_vis = nr_vis + 1;
            ind_vis(nr_vis) = almanac(k,1);
            if nr_vis == 1
               htransp = [ulos(:,k)' 1.]';
            else
               htransp = [htransp [ulos(:,k)' 1.]'];
            end
         end
      end

%  Determine the slopemax for each location

      if  nr_vis <= 4
         slmax = 1.e+10;
         svcrit = 0;
      else
         clear  slope  a  b 
         a = inv(htransp * htransp') * htransp;
         b = htransp' * a;
         for k = 1:nr_vis
            slope(k) = sqrt( (a(1,k) * a(1,k) + a(2,k) * a(2,k)) / (1. - b(k,k)) );
         end
         slmax = max(slope);

%  Determine the index of the critical satellite (satellite with the maximum
%  slope)

          for k = 1:nr_vis
            if (slmax == slope(k))
               svcrit = ind_vis(k);
            end
         end
      end

%  Determine and store horizontal protection level

      if  nr_vis <= 4
         prot_rad = 1.e+10;
      else
         prot_rad = slmax * sigma * pbiasb(nr_vis - 4);
      end
      hpl(kt) = prot_rad;
      kt = kt + 1;

%  Determine and count the availability/unavailability repartition based on 
%  number of visible satellites, and store available/unavailable cases
% 
      if prot_rad > prot_lim
         count_una = count_una + 1;
         inad(nr_vis) = inad(nr_vis) + 1;   
         if pf2 ~= 0
            fprintf(pf2,'%3.0f %3.0f %3.0f %3.0f %10.3f %10.3f\n',...
                         kkk,kk,nr_vis,svcrit,prot_rad,prot_lim);
         end
      else
         count_a = count_a + 1;
         hplav(count_a) = prot_rad;
      end  

%  Count the satellite visibility repartition

      count(nr_vis) = count(nr_vis) + 1;

   end      %  end of the location loop (index kk)

   tsim = tsim + t_incr;

end      %  end of the time loop (index kkk)

ad_per = count_a * 100. / nr_sp_time;
inad_per = count_una * 100. / nr_sp_time;

fprintf(pf1,'\n****************************************************************\n');
fprintf(pf1,'\n   RAIM  FAULT  DETECTION  AVAILABILITY  TEST  SUMMARY  \n\n');
fprintf(pf1,'%s\n',pof_type);
fprintf(pf1,'-  horizontal alarm limit (HAL) (meters): %9.3f\n',prot_lim);
fprintf(pf1,'-  GPS receiver mode:  %s\n',gpstype);
fprintf(pf1,'-  almanac: %s \n',ff);
fprintf(pf1,'-  number of satellites in constellation = %3.0f \n',nr_sat);
fprintf(pf1,'-  number of locations:  %4.0f\n',npt);
fprintf(pf1,'-  number of time samples: %4.0f\n',t_sample);
fprintf(pf1,'-  time increment (seconds):  %5.0f\n',t_incr);
fprintf(pf1,'-  std. dev. of pseudorange error (meters): %8.3f\n',sigma);
fprintf(pf1,'-  elevation mask angle limit (degrees): %8.4f  \n',ele_limd);
fprintf(pf1,'-  probability of false alarm: %14.7e\n',pfa);
fprintf(pf1,'-  probability of missed detection: %14.7e\n',pmd);
fprintf(pf1,'-  total space-time sample points: %6.0f\n\n',nr_sp_time);
fprintf(pf1,'RAIM fault detection available   =  %5.0f,  percentage: %8.5f\n', ...
             count_a,ad_per);
fprintf(pf1,'RAIM fault detection unavailable =  %5.0f,  percentage: %8.5f\n\n', ...
             count_una,inad_per);
fprintf(pf1,'Repartition of RAIM fault detection available/unavailable cases:\n');
fprintf(pf1,'- when < 5 visible satellites: %4.0f available, %4.0f unavailable\n',...
               (count(4)-inad(4)),inad(4));
for k = 5:12
   fprintf(pf1,'- when %3.0f visible satellites: %4.0f available, %4.0f unavailable\n',...
               k,(count(k)-inad(k)),inad(k));
end
fprintf(pf1,'\n****************************************************************\n');
fprintf(pf1,'\n             SATELLITE  VISIBILITY  STATISICS    \n');
fprintf(pf1,'\nNumber of cases with < 5 visible satellites  = %4.0f\n',count(4));
for k = 5:12
   fprintf(pf1,'Number of cases with %3.0f visible satellites  = %4.0f\n',k,count(k));
end   
fprintf(pf1,'\n****************************************************************\n');

%  Execute a bar graph presenting number of cases versus number of
%  visible satellites

disp('  ');
disp('Bar plot - Number of cases versus Number of visible satellites'); 
answer9 = input('Do you want to execute the bar graph? (y/n)[y] ','s');
if isempty(answer9)
   answer9 = yes;
end   
if (strcmp(answer9,yes) == 1)
   x = 4:1:12;
   y = [count(4),count(5),count(6),count(7),count(8),count(9),count(10),...
         count(11),count(12)];
   figure(1)
   bar(x,y); grid, ...
   xlabel('Number of visible satellites');
   ylabel('Number of cases');
   title('Satellite visibility repartition');
end

%  Execute a histogram plot presenting the repartition of horizontal protection 
%  levels of all available cases

disp('  ');
disp('Histogram plot - Number of cases versus Horizontal protection level'); 
answer10 = input('Do you want to execute the histogram plot? (y/n)[y] ','s');
if isempty(answer10)
   answer10 = yes;
end   
if (strcmp(answer10,yes) == 1)
   xb = 0:50:(prot_lim+150);
   xmin = min(hplav);
   xmax = max(hplav);
   xmean = mean(hplav);
   xmedian = median(hplav);
   %xtext = ['min = ' num2str(xmin) ', max = ' num2str(xmax) ', median = '];
   %xtext = [xtext num2str(xmedian) ', mean = ' num2str(xmean) ' meters'];
   xtext = ['min = ' num2str(xmin) ', max = ' num2str(xmax) ];
   xtext = [xtext ', mean = ' num2str(xmean) ' meters'];
   figure(2)
   hist(hplav,xb), grid, ...
   xlabel('Horizontal protection level, in meters');
   ylabel('Number of cases');
   title([pof_mode 'Horizontal protection level histogram of available cases']);
   v = axis;
   axis([ v(1) v(2) v(3)  (v(4)+50) ]);
   text(0,v(4)+25,xtext)
end

disp('  ');
disp('End of the program  XRAIMDA');
disp('  ');