xelaz.m

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

%                                xelaz.m
%  Scope:  This Matlab program determines elevation and azimuth angles for
%          specified users and a specified time interval, for all satellites 
%          in-view, and saves into a specified file; WGS-84 constants are used.
%  Usage:  xelaz1
%  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
%           - GPS time of week (simulation time); the default value is 0.
%           - time step; the default value is 120 seconds
%           - number of time samples; the default value is 360
%           - name of the output file; the default is xelaz1.out
%           - other parameters are initialized by default
%  Outputs: - a file containing for each time step and for each visible
%             satellite the following (the default file is xelaz1.out):
%             1) time step number, 2) user's position index, 3) number of
%             visible satellites, 4) id of visible satellite, 5) elevation
%             angle (rad.), 6) azimuth angle (rad.).
%             The output file can be used as input file for the program
%             xpgelaz.m, which executes plots related to above-listed
%             quantities.
%  Remarks:   If the structures of the input/output data files are changed
%             then the statements related to reading and/or writing should
%             also be changed.
%  External Matlab macros used: elevaz, svpalm, tecefgd, tgdecef, vecefenu, 
%                               wgs84con
%  Last update:  11/12/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

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

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

%  Read the input file containing satellites data

disp('Enter almanac data - the default is data file wk749.dat');
answer = input('Do you want to use the default data file? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,yes) == 1)
   fname = 'wk749.dat';
else
   fname = input('Specify the name of the input file (with extension) -->  ','s');
   disp('  ');
end
almanac = load(fname);
nr_sat = size(almanac,1);        %  number of satellites used
 
%  Read the geographic location(s) for the test

disp('Enter geographic location data  ');
answer = input('Do you want to enter data from keyboard? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,yes) == 1)
   lat(1) = input('Specify the latitude, in radians -->  ');
   disp('  ');
   long(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 data file - the default is xelaz1.dat');
   answer = input('Use the default data file? (y/n)[n] --> ','s');
   disp('  ');
   if (strcmp(answer,yes) == 1)
      fff = 'xelaz1.dat';
   else
      fff = input('Specify the filename (with extension) --> ','s');
      disp('  ');
   end
   hh = load(fff);
   npt = size(hh,1);
   lat   = hh(1:npt,1) * deg_to_rad;         %  in radians
   long  = hh(1:npt,2) * deg_to_rad;         %  in radians
   alt   = hh(1:npt,3);                      %  in meters
end

%  Determine the geographic location(s) in ECEF

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

%  Specify the elevation angle limit

disp('Enter elevation angle limit - the default value is 5. degrees');
answer = input('Do you want to use the default value? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,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 GPS time of week (initial simulation time)

disp('Enter simulation time (GPS time of week) - the default value is 0. ');
answer = input('Do you want to use the default value? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,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 ');
answer = input('Do you want to use the default value? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,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 180 ');
answer = input('Do you want to use the default value? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,yes) == 1)
   t_sample = 180;           %  number of time samples
else
   t_sample = input('Specify the number of time samples -->  ');
   disp('  ');
end

%  Specify the output filename

disp('Enter the name of the output file ');
answer = input('Do you want to use the default xelaz1.out? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer,yes) == 1)
   pf1 = 'xelaz1.out';
else
   pf1 = input('Specify the output filename (with extension) -->  ','s');
   disp('  ');
end

%  Execute the main loop computation

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

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

   for kk = 1:npt            %  cycle for all locations

      user = user_ecef(:,kk);

      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(k,1:3) = psat_temp';          %  save satellite position

%        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 and save the index, elevation
%     angle and azimuth angle for all visible satellites for a specific location

      nr_vis(kk) = 0;
      for k = 1:nr_sat,
         if  elangle(k) >= ele_lim
            nr_vis(kk) = nr_vis(kk) + 1;
            ind_vis(nr_vis(kk),kk) = almanac(k,1);
            el_vis(nr_vis(kk),kk) = elangle(k);
            az_vis(nr_vis(kk),kk) = azangle(k);
         end
      end

%     For each geographic location save the following computed quantities into 
%     an external file: kkk, kk, nr_vis, ind_vis, el_vis, az_vis

      for i = 1:nr_vis(kk)
         fprintf(pf1,'%5.0f %5.0f %5.0f %5.0f %18.7f %18.7f\n',...
                 kkk,kk,nr_vis(kk),ind_vis(i,kk),el_vis(i,kk),az_vis(i,kk));
      end

   end                       %  end of the location loop (index kk)

   tsim = tsim + t_incr;

end                          %  end of the time loop (index kkk)

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