xsvpcomp.m

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

%                                xsvpcomp.m
%  Scope:   This MATLAB program computes RSS between ECEF satellite position 
%           based on satellite ephemeris data and the corresponding ECEF 
%           satellite position based on almanac data; WGS-84 constants are used.
%  Usage:   xsvcomp
%  Inputs:  - tsim, GPS system time at time of transmission, i.e. GPS time
%             corrected for transit time (range/speed of light), in seconds
%           - satellite ephemeris data from a specified data file; default 
%             data are given. The ephemeris satellite data are specified in
%             the following order:
%             - satellite number
%             - toe, reference time ephemeris, in seconds         
%             - smaxis (a), satellite semi-major axis, in meters
%             - ecc (e), satellite eccentricity  
%             - izero (I_0), inclination angle at reference time, in radians 
%             - razero (OMEGA_0), right ascension at reference time, in radians 
%               (longitude of ascending node of orbit plane at weekly epoch)
%             - argper (omega), argument of perigee, in radians
%             - mzero (M_0), mean anomaly at reference time, in radians      
%             - radot (OMEGA_DOT), rate of right ascension, in radians/second 
%             - deln (delta_n), mean motion difference from computed value, in 
%               radians/second 
%             - idot (I_DOT), rate of inclination angle, in radians/second 
%             - cic, amplitude of the cosine harmonic correction term to the
%               angle of inclination, in radians
%             - cis, amplitude of the sine harmonic correction term to the
%               angle of inclination, in radians
%             - crc, amplitude of the cosine harmonic correction term to the
%               orbit radius, in meters
%             - crs, amplitude of the sine harmonic correction term to the
%               orbit radius, in meters
%             - cuc, amplitude of the cosine harmonic correction term to the
%               argument of latitude, in radians
%             - cus, amplitude of the sine harmonic correction term to the
%               argument of latitude, in radians
%           - satellite almanac data from a specified data file; default 
%             data are given. The almanac satellite data are specified in
%             the following order:
%             - satellite number
%             - toa, reference time applicability, in seconds         
%             - smaxis (a), satellite semi-major axis, in meters
%             - ecc (e), satellite eccentricity  
%             - izero (I_0), inclination angle at reference time, in radians 
%             - razero (OMEGA_0), right ascension at reference time, in radians 
%               (longitude of ascending node of orbit plane at weekly epoch)
%             - argper (omega), argument of perigee, in radians
%             - mzero (M_0), mean anomaly at reference time, in radians      
%             - radot (OMEGA_DOT), rate of right ascension, in radians/second 
%           - name of the input data files (optional)
%           - number of time steps (integer)
%           - time step value, in seconds; initial time is tsim
%           - name of the output files (optional)
%  Outputs: - plot of the RSS between the satellite ephemeris position and 
%             satellite almanac position
%           - 2 output data files (optional) storing the ECEF satellite position 
%             based on ephemeris/almanac; each record contains the following data:
%             time step number, satellite ECEF x-component, satellite ECEF 
%             y-component, satellite ECEF z-component
%  External Matlab macros used:  svpalm, svpeph, wgs84con
%  Last update:  11/09/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';

%  Constants and initialization

wgs84con;
% global constants used: gravpar, rot_rate

%  Specify the ephemeris/almanac data for the satellite

disp('  ');
answer1 = input('Do you want to use the default data? (y/n)[y] ','s');
if isempty(answer1)
   answer1 = yes;
end   
disp('  ');
if (strcmp(answer1,yes) == 1)
   tsim   =  375300.;                %  GPS system time at time of transmission,
                                     %  i.e. GPS time corrected for transit time 
                                     %  (range/speed of light), in seconds
   eph(1) = 1;                       %  default value                                  
   eph(2) =  381600.;                %  reference time ephemeris (toe), in seconds     
   eph(3) =  26560092.2173461914;    %  semi-major axis (a), in meters
   eph(4) =  0.008817672729492;      %  satellite eccentricity
   eph(5) =  0.95694279737716;       %  inclination angle at reference time (I_0), 
                                     %  in radians
   eph(6) =  1.08030945307641;       %  right ascension at reference time (OMEGA_0), 
                                     %  in radians 
   eph(7) =  -0.65993493835624;      %  argument of perigee (omega), in radians
   eph(8) =  1.295751523149107;      %  mean anomaly at reference time (M_0), in 
                                     %  radians
   eph(9) =  -7.840326701505655e-09; %  rate of right ascension (OMEGA_DOT), in 
                                     %  radians/second
   eph(10) =  0.0e-9;                %  mean motion difference from computed 
                                     %  value (delta_n), in radians/second
   eph(11) =  0.0e-9;                %  rate of inclination angle (I_DOT), in 
                                     %  radians/second
   eph(12) =   3.911554601483828e-008;   %  cic,  in radians
   eph(13) =   9.313226288825713e-009;   %  cis,  in radians
   eph(14) =  172.0625   ;               %  crc,  in meters  
   eph(15) =  21.90625    ;              %  crs,  in meters  
   eph(16) =   3.706663761359739e-007;   %  cuc,  in radians
   eph(17) =   3.371387574121309e-006;   %  cus,  in radians
   
   alm(1:9) = eph(1:9);              % almanac data <= truncated ephemeris data
else
   tsim = input('Specify the tsim, e.g. 238205., --> ');
   disp('  ');
   fname1 = input('Specify the input data filename for the ephemeris data --> ','s');
   tt = load(fname1);
   [nrow,ncol] = size(tt);
   if  nrow ~= 17
      disp('XSVPCOMP - Error: check the input ephemeris data ');
      disp('  ');
      break;
   end
   eph = tt(1,:);    %  for the first satellite
   
   disp('  ');
   fname2 = input('Specify the input data filename for the almanac data --> ','s');
   ttt = load(fname2);
   [nrow,ncol] = size(ttt);
   if  nrow ~= 9
      disp('XSVPCOMP - Error: check the input almanac data ');
      disp('  ');
      break;
   end
   alm = ttt(1,:);   %  for the first satellite (it is assumed that it is the same,
                     %  otherwise a selection can be inserted)
end

nrstep = input('Specify number of steps  -- > ');
disp('  ');
timestep = input('Specify the time-step, in seconds, e.g. 10. --> ');
disp('  ');
tsiminit = tsim;

edata = eph(2:17);
adata = alm(2:9);

for kk = 1:nrstep

%  Compute the ECEF satellite position based on ephemeris and almanac data 
%  (truncated ephemeris data)

   pse = svpeph(tsim,edata);
   psekeep(kk,:) = pse;
   psa = svpalm(tsim,adata);
   psakeep(kk,:) = psa;
   svpdif(kk) = norm(pse-psa);
   tsim = tsim + timestep;

end

%  Execute the plot 

time = 1:nrstep;
y = svpdif';
aa = 'RSS Satellite position difference based on ephemeris and almanac';
aa = [aa '  versus  Time step number']; 
bb = 'RSS Satellite position difference, in meters'; 
cc = ['Time step number;  initial time (tow) = ' num2str(tsiminit) ];
cc = [cc ' second(s) and time step = ' num2str(timestep) ' second(s)'];
ymean = mean(y);
ystd = std(y);
yrms = rms(y);
temp1 = ['mean = ',num2str(ymean)];
temp2 = ['st.dev. = ',num2str(ystd)];
temp3 = ['rms = ',num2str(yrms)];
g = [temp1,' ;   ',temp2,' ;   ',temp3];
disp('Place the text with statistics on the graph; ');
temp = input('      Press <Enter> or <Return> key to start ... ');
figure(1)
plot(time,y),grid, ...
title(aa), ylabel(bb),xlabel(cc)
gtext(g)           %  mouse placement of the text on a graph
disp('   ');

%  Save the generated data into a specified file or displayed on screen

temp = input('Press <Enter> or <Return> key to continue ... ');
disp('  ');
answer2 = input('Do you want to save the generated data? (y/n)[y] ','s');
if  isempty(answer2)
   answer2 = yes;
   disp('  ');
end

if strcmp(answer2,yes) == 1

%  Specify the output files and save the data (time, and x/y/z satellite
%  position at each time step)

   f3 = input('Specify the output filename based on ephemeris computation --> ','s');
   disp('  ');
   f4 = input('Specify the output filename based on almanac computation --> ','s');
   for kk = 1:nrstep
      fprintf(f3,'%d %f %f %f\n',kk,psekeep(kk,1),psekeep(kk,2),psekeep(kk,3));
   end
   for kk = 1:nrstep
      fprintf(f4,'%d %f %f %f\n',kk,psakeep(kk,1),psakeep(kk,2),psakeep(kk,3));
   end

end

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