📄 xmapdop.m
字号:
% xmapdop.m
% Scope: This Matlab program determines repartition of the number of visible
% satellites and the corresponding dops for a specified geographical
% area.
% Usage: xmapdop
% 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
% - geographic locations for the test; the default is data file
% map01.dat. The first record of the file should contain the latitude
% limits while the second record should contain the longitude limits,
% both in degrees
% - elevation angle limit; the default value is 5 degrees
% - simulation time (GPS time of week), in seconds; the default value
% is 0.
% - selection of xdop to be plotted
% - number of time samples and time step value, in seconds (if more
% than one time step is selected)
% - name of the output file for satellite visibility statistics (optional)
% - other parameters are initialized by default; the default latitude-
% longitude grid is 25 by 41 (default)
% Outputs: - plot presenting selected dop value and/or number of visible
% satellites versus latitude-longitude location
% - table containing test summary and satellite visibility
% statistics (optional)
% Remarks: If more than one step then only the results from the last time sample
% are plotted.
% External Matlab macros used: dop1, eleva, svpalm, tgdecef, uverv, vecefenu,
% wgs84con, and world map data file worldmap.mat (optional)
% Last update: 06/26/00
% Copyright (C) 1998-00 by LL Consulting. All Rights Reserved.
clear
close all
yes = 'y';
deg_to_rad = pi / 180.; % degrees to radians transformation
% Read the input file containing satellites data
disp(' ');
disp('Enter almanac data file - 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 --> ','s');
disp(' ');
end
almanac = load(ff);
nr_sat = size(almanac,1); % number of satellites used
% Read the geographic location for the test
disp('Enter geographic location data - the default is the data file map01.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)
hh = load('map01.dat');
ylat(1) = hh(1,1);
ylat(2) = hh(1,2);
xlon(1) = hh(2,1);
xlon(2) = hh(2,2);
else
ylat(1) = input('Enter the south most latitude desired, in degrees --> ');
ylat(2) = input('Enter the north most latitude desired, in degrees --> ');
xlon(1) = input('Enter the west most longitude desired, in degrees --> ');
xlon(2) = input('Enter the east most longitude desired, in degrees --> ');
end
ylat_lim = ylat * deg_to_rad; % in radians
xlon_lim = xlon * deg_to_rad; % in radians
% Determine the longitude-latitude grid (default is 25 by 41)
ylat_nr = 25; % default value (can be changed)
xlon_nr = 41; % default value (can be changed)
ylat_incr = (ylat_lim(2) - ylat_lim(1))/(ylat_nr - 1);
xlon_incr = (xlon_lim(2) - xlon_lim(1))/(xlon_nr - 1);
ylatgrid(1:ylat_nr,1) = (ylat_lim(1):ylat_incr:ylat_lim(2))';
xlongrid(1:xlon_nr,1) = (xlon_lim(1):xlon_incr:xlon_lim(2))';
user = zeros(3,1);
% Compute the geographic locations in ECEF
for kylat = 1:ylat_nr
for kxlon = 1:xlon_nr
temp = tgdecef(ylatgrid(kylat),xlongrid(kxlon),0.);
uecefx(kxlon,kylat) = temp(1);
uecefy(kxlon,kylat) = temp(2);
uecefz(kxlon,kylat) = temp(3);
end
end
% 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 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 simulation time (GPS time of week) --> ');
end
tsim = tweek;
% Specify what xdop is plotted
disp(' ');
disp('Select xdop to be plotted');
disp('Enter 1 for gdop');
disp('Enter 2 for pdop');
disp('Enter 3 for hdop');
disp('Enter 4 for vdop');
disp('Enter 5 for tdop');
ind_dop = input('Make selection --> ');
if ind_dop == 1
xfirst = 'G';
elseif ind_dop == 2
xfirst = 'P';
elseif ind_dop == 3
xfirst = 'H';
elseif ind_dop == 4
xfirst = 'V';
else
xfirst = 'T';
end
% Select number of time samples and time step
disp(' ');
disp('Enter number of time samples - the default value is 1 ');
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 = 1; % 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, in seconds - the default value is 3600. ');
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 = 3600.; % default - time step
else
t_incr = input('Enter time step, in seconds --> ' );
end
end
% Specify the output file for visibility statistics
disp(' ');
answer7 = input('Do you want to save the visibility statistics? (y/n)[y] ','s');
if isempty(answer7)
answer7 = yes;
end
if (strcmp(answer7,yes) == 1)
pf1 = input('Specify the filename for satellite visibility statistics --> ','s');
disp(' ');
end
% Initialization of counters for number of visible satellites
if (strcmp(answer7,yes) == 1)
count4 = 0.;
count5 = 0.;
count6 = 0.;
count7 = 0.;
count8 = 0.;
count9 = 0.;
count10 = 0.;
count11 = 0.;
count12 = 0.;
end
% Start the main computation loop - time loop
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
% Compute satellite position in ECEF
temp = almanac(k,2:9);
psat_temp = (svpalm(tsim,temp))';
psat(1:3,k) = psat_temp;
end % end of the satellite position loop (index k)
for kylat = 1:ylat_nr % cycle for all locations (latitude)
for kxlon = 1:xlon_nr % cycle for all locations (longitude)
user(1) = uecefx(kxlon,kylat);
user(2) = uecefy(kxlon,kylat);
user(3) = uecefz(kxlon,kylat);
for k = 1:nr_sat, % cycle for all satellites
% Compute elevation angle
temp2 = psat(:,k);
[temp3,temp1] = eleva(user,temp2);
eangle(k) = temp3;
ulos(:,k) = vecefenu(temp1,ylatgrid(kylat),xlongrid(kxlon));
end % end of the satellite loop (index k)
% Determine the number of visible satellites and save the index and
% elevation angle for all visible satellites for a specific location
clear temp
clear htransp
nr_vis(kxlon,kylat) = 0;
for k = 1:nr_sat,
if eangle(k) >= ele_lim
nr_vis(kxlon,kylat) = nr_vis(kxlon,kylat) + 1;
temp4 = eangle(k);
if nr_vis(kxlon,kylat) == 1
htransp = ulos(:,k);
else
temp = ulos(:,k);
htransp = [htransp temp];
end
end
end
clear hhtt
hhtt = htransp';
% Compute and save the selected xdop
dops = dop1(hhtt);
xdop(kxlon,kylat) = dops(ind_dop);
% Determine the number of cases with less than 5, 5, 6, 7, 8, 9, 10, 11,
% and greater than 11 visible satellites
if (strcmp(answer7,yes) == 1)
nrtemp = nr_vis(kxlon,kylat);
if nrtemp <= 4
count4 = count4 + 1;
elseif nrtemp == 5
count5 = count5 + 1;
elseif nrtemp == 6
count6 = count6 + 1;
elseif nrtemp == 7
count7 = count7 + 1;
elseif nrtemp == 8
count8 = count8 + 1;
elseif nrtemp == 9
count9 = count9 + 1;
elseif nrtemp == 10
count10 = count10 + 1;
elseif nrtemp == 11
count11 = count11 + 1;
else
count12 = count12 + 1;
end
end
end % end of the location loops (index kxlon)
end % end of the location loops (index kylat)
tsim = tsim + t_incr;
end % end of the time loop (index kkk)
% Execute xdop and visible satellite graphs
xlongrid = xlongrid / deg_to_rad;
ylatgrid = ylatgrid / deg_to_rad;
figure(1)
cc = contour(xlongrid,ylatgrid,xdop');...
xlabel('Longitude, in degrees'),...
ylabel('Latitude, in degrees'),...
twk = num2str(tsim);
title([xfirst 'DOP versus longitude-latitude grid, for ' ff ' at time-of-week = ' twk]);
clabel(cc);
disp(' ');
answer8 = input('Include the geographic map? (y/n)[y] --> ','s');
if isempty(answer8)
answer8 = yes;
end
if (strcmp(answer8,yes) == 1)
hold on
westlon = xlon(1);
eastlon = xlon(2);
lowlat = ylat(1);
highlat = ylat(2);
load worldmap
k = 1;
[nrow,ncol] = size(long);
for i = 1:nrow
if (lat(i) < highlat) & (lat(i) > lowlat)
if (long(i) > westlon) & (long(i) < eastlon)
maplat(k) = lat(i);
maplon(k) = long(i);
k = k + 1;
end
end
end
plot(maplon,maplat,'.');
end
disp(' ');
answer9 = input('Plot the visible satellite contour? (y/n)[y] --> ','s');
if isempty(answer9)
answer9 = yes;
end
disp(' ');
if (strcmp(answer9,yes) == 1)
% Graph number of visible satellites
figure(2)
v = [4 5 6 7 8 9 10 11 12]';
ccc = contour(xlongrid,ylatgrid,nr_vis',v);
% Note: The array nr_vis should have different elements!
xtitle = 'Number of visible satellites versus longitude-latitude grid, for ';
title([xtitle ff ' at time-of-week = ' twk]);
xlabel('Longitude, in degrees'),...
ylabel('Latitude, in degrees'),...
clabel(ccc);
answer10 = input('Include the geographic map? (y/n)[y] --> ','s');
if isempty(answer10)
answer10 = yes;
end
if (strcmp(answer10,yes) == 1)
hold on
westlon = xlon(1);
eastlon = xlon(2);
lowlat = ylat(1);
highlat = ylat(2);
load worldmap
k = 1;
[nrow,ncol] = size(long);
for i = 1:nrow
if (lat(i) < highlat) & (lat(i) > lowlat)
if (long(i) > westlon) & (long(i) < eastlon)
maplat(k) = lat(i);
maplon(k) = long(i);
k = k + 1;
end
end
end
plot(maplon,maplat,'.');
end
end
% Save satellite visibility statistics (optional)
if (strcmp(answer7,yes) == 1)
nr_sp_time = ylat_nr * xlon_nr * t_sample;
fprintf(pf1,'***********************************************************************\n\n');
fprintf(pf1,' XMAPDOP - RESULTS \n\n');
fprintf(pf1,'************************************************************************\n');
fprintf(pf1,'Almanac: %s\n',ff);
fprintf(pf1,'Time of week: %9.2f\n', tweek);
fprintf(pf1,'Geographic area:\n');
fprintf(pf1,'- longitude (in degrees): %6.2f to %6.2f\n', xlon(1), xlon(2));
fprintf(pf1,'- latitude (in degrees): %6.2f to %6.2f\n', ylat(1), ylat(2)); fprintf(pf1,'Number of space-time points analyzed = %6.0f\n',nr_sp_time);
fprintf(pf1,'***********************************************************************\n');
fprintf(pf1,'\n SUMMARY - SATELLITE VISIBILITY STATISTICS \n');
fprintf(pf1,'\nNumber of cases with <5 visible satellites = %4.0f, percentage = %6.2f\n',count4, (100.*count4/nr_sp_time) );
fprintf(pf1,'Number of cases with 5 visible satellites = %4.0f, percentage = %6.2f\n',count5, (100.*count5/nr_sp_time));
fprintf(pf1,'Number of cases with 6 visible satellites = %4.0f, percentage = %6.2f\n',count6, (100.*count6/nr_sp_time));
fprintf(pf1,'Number of cases with 7 visible satellites = %4.0f, percentage = %6.2f\n',count7, (100.*count7/nr_sp_time));
fprintf(pf1,'Number of cases with 8 visible satellites = %4.0f, percentage = %6.2f\n',count8, (100.*count8/nr_sp_time));
fprintf(pf1,'Number of cases with 9 visible satellites = %4.0f, percentage = %6.2f\n',count9, (100.*count9/nr_sp_time));
fprintf(pf1,'Number of cases with 10 visible satellites = %4.0f, percentage = %6.2f\n',count10, (100.*count10/nr_sp_time));
fprintf(pf1,'Number of cases with 11 visible satellites = %4.0f, percentage = %6.2f\n',count11, (100.*count11/nr_sp_time));
fprintf(pf1,'Number of cases with >11 visible satellites = %4.0f, percentage = %6.2f\n\n',count12, (100.*count12/nr_sp_time));
fprintf(pf1,'***********************************************************************\n');
end
disp(' ');
disp('End of the program XMAPDOP ');
disp(' ');
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -