xraimda.m

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

%                                   xraimda.m
%  Scope:  This Matlab program determines the Receiver Autonomous Integrity 
%          Monitoring (RAIM) fault detection availability based on a constant 
%          false alarm rate algorithm [1], [2]. WGS-84 constants are used.
%  Usage:  xraimda 
%  Inputs:   - almanac data; the default is datafile svprime.dat. Each record 
%              contains the nine data related to a satellite in the following order: 
%              1) satellite id, 
%              2) satellite time of applicability (toa), in seconds,
%              3) satellite semi-major axis (a), in meters, 
%              4) satellite eccentricity (e),
%              5) satellite orbital inclination (I_0), in radians,
%              6) satellite right ascension at toa (OMEGA_0), in radians, 
%              7) satellite argument of perigee (omega), in radians,
%              8) satellite mean anomaly (M_0), in radians,
%              9) satellite rate of right ascension (OMEGA_DOT), radians/second
%            - geographic locations for the test; the default is datafile
%              mops24.dat. Each record contains the data related to a location
%              in the following order: latitude (degrees), longitude (degrees), 
%              altitude (meters)
%            - elevation mask angle limit in degrees; the default is 5. degrees
%            - GPS receiver type: Standard Position System (SPS) with SA on,
%              SPS without SA, or user's selected standard deviation of 
%              measurement error, in meters
%            - phase-of-flight (Non-Precision Approach, Terminal, En_Route, Oceanic
%              and user defined HAL)
%            - probability of false alarm, probability of missed detection, 
%              and the corresponding pbias thresholds; default values are 
%              provided
%            - initial simulation time (GPS time of week); the default value is 0.
%            - number of time steps; default value is 48 [1]
%            - time increment, in seconds; default value is 1800 [1]
%            - name of the statistics output file, the default is the screen
%            - name of the output file containing unavailable cases (optional)
%            - selection of the bar graph of visible satellites (optional)
%            - selection of the histogram plot of horizontal protection level (HPL)
%              repartition (optional)
%            - other parameters are initialized by default
%  Outputs:  - table containing RAIM fault detection availability test summary
%              and satellite visibility statistics (see file pf1 or on sceen)
%            - file containing data for RAIM unavailable cases (time, location,
%              number of visible satellites, critical satellite index, horizontal
%              protection level and horizontal alarm limit) (see file pf2)
%            - bar graph presenting repartition of visible satellites (optional)
%            - histogram plot presenting repartition of horizontal protection level 
%              for all available cases (optional)
%  External Matlab macros used: eleva, svpalm, tgdecef, uverv, vecefenu, wgs84con
%  References: 
%          [1] Anonym., Minimum operational performance standards for airborne
%              supplemental navigation equipment using Global Positioning System
%              (GPS). Document No. RTCA/DO-208, July 1991, Prepared by SC-159.
%          [2] Brown, R. G., A baseline RAIM scheme and a note on the
%              equivalence of three RAIM methods. Proceedings of the National
%              Technical Meeting, Institute of Navigation, San Diego, CA,
%              Jan. 27-29, 1992, pp. 127-137.
%          [3] Brown, R. G., GPS RAIM: Calculation of thresholds and protection 
%              radius using Chi-square methods - A geometric approach. RTCA 
%              Paper No. 491-94/SC159-584, November 1994.
%  Last update: 11/23/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';

%  Constants

deg_to_rad = pi / 180.; %  degrees to radians transformation

%  Read the input file containing satellites data

disp('  ');
disp('Enter almanac data file - the default is datafile svprime.dat');
answer1 = input('Do you want to use the default datafile? (y/n)[y] ','s');
if isempty(answer1)
   answer1 = yes;
end   
if (strcmp(answer1,yes) == 1)
   ff = 'svprime.dat';
else
   ff = input('Specify the input filename (with extension) -->  ','s');
end
almanac = load(ff);
[nr_sat,ncol] = size(almanac);         %  number of satellites used
if  ncol ~= 9
   disp('Error - XRAIMDA; check the almanac input file');
   disp('  ');
   return
end   

%  Read the geographic location for the test

disp('  ');
disp('Enter geographic location data - the default is the datafile mops24.dat');
answer2 = input('Do you want to use the default datafile? (y/n)[y] ','s');
if isempty(answer2)
   answer2 = yes;
end   
if (strcmp(answer2,yes) == 1)
   fff = 'mops24.dat';
else
   fff = input('Specify the input filename (with extension) --> ','s');
end
hh = load(fff);
[npt,ncol] = size(hh);     %  npt - number of locations
if  ncol ~= 3
   disp('Error - XRAIMDA; check the location input file');
   disp('  ');
   return
end   
 
lat   = hh(1:npt,1) * deg_to_rad;         %  in radians
lon   = hh(1:npt,2) * deg_to_rad;         %  in radians
alt   = hh(1:npt,3);                      %  in meters

%  Compute the geographic locations in ECEF

for kk = 1:npt
   temp = tgdecef(lat(kk),lon(kk),alt(kk));
   user_ecef(:,kk) = temp;
end

%  Specify the elevation mask angle limit

disp('  ');
disp('Enter elevation mask angle limit - the default value is 5. degrees');
answer3 = input('Do you want to use the default value? (y/n)[y] ','s');
if isempty(answer3)
   answer3 = yes;
end   
if (strcmp(answer3,yes) == 1)
   ele_limd = 5.;
else
   ele_limd = input('Specify the elevation mask angle limit (in degrees) --> ');
end
ele_lim = ele_limd * deg_to_rad;    %  elevation angle limit in radians

%  Specify GPS receiver type and select the corresponding sigma (in meters)

disp('  ');
disp('Enter GPS receiver type :');
disp('      1  -->  Standard Position System (SPS) with SA On');
disp('      2  -->  Standard Position System (SPS) with SA Off');
disp('      3  -->  User defined Std. Dev. of measurement errors');
gps_type = input('Make the selection  -->  ');
if  (gps_type ~= 1) & (gps_type ~= 2) & (gps_type ~= 3)
   disp('Error - XRAIMDA; check the input for gps_type ');
   disp('  ');
   return
end
if  gps_type == 1
   sigma = 33.3;
   gpstype = 'SPS with SA On';
elseif  gps_type == 2
   sigma = 6.5;
   gpstype = 'SPS with SA Off';
elseif  gps_type == 3
   sigma = input('Specify Std. Dev. of measurement errors --> ');
   gpstype = 'User defined';  
end

%  Select phase-of-flight and the corresponding horizontal alarm limit

disp('  ');
disp('Select phase-of-flight:');
disp('      1  -->  Non-Precision Approach (NPA) ');
disp('      2  -->  Terminal  ');
disp('      3  -->  En Route ');
disp('      4  -->  Oceanic ');
disp('      5  -->  User defined ');
pof = input('Make the selection  -->  ');
if  (pof ~= 1) & (pof ~= 2) & (pof ~= 3) & (pof ~= 4) & (pof ~= 5)
   disp('Error - XRAIMDA; check the input for phase-of-flight ');
   disp('  ');
   return
end
if  pof == 1
   prot_lim = 555.6;        % HAL for  NPA  (in meters)
   pof_type = 'For Non-Precision Approach with All-In-View satellites:';
   pof_mode = 'Non-Precision Approach phase of flight - ';
elseif  pof == 2
   prot_lim = 1852.;        % HAL for  Terminal  (in meters)
   pof_type = 'For Terminal with All-In-View satellites:';
   pof_mode = 'Terminal phase of flight - ';
elseif  pof == 3
   prot_lim = 3704.;        % HAL for  En Route (in meters)
   pof_type = 'For En Route with All-In-View satellites:';
   pof_mode = 'En Route phase of flight - ';   
elseif  pof == 4
   prot_lim = 7408.;        % HAL for  Oceanic  (in meters)
   pof_type = 'For Oceanic with All-In-View satellites:';
   pof_mode = 'Oceanic phase of flight - ';   
elseif  pof == 5
   prot_lim = input('Specify Horizontal Alarm Limit (HAL), in meters --> ');
   pof_type = 'For User defined HAL with All-In-View satellites:';
   pof_mode = ['For HAL = ' num2str(prot_lim) ' meters - '];   
end

%  Select the tresholds table pbiasb for specified probability of false 
%  alarm (pfa) and probability of missed detection (pmd)

disp('  ');
disp('Enter values for pfa, pmd and pbiasb:');
disp('      1  -->  For pfa = 6.6667e-5  and  pmd = 0.001 ');
disp('      2  -->  For pfa = 3.3333e-7  and  pmd = 0.001 ');
disp('      3  -->  Enter the user defined values ');
pfa_pmd = input('Make the selection  -->  ');
if  (pfa_pmd ~= 1) & (pfa_pmd ~= 2) & (pfa_pmd ~= 3)
    error('Error - XRAIMDA; check the input for pfa_pmd ');
end
if  pfa_pmd == 1
   pfa = 6.6667e-5;
   pmd = 0.001;
   pbiasb = [7.0781 7.3883 7.6091 7.7880 7.9413 8.0770 8.1994 8.3116]; 
elseif pfa_pmd == 2
   pfa = 3.3333e-7;
   pmd = 0.001;
   pbiasb = [8.1937 8.4785 8.6874 8.8597 9.0090 9.1423 9.2634 9.3759];
else
   pfa = input('Enter probability of false alarm (pfa) --> ');
   pmd = input('Enter probability of missed detection (pmd) --> ');
   pbiasb(1) = input('Enter the value for  5 satellites, pbiasb(1) -->  ');
   pbiasb(2) = input('Enter the value for  6 satellites, pbiasb(2) -->  ');
   pbiasb(3) = input('Enter the value for  7 satellites, pbiasb(3) -->  ');
   pbiasb(4) = input('Enter the value for  8 satellites, pbiasb(4) -->  ');
   pbiasb(5) = input('Enter the value for  9 satellites, pbiasb(5) -->  ');
   pbiasb(6) = input('Enter the value for 10 satellites, pbiasb(6) -->  ');
   pbiasb(7) = input('Enter the value for 11 satellites, pbiasb(7) -->  ');
   pbiasb(8) = input('Enter the value for 12 satellites, pbiasb(8) -->  ');
end

%  Select initial simulation time (GPS time of week)

disp('  ');
disp('Enter initial simulation time (GPS time of week) - the default value is 0.');
answer4 = input('Do you want to use the default value? (y/n)[y] ','s');
if isempty(answer4)
   answer4 = yes;
end   
if (strcmp(answer4,yes) == 1)
   tsim = 0.;
else
   tsim = input('Specify initial simulation time (GPS time of week) -->  ');
end

%  Select number of time samples and time increment

disp('  ');
disp('Enter number of time samples - the default value is 48');
answer5 = input('Do you want to use the default value? (y/n)[y] ','s');
if isempty(answer5)
   answer5 = yes;
end   
if (strcmp(answer5,yes) == 1)
   t_sample = 48;          %  number of time samples