xionoc.m

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

%                                xionoc.m
%  Scope:  This Matlab program determines iono corrections for a specified 
%          satellite, time interval, and a selected location by using the
%          Klobuchar model; several representative plots are available.
%          WGS-84 constants are used.
%  Usage:  xionoc 
%  Inputs:   - name of the input almanac data file; the default is data file 
%              wk749.dat. Each record contains the data related to a specific
%              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,
%            - name of the output file (optional); the default name is xionoc1.out
%            - other parameters are initialized by default
%  Outputs:  - iono correction data into a specified file (optional); each record
%              contains iono correction data for all satellites at a time sample 
%              (the value 0. is recorded when the satellite is not visible)
%            - plot of number of visible satellites versus time sample number,
%              when time sample number is greater than 1
%            - 3-D plot of iono correction versus time sample number and 
%              satellite id, when time sample number is greater than 1
%            - plot of iono correction versus satellite id for a selected time
%              sample number 
%            - plot of iono correction versus time number for a selected 
%              satellite, when time sample number is greater than 1
%  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: eleva, ionoc, ionocon, svpalm, tecefgd, 
%                               tgdecef, uverv, vecefenu, wgs84con, 
%  Last update:  12/23/00
%  Copyright (C) 1999-00 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';

%  Initialize the iono constants alpha and beta

ionocon

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

%  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

%  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_to_rad;         %  in radians
lon_rad   = lon * deg_to_rad;         %  in radians

%  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_to_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
disp('  ');

%  Start main computation loop - time loop

nr_vis = zeros(1,t_sample);
t_iono = zeros(nr_sat,t_sample);
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 iono correction
%  for each visible satellite; if the satellite is not visible the iono 
%  correction is set to 0.

   for k = 1:nr_sat,
      if  eangle(k) >= ele_lim
         nr_vis(kkk) = nr_vis(kkk) + 1;
         vis_sat(k,kkk) = k;
         el = eangle(k);
         az = aangle(k); 
         ionocorr = ionoc(lat_rad,lon_rad,el,az,tsim,alpha,beta);
         t_iono(k,kkk) = ionocorr;
      else
         vis_sat(k,kkk) = 0;
         t_iono(k,kkk) = 0.;
      end
   end

   tsim = tsim + t_incr;

end      %  end of the time loop (index kkk)

%  Save the iono correction data (optional)

disp('  ');
answer7 = input('Do you want to save the iono correction data? (y/n)[n] ','s');
if isempty(answer7)
   answer7 = 'n';
end
if (strcmp(answer7,yes) == 1)
   disp('  ');
   answer8 = input('Do you want to use the default file, xionoc1.out? (y/n)[y] ','s');
   if isempty(answer8)
      answer8 = yes;
   end
   if (strcmp(answer8,yes) == 1)
      fout = 'xionoc1.out';         %  default file name 
   else
      disp('  ');
      fout = input('Enter the name of the output file (with extension) -->  ','s'); 
   end
   for k = 1:t_sample
      for k = 1:nr_sat
         fprintf(fout,'%10.5f ',t_iono(k,kkk));     
      end
      fprintf(fout,'\n');
   end
end

%  Execute the plots

time = 1:t_sample;
ind_sv_max = max(ind_sv);
t_iono_sv = zeros(ind_sv_max,t_sample);
ind_sv_sv = zeros(1,ind_sv_max);
for k = 1:nr_sat
   t_iono_sv(ind_sv(k),:) = t_iono(k,:);
   ind_sv_sv(ind_sv(k)) = ind_sv(k);
end   

%  Number of visible satellites versus time sample

if (t_sample > 1)
   figure(1)       
   plot(time,nr_vis,'*')
   grid, axis([1 t_sample 3 12]);
   aa =['Time sample number (time step = ' num2str(t_incr) ' seconds, '];
   aa = [aa 'start time tow = ' num2str(tweek) ' seconds)'];
   xlabel(aa);
   ylabel('Number of visible satellites');
   cc = 'Number of visible satellites  versus  Time sample number  (';
   cc = [cc  'almanac  ' ff ')'];
   title(cc);
else
   fprintf('\nNumber of visible satellites = %3.0f\n',nr_vis(1));
end   

%  Iono correction versus time sample number and satellite id

if (t_sample > 1)
   figure(2)
   surf(t_iono_sv);
   axis([0 t_sample 0 ind_sv_max 0 16.]);
   xlabel('Time sample number');
   ylabel('Satellite number');
   zlabel('Iono correction, in meters');
   title('Iono correction versus Time sample and Satellite number');
else
   for k = 1:nr_sat
      if  (t_iono(k,1) ~= 0.0)
         fprintf('\nIono correction for satellite %3.0f is %11.4f meters\n',ind_sv(k),t_iono(k,1));
      end  
   end 
end   

%  Iono correction versus satellite id for a selected time sample

fprintf('\nNumber of time sample available = %3.0f\n',t_sample);
sel_sample = input('Select the time sample number to be analyzed --> ');
figure(3)
plot(ind_sv,t_iono(:,sel_sample),'*');
grid, axis([1 ind_sv_max -1 16]);
xlabel(['Satellite number  (almanac ' ff ')']);
ylabel('Iono correction, in meters  (0 if satellite is not visible)');
title(['Iono correction versus Satellite number, for Time sample # ' num2str(sel_sample)]);

%  Iono correction versus time sample for a selected satellite

fprintf('\nThe following satellites are in the almanac:\n');
for k = 1:nr_sat
   fprintf('%3.0f',ind_sv(k));
end  
disp(' ');
sel_sat = input('Select the satellite id --> ');
if (t_sample > 1)
   figure(4)
   plot(time,t_iono_sv(sel_sat,:));
   hold on
   plot(time,t_iono_sv(sel_sat,:),'*');
   grid
   xlabel(aa);
   ylabel('Iono correction, in meters');
   cc = ['Iono correction  versus  Time sample number, for satellte # ' num2str(sel_sat)];
   cc = [cc  ' (almanac  ' ff ')'];
   title(cc);
else
   fprintf('\nIono correction for satellite %3.0f is %11.4f meters\n',sel_sat,t_iono_sv(sel_sat,1));
end   
disp('  ');
disp('End of the program  XIONOC');
disp('  ');