xuercor.m

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

%                                xuercor.m
%  Scope:   This MATLAB program generates and plots the magnitude of the user earth 
%           rotation correction vector (Coriolis factor) for a specified longitude/
%           latitude grid. WGS-84 constants are used.
%  Usage:   xuercor
%  Inputs:  - name of the input almanac data file; the default is data file 
%             wk749.dat. Each record contains the 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
%           - elevation mask angle, in degrees; the default value is 5 degrees
%           - simulation time, time of week, in seconds; the default value is 0.
%           - index of the selected satellite (prn)
%           - longitude/latitude/altitude grid by default (longitude and latitude by
%             10 degrees from 0 to 360, and -90 to 90, respectively, and altitude is
%             always 0 meters); the grid can be easily modified
%           - option to save the generated data in the file ercor.mat
%  Outputs: - generated earth rotation correction vector magnitude for the
%             longitude/latitude/altitude grid into the file ercor.mat (optional)
%           - surface plot of the magnitude of the earth rotation correction vector
%             magnitude for the longitude/latitude/altitude grid 
%  External Matlab macros used:  eleva, svpalm, tgdecef, uercor, uverv, wgs84con     
%  References:  [1] Anon., Navstar GPS Space Segment/Navigation User Interfaces, 
%                   ICD-GPS-200C, October 10, 1993
%               [2] DiEsposti, R., Time-dependency and coordinate system issue in GPS
%                   measurement models. ION GPS-2000, Salt Lake City, UT, September
%                   19-22, 2000
%  Last update:  09/29/00
%  Copyright (C) 1998-00 by LL Consulting. All Rights Reserved.

clear
yes = 'y';
close all

%  Constants
deg_to_rad = pi / 180.; %  degrees to radians transformation
wgs84con;               %  global variable used:  c_speed

%  Read the input file containing almanac data

disp('  ');
disp('Enter almanac data file - the default is datafile wk749.dat');
answer1 = input('Do you want to use the default datafile? (y/n)[y] ','s');
if isempty(answer1)
   answer1 = yes;
end   
if (strcmp(answer1,yes) == 1)
   ff = 'wk749.dat';
else
   ff = input('Specify the input filename (with extension) -->  ','s');
end
almanac = load(ff);
[nsat,ncol] = size(almanac);         %  number of satellites used
if  ncol ~= 9
   disp('Error - XUERCOR; check the almanac input file');
   disp('  ');
   return
end   

%  Specify the elevation mask angle limit

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

%  Select simulation time (GPS time of week)

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

% Select the satellite index

disp('  ');
disp('Satellites in the almanac:');
for k = 1:nsat  
   fprintf('%3.0f ',almanac(k,1));
end
disp('  ');
prn = input('Select satellite index (prn) --> ');
for k = 1:nsat
   if  (almanac(k,1) == prn)
      satind = k;
      break
   end   
end   

%  Select the grid from 10 to 10 degrees for latitude/longitude, and altitude = 0.

latdeg = -90:10:90;                                  %  latitude in degrees                
londeg = 0:10:350;                                   %  longitude in degrees
alt = 0.;                                            %  altitude above ellipsoid, in meters
latrad   = latdeg * deg_to_rad;                      %  latitude in radians
lonrad   = londeg * deg_to_rad;                      %  longitude in radians

nlat = size(latdeg,2);
nlon = size(londeg,2);
npt = nlat*nlon;                                     %  total number of points in the grid
user_ecef = zeros(3,npt);
xcor = zeros(nlon,nlat);

%  Compute the user positions in ECEF

kk = 0;
for klat = 1:nlat
   for klon = 1:nlon
      kk = kk + 1;
      temp = tgdecef(latrad(klat),lonrad(klon),alt);
      user_ecef(:,kk) = temp; 
   end   
end
      
%  Determine selected satellite position in ECEF

temp = almanac(satind,2:9);
psat = (svpalm(tsim,temp))';

%  Main loop for earth rotation correction vector computation

for kk = 1:npt          %  cycle for all locations
   
   klon = 1 + mod(kk-1,nlon);
   klat = 1 + fix((kk-1)/nlon);
   upos = user_ecef(:,kk);
   
%  Compute elevation angle

   [elangle,temp1] = eleva(upos,psat);
                
%  Determine the earth rotation correction vector magnitude when the satellite is visible

   if  elangle >= ele_lim
      dt = norm(upos-psat) / c_speed;
      xercor = uercor(upos,dt);
      xcor(klon,klat) = norm(xercor);
   else
      xcor(klon,klat) = 0.;
   end

end                    %  end of the location loop (index kk)

%  Save the generated data into a mat-file (optional)

disp('  ');
answer4 = input('Do you to save the generated data into mat-file? (y/n)[n] ','s');
if isempty(answer4)
   answer4 = 'n';
end   
if (strcmp(answer4,yes) == 1)
   disp('  ');
   disp('The magnitude of the Earth Rotation Correction Vector is saved in ercor.mat');
   save   ercor  xcor
end

%  Execute the graph of the magnitude of the earth rotation correction vector

[xx,yy] = meshgrid(latdeg,londeg);
aa = 'User Earth Rotation Correction Vector Magnitude (meters), for satellite ';
figure(1)
surf(xx,yy,xcor), ...
axis([-90 90 0 360 0 50]),...   
title([aa num2str(prn)]),...
ylabel('Longitude, in degrees'),...
xlabel('Latitude, in degrees'),...
zlabel('Earth Rotation Correction Vector Mganitude, in meters'),...
colormap([1 1 0; 0 1 1])

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