xpelaza.m

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

%                                xpelaza.m
%  Scope:  This Matlab program executes the plots of azimuth-elevation and
%          number of visible satellites for a specified user's position, almanac
%          and time interval, for all satellites in-view; WGS-84 constants are used. 
%  Usage:  xpelaza
%  Inputs:  - reduced almanac data; the default is data file wk749.dat
%             To create an input data file each row contains the following 
%             quantities, for each satellite: 1) satellite id, 2) time of 
%             applicability (toa), in seconds, 3) semi-major axis (a), in meters, 
%             4) eccentricity (e), 5) inclination angle at reference time (I_0),
%             in radians, 6) right ascension at reference time (OMEGA_0), in 
%             radians, 7) argument of perigee (omega), in radians 8) mean anomaly
%             at reference time (M_0), in radians, 9) rate of right ascension
%             (OMEGA_DOT), in radians/second
%           - geographic location(s) for the test; the default is data file
%             xelaz1.dat. The data can be entered from a selected input file or
%             from keyboard by specifying latitude in radians, longitude in
%             radians, and altitude above the ellipsoid in meters
%           - elevation mask angle limit; the default value is 5. degrees
%           - simulation time (GPS time of week); the default value is 0.
%           - time step; the default value is 120 seconds
%           - number of time samples; the default value is 60
%           - name of the output file; the default is xelaz1.out
%           - other parameters are initialized by default
%  Outputs: - azimuth-elevation plot, for all satellites in view, and for 
%             selected 
%           - number of visible satellites plot
%  External Matlab macros used: elevaz, svpalm, tecefgd, tgdecef, vecefenu, 
%                               wgs84con
%  Last update:  12/23/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

clear 
close all

yes = 'y';
disp('  ');

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

%  Read the input file containing satellites data

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)
   fname = 'wk749.dat';
else
   fname = input('Specify the name of the almanac file (with extension) -->  ','s');
   disp('  ');
end
almanac = load(fname);
nr_sat = size(almanac,1);        %  number of satellites in the almanac
fprintf('There are %3.0f satellites in the %s almanac\n\n',nr_sat,fname); 
 
%  Read the geographic location(s) for the test

disp('Enter geographic location data  ');
answer2 = input('Do you want to use a data file? (y/n)[y] --> ','s');
if isempty(answer2)
   answer2 = yes;
end   
disp('  ');
if (~strcmp(answer2,yes) == 1)
   lat(1) = input('Specify the latitude, in radians -->  ');
   disp('  ');
   lon(1) = input('Specify the longitude, in radians -->  ');
   disp('  ');
   alt(1) = input('Specify the altitude above ellipsoid, in meters --> ');
   disp('  ');
   npt = 1;
else
   disp('Enter name of the location data file - the default is xelaz1.dat');
   answer3 = input('Use the default data file? (y/n)[y] --> ','s');
   if isempty(answer3)
      answer3 = yes;
   end   
   disp('  ');
   if (strcmp(answer3,yes) == 1)
      fff = 'xelaz1.dat';
   else
      fff = input('Specify the filename (with extension) --> ','s');
      disp('  ');
   end
   hh = load(fff);
   if (size(hh,1) ~= 1)
      hh_temp(1,:) = hh;
      clear hh;
      hh = hh_temp;
   end   
   % select first location
   lat   = hh(1) * deg_to_rad;         %  in radians
   lon   = hh(2) * deg_to_rad;         %  in radians
   alt   = hh(3);                      %  in meters
end

%  Determine the geographic location(s) in ECEF

user = tgdecef(lat,lon,alt);

%  Specify the elevation mask angle limit

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

%  Select initial simulation time (GPS time of week)

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

%  Initialize time step and number of time samples

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

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

%  Execute the main loop computation and execute the plot

kk = 0;
for kkk = 1:t_sample         %  cycle for all time samples
   
   clear elangle azangle
   for k = 1:nr_sat,        %  cycle for all satellites  

%     Compute satellite position in ECEF

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

%     Compute elevation and azimuth angles

      [el_temp,az_temp,u_temp,r_temp] = elevaz(user,psat_temp);
      elangle(k) = el_temp;
      azangle(k) = az_temp;

   end                    %  end of the satellite loop (index  k)

%  Determine the number of visible satellites, elevation and azimuth angles

   clear el_vis az_vis
   nr_vis = 0;
   for k = 1:nr_sat,
      if  elangle(k) >= ele_lim
         nr_vis = nr_vis + 1;
         el_vis(nr_vis) = elangle(k);
         az_vis(nr_vis) = azangle(k);
         if  (kkk == 1)
            ind_vis_first(nr_vis) = almanac(k,1);
         end   
      end
   end
   nrvis(kkk) = nr_vis;
   
   dist = (90.*ones(1,nr_vis) - el_vis*rad_to_deg)/100.; 
   for i = 1:nr_vis,
      kk = kk + 1;
      svx(kk) = dist(i)*cos(az_vis(i));
      svy(kk) = dist(i)*sin(az_vis(i));
   end  
   if  (kkk == 1)
      svx_first = svx;
      svy_first = svy;
      nr_vis_first = nr_vis;
   end   
   
   tsim = tsim + t_incr;

end                          %  end of the time loop (index kkk)

%  Execute the azimuth-elevation plot

figure(1)
aa = ' Time sequences   (Time step = ';
aa = [num2str(t_sample)  aa  num2str(t_incr)  ' seconds)'];
cc = ['Azimuth-Elevation plot for ' aa];
h0 = [0  0]; 
h9 = [-0.9  0.9];
an1 = [0:120]*pi/60;
an2 = [0:100]*pi/50;
an3 = [0:80]*pi/40;
an4 = [0:60]*pi/30;
an5 = [0:32]*pi/16;
an6 = [0:24]*pi/12;
plot(0.9*sin(an1),0.9*cos(an1),'.',0.8*sin(an1),0.8*cos(an1),'.',...
     0.7*sin(an2),0.7*cos(an2),'.',0.6*sin(an2),0.6*cos(an2),'.',...
     0.5*sin(an3),0.5*cos(an3),'.',0.4*sin(an3),0.4*cos(an3),'.',...
     0.3*sin(an4),0.3*cos(an4),'.',0.2*sin(an5),0.2*cos(an5),'.',...
     0.1*sin(an6),0.1*cos(an6),'.',h9,h0,'-k',h0,h9,'-k')
hold on
plot(svx,svy,'.k'); 
hold on
plot(svx_first,svy_first,'x'); 
title(cc)
axis('square');
axis('off');
text(1.,0,'East')
text(-1.15,0,'West')
text(-0.07,0.96,'North')
text(-0.07,-0.96,'South')
text(0.13,-0.22,'60')
text(0.31,-0.46,'30')
text(0.55,-0.68,'0')
for k = 1:nr_vis_first
   text(svx_first(k),svy_first(k),num2str(ind_vis_first(k)))
end 
xlab = ['Location:  lat = ' num2str(lat) ' rad, lon = ' num2str(lon) ];
xlab = [xlab  ' rad, alt = ' num2str(alt) ' m' ];
text(-1.15,-1.15,xlab)
text(0.,-1.05,['Simulation start time (tow) = ' num2str(tow) ' sec.']);
text(-1.15,-01.05,['Almanac = ' fname ]);

% Execute the number of visible satellites plot

x = 1:1:t_sample;              
figure(2)
f2 = ['Time sequence  (Time step = ' num2str(t_incr)];
f2 = [f2 ' sec, Simulation start time (tow) = ' num2str(tow) ' sec'];
f2 = [f2 ', almanac = ' fname ')'];
xlab = [' (User: lat=' num2str(lat) 'rad,lon=' num2str(lon) ];
xlab = [xlab  'rad,alt=' num2str(alt) 'm)' ];
cc = ['Number of visible satellites versus Time sequence ' xlab];
ee = 'Number of visible satellites';
plot(x,nrvis), ...
title(cc), ylabel(ee), xlabel(f2), grid
axis([0 length(x) 4 12])

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