xpradr.m

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

%                                xpradr.m
%  Scope:   This Matlab program determines pseudorange and accumulated delta range for
%           a specified satellite. WGS-84 constants are used.
%  Usage:   xpradr 
%  Inputs:   - name of the input almanac data file; the default is data file 
%              wk749.dat. Each record contains the eight 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 location for the test (latitude in degrees, longitude
%              in degrees, altitude in meters); the default is data file locat1.dat
%            - elevation angle limit; the default value is 5 degrees
%            - initial simulation time (GPS time of week), in seconds; the 
%              default value is 0.
%            - number of time samples and time step value (if more than one time 
%              step); the default values are 1 and 300, respectively,
%            - selection of value for model flags as follows
%              flags(1) - iono error flag (0 value means iono is not used)
%              flags(2) - tropo error flag (0 value means tropo is not used)
%              flags(3) - SA error flag (0 value means SA error is not used)
%              flags(4) - multipath flag (0 value means multipath is not used)
%              flags(5) - user clock bias flag (0 value means clock bias is not used)
%              flags(6) - earth rotation error flag (0 value means it is not used)
%              flags(7) - pseudorange noise flag (0 value means it is not used)
%              flags(8) - carrier phase flag (0 value means is not generated)
%              flags(9) - delta range noise flag (0 value means is not used)
%            - selection of a satellite index to be plotted (pseudorange and 
%              accumulated delta range)
%            - other parameters are initialized by default
%  Outputs:  - plots of number of visible satellites and repartition of visible 
%              satellites versus time sample number,
%            - plot of all pseudoranges versus time sample number
%            - plot of all accumulated delta ranges versus time sample number
%            - plot of pseudorange and accumulated delta range for a selected satellite
%  References:
%           [1] Parkinson, B. W., Spilker, J. J., Editors, Global Positioning 
%               System: Theory and Applications. AIAA, 1996 (Chapter 12, 
%               Ionospheric effects on GPS by J. A. Klobuchar).
%           [2] Farrell, J., Barth, M., The Global Positioning System and 
%               Inertial Navigation. McGraw Hill, 1998, Appendix E.
%  External Matlab macros used: convcon, eleva, genrn, gpscon, ionoc, ionocon, svpalm, 
%                               tecefgd, tgdecef, uverv, vecefenu, wgs84con, 
%  Last update:  01/19/01
%  Copyright (C) 1999-01 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';

global  alpha  beta  tgeoid  sa  mpath  ucbias  ambig  prnoise  adrnoise  lambda

%  Initialize constants

convcon;
%  global constants used:   deg_rad
gpscon;
%  global constants used:   lambdaL1
lambda = lambdaL1;

%  Read the input file containing satellites data (almanac)

disp('  ');
disp('Enter almanac data - the default is data file wk749.dat');
answer1 = input('Do you want to use the default data file? (y/n)[y] --> ','s');
if isempty(answer1)
   answer1 = yes;
end
disp('  ');
if (strcmp(answer1,yes) == 1)
   ff = 'wk749.dat';
else
   ff = input('Specify the input filename (with extension) -->  ','s');
   disp('  ');
end
almanac = load(ff);
[npt,nmax] = size(almanac);
nr_sat = npt;                %  number of satellites in the almanac
ind_sv = almanac(:,1);       %  id of the svs in the almanac
fprintf('Satellites available in the almanac: \n');
for k = 1:nr_sat
   fprintf('%4.0f',ind_sv(k));
end
fprintf('\n\n');

%  Read the geographic location for the test

disp('Enter geographic location data - the default is the data file  locat1.dat');
answer2 = input('Do you want to use the default data file? (y/n)[y] ','s');
if isempty(answer2)
   answer2 = yes;
end
if (strcmp(answer2,yes) == 1)
   fff = 'locat1.dat';
   hh = load(fff);
   [npt,nmax] = size(hh);
   lat = hh(1);      %  only the first point is selected 
   lon = hh(2);
   alt = hh(3);
else
   lat = input('Enter latitude in degrees  -->  ');
   lon = input('Enter longitude in degrees -->  ');
   alt = input('Enter altitude in meters -->  ');
end

lat_rad   = lat * deg_rad;         %  in radians
lon_rad   = lon * deg_rad;         %  in radians

loc_lla(1) = lat_rad;
loc_lla(2) = lon_rad;
loc_lla(3) = alt;

%  Compute the geographic locations in ECEF

loc_xyz = tgdecef(lat_rad,lon_rad,alt);

%  Specify the elevation angle limit

disp('  ');
disp('Enter elevation 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 angle limit (in degrees) -->  ');
end
ele_lim = ele_limd * deg_rad;   %  elevation angle limit in radians


%  Select initial simulation time (GPS time of week)

disp('  ');
disp('Enter initial simulation time (GPS time of week), in seconds - the default 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)
   tweek = 0.;
else
   tweek = input('Specify initial simulation time (GPS time of week) -->  ');
end
tsim = tweek;

%  Select number of time samples and time step 

disp('  ');
disp('Enter number of time samples - the default value is  10 ');
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 = 10;           %  default - number of time samples
else
   t_sample = input('Enter number of time samples  -->  ' ); 
end

t_incr = 0.;         
if t_sample > 1
   disp('  ');
   disp('Enter time step - the default value is 300. ');
   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 = 300.;         %  default - time step 
   else
      t_incr = input('Enter time step, in seconds  -->  ' ); 
   end
end

% Select model flags and initialize the required data

disp('  ');
disp('Enter selection of flags - the default values are "1" (used) for all');
answer7 = input('Do you want to use the default values? (y/n)[y] ','s');
if isempty(answer7)
   answer7 = yes;
end
if (strcmp(answer7,yes) == 1)
   flags = ones(1,9);            %  default - all terms used
else
   flags(1) = input('Enter  "1"  if iono is used, else "0" --> ');
   flags(2) = input('Enter  "1"  if tropo is used, else "0" --> ');
   flags(3) = input('Enter  "1"  if SA is used, else "0" --> ');
   flags(4) = input('Enter  "1"  if multipath is used, else "0" --> ');
   flags(5) = input('Enter  "1"  if user clock bias is used, else "0" --> ');
   flags(6) = input('Enter  "1"  if earth rotation term is used, else "0" --> ');
   flags(7) = input('Enter  "1"  if pseudorange noise is used, else "0" --> ');
   flags(8) = input('Enter  "1"  if carrier phase is generated, else "0" --> ');
   flags(9) = input('Enter  "1"  if delta range noise is used, else "0" --> '); 
end

if flags(1) == 1
   ionocon
end   
if flags(2) == 1
   tgeoid = load('tgeoid84.dat');
end   
if flags(3) == 1
   load safile
end   
if flags(4) == 1
   load mpfile
end   
if flags(5) == 1
   load ucbfile
end   
if flags(7) == 1
   %  default values for maximum number of satellite = 32
   prnoise = zeros(t_sample,32);
   for kkk = 1:t_sample
      prnoise(kkk,:) = (genrn(32,0.,1.))';
   end   
end   
if flags(8) == 1
   load ambfile
end   
if flags(9) == 1
   %  default values
   adrnoise = zeros(t_sample,32);
   for kkk = 1:t_sample
      adrnoise(kkk,:) = (genrn(32,0.,0.1))';
   end   
end   

%  Start main computation loop - time loop

nr_vis = zeros(1,t_sample);
t_pr = zeros(t_sample,32);        %  32 = maxim satellite id
t_adr = zeros(t_sample,32);
vis_sat = zeros(t_sample,nr_sat);
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  
      temp = almanac(k,2:9);
      psat_temp = (svpalm(tsim,temp))';
      psat(1:3,k) = psat_temp;
   end        

%  Compute elevation and azimuth angles

   for k = 1:nr_sat,      %  cycle for all satellites  
      temp2 = psat(:,k);
      [temp3,temp1,ulos,range] = elevaz(loc_xyz,temp2);
      eangle(k) = temp3;
      aangle(k) = temp1;
   end                    

%  Determine the number of visible satellites and the pseudorange and accumulated 
%  delta range for each visible satellite; if the satellite is not visible the 
%  values are 0.

   for k = 1:nr_sat,
      temp2 = psat(:,k);
      if  eangle(k) >= ele_lim
         nr_vis(kkk) = nr_vis(kkk) + 1;
         vis_sat(kkk,k) = ind_sv(k);
         uelaz(1) = eangle(k);
         uelaz(2) = aangle(k); 
         [pr,adr] = pradr(loc_xyz,loc_lla,uelaz,temp2,tsim,kkk,k,flags);
         t_pr(kkk,ind_sv(k)) = pr;
         t_adr(kkk,ind_sv(k)) = adr;
      end
   end

   tsim = tsim + t_incr;

end      %  end of the time loop (index kkk)

%  Execute the plots

time = 1:t_sample;

%  Number of visible satellites versus time sample

if (t_sample > 1)
   figure(1) 
   subplot(211)
   plot(time,nr_vis,'*')
   grid, axis([1 t_sample 3 12]);
   aa1 =['Time sample number (time step = ' num2str(t_incr) ' seconds, '];
   aa1 = [aa1 'start time tow = ' num2str(tweek) ' seconds)'];
   %xlabel(aa1);
   bb1 = 'Number of visible satellites';
   ylabel(bb1);
   cc1 = [ bb1 '  versus  Time sample number  (almanac  ' ff ')'];
   title(cc1);
   
   subplot(212)
   plot(vis_sat,'*');
   grid
   xlabel(aa1)
   bb2 = 'Repartition of visible satellites';
   ylabel(bb2);
   cc2 = [bb2 '  versus  Time sample number  (almanac  ' ff ')'];
   title(cc2);
   axis([1 t_sample 1  32 ]);
   
   figure(2)
   plot(t_pr);
   xlabel(aa1)
   ylabel('Pseudoranges, in meters');
   title('Pseudoranges  versus  Time sample number');
   grid
   
   figure(3)
   plot(t_adr)
   xlabel(aa1)
   ylabel('Accumulated deltaranges, in meters');
   title('Accumulated deltaranges  versus  Time sample number');
   grid
   
   disp('  ');
   disp('Plot pseudorange and accumulated delta range for a specific satellite');
   answer8 = input('Do you want to execute the plot? (y/n)[y] ','s');
   if isempty(answer8)
      answer8 = yes;
   end
   if (strcmp(answer8,yes) == 1)
      sat_select = input('Enter the satellite id --> ');
      t_pr_sat = t_pr(:,sat_select);
      t_adr_sat = t_adr(:,sat_select);
      ee = num2str(sat_select);
      
      figure(4)
      subplot(211)
   	plot(time,t_pr_sat);
      ylabel('Pseudorange, in meters');
      title(['Pseudorange, for satellite ' ee '  versus  Time sample number']);
      grid

   	subplot(212)
   	plot(time,t_adr_sat);
	   xlabel(aa1)
   	ylabel('Accumulated deltarange, in meters');
      title(['Accumulated deltarange, for satellite ' ee '  versus  Time sample number']);  
      grid
   end
else
   fprintf('\nNumber of visible satellites = %3.0f\n',nr_vis(1));
end   

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