plotpass.m

此功能包用于各种GPS坐标和时间的转换

function p_handle = plotpass(t,plot_data,prn_pass,plot_title,y_label,y_scale)

% p_handle = plotpass(t,plot_data,prn_pass,plot_title,y_label,y_scale);
%
% Function to plot y data as a function of time for a single observer
% and return the plot handles.  Each arc generated by a satellite pass is colored
% and labeled with the prn-pass number. If all the prn numbers are 0, then no prn
% numbers will be added to the plot. New figures are created to plot all 
% plot_data. Zoom capability is added to the figure. See help on ZOOM for more
% information.
%
% Input:
%   t          - time in GPS format (nx2) [GPS_week GPS_sec]
%                GPS_week values may range from 0 to 3640
%                GPS_sec values may range from 0 to 604800
%   plot_data  - data to be plotted (nxm). each column is a different
%                data set, i.e. [azimuth, elevation, range]
%   prn_pass   - [satellite number, pass number], nx1 or nx2 (optional),
%                pass numbers are computed by PASSDATA. 
%                if not provided, a single line is drawn
%   plot_title - title of plot (1xj string or mxj string) (optional) 
%                default = 'Data Plot'
%   y_label    - label for the y-axis, (1xk string or mxk string) (optional) 
%                default = 'Y Data'
%   y_scale    - min and max limits of y scale (1x2 or mx2) (optional). 
%                units match plot_data.
%                if provided, data will be truncated when passing from
%                one limit to the other. Example: Azimuth data that goes
%                from 360 to 0 degrees during a pass will not have a vertical
%                line connecting those 2 points.
%   Note: if only 1 plot_title, y_label, or y_scale is provided for multiple
%         sets of plot_data, then it will be used for all plots.
% Output:
%   p_handle   - graphics handles to the figures (mx1)
%
% See also PLOTSKY

% Written by: Maria Evans 5/4/98
% Copyright (c) 1998 by Constell, Inc.

% functions called: ERR_CHK, GPS2UTC, UTC2GPS

%%%%% BEGIN VARIABLE CHECKING CODE %%%%%
% declare the global debug variable
global DEBUG_MODE

% Initialize the output variables
p_handle=[];

% Check the number of input arguments and issues a message if invalid
msg = nargchk(2,6,nargin);
if ~isempty(msg)
  fprintf('%s  See help on PLOTPASS for details.\n',msg);
  fprintf('Returning with empty outputs.\n\n');
  return
end

if nargin < 3,
  prn_pass = zeros(size(plot_data,1),2);
end % if nargin < 3

% Check to make sure the title string is a valid size
if nargin < 4
  plot_title = ['Data Plot'];
end % if nargin < 4

% Check to make sure the y-label string is a valid size
if nargin < 5
  y_label = ['Y Data'];
end % if nargin < 5

% Get the current Matlab version
matlab_version = version;
matlab_version = str2num(matlab_version(1));

% If the Matlab version is 5.x and the DEBUG_MODE flag is not set
% then set up the error checking structure and call the error routine.
if matlab_version >= 5.0                        
  estruct.func_name = 'PLOTPASS';

  % Develop the error checking structure with required dimension, matching
  % dimension flags, and input dimensions.
  estruct.variable(1).name = 't';
  estruct.variable(1).req_dim = [901 2];
  estruct.variable(1).var = t;
  
  estruct.variable(2).name = 'plot_data';
  estruct.variable(2).req_dim = [901 1];
  estruct.variable(2).var = plot_data(:,1);
  
  estruct.variable(3).name = 'prn_pass';
  estruct.variable(3).req_dim = [901 1; 901 2;];
  estruct.variable(3).var = prn_pass;

  if nargin==6,
    estruct.variable(4).name = 'y_scale';
    estruct.variable(4).req_dim = [1 2; size(plot_data,2) 2;];
    estruct.variable(4).var = y_scale;
  end;
  
  % Call the error checking function
  stop_flag = err_chk(estruct);
  
  if stop_flag == 1           
    fprintf('Invalid inputs to %s.  Returning with empty outputs.\n\n', ...
             estruct.func_name);
    return
  end % if stop_flag == 1
end % if matlab_version >= 5.0 & isempty(DEBUG_MODE) 

%%%%% END VARIABLE CHECKING CODE %%%%%

%%%%% BEGIN ALGORITHM CODE %%%%%

% convert GPS time to linear time past an epoch
week_min = min(t(:,1));

% seconds past a minimum week
t_lin = (t(:,1) - week_min) * 604800 + t(:,2);

% sort the data in time (linear time)
[t_sort, I_sort] = sort(t_lin);

t_lin = t_lin(I_sort);
t = t(I_sort,:);
plot_data = plot_data(I_sort,:);
prn = prn_pass(I_sort,1);
if size(prn_pass,2)==2,
  pass = prn_pass(I_sort,2);
else,
  pass = ones(size(prn));
end;

% compute the starting UTC time tag
start_utc = gps2utc(t(1,:));
  
% total duration (seconds)
duration = max(t_lin) - min(t_lin);

% set the time scale to be used for the plotting (minutes, hours, or days)
if duration <= 3600        % shorter than an hour
  min_flag = 1;            % flag to use the minute time scale
  hour_flag = 0;           % flag not to use the hour time scale
  day_flag = 0;            % flag not to use the day time scale 
  
  % round to the start of a second
  start_utc(6) = 0;        % 0 seconds into the minute
  
  % convert to GPS time
  start_time_gps = utc2gps(start_utc);        
  
  t_lin_start = (start_time_gps(:,1) - week_min) * 604800 + start_time_gps(:,2);
  
  t_plot = (t_lin - t_lin_start) / 60;     % plotting times in minutes
  
elseif duration > 3600 & duration <= 86400  % between an hour and a day
  min_flag = 0;            % flag not to use the minute time scale
  hour_flag = 1;           % flag to use the hour time scale
  day_flag = 0;            % flag not to use the day time scale

  % round to the start of a minute
  start_utc(5) = 0;        % 0 minutes into the hour
  start_utc(6) = 0;        % 0 seconds into the minute
  
  % convert to GPS time
  start_time_gps = utc2gps(start_utc);
  
  t_lin_start = (start_time_gps(:,1) - week_min) * 604800 + start_time_gps(:,2);

  t_plot = (t_lin - t_lin_start) / 3600;   % plotting times in hours

elseif duration > 86400    % longer than a day
  min_flag = 0;            % flag not to use the minute time scale
  hour_flag = 0;           % flag not to use the hour time scale
  day_flag = 1;            % flag to use the day time scale

  % round to the beginning of the day
  start_utc(4) = 0;        % 0 hours into the day
  start_utc(5) = 0;        % 0 minutes into the day
  start_utc(6) = 0;        % 0 seconds into the day
  
  % convert to GPS time
  start_time_gps = utc2gps(start_utc);
  
  t_lin_start = (start_time_gps(:,1) - week_min) * 604800 + start_time_gps(:,2);

  t_plot = (t_lin - t_lin_start) / 86400;  % plotting times in days

end % if duration <= 3600 

% Create the x axis label
year_2_digit = start_utc(1)-1900;
if (year_2_digit >= 100)
   year_2_digit = 0;
end % if

start_time_string = sprintf('%d/%d/%02d %d:%02d:%02d',  ...
                    start_utc(2),start_utc(3),year_2_digit,start_utc(4:6));

if min_flag == 1        % using the minute time scale
  xlabel_text = sprintf('Minutes past %s ',start_time_string);

elseif hour_flag == 1   % using the hour time scale
  xlabel_text = sprintf('Hours past %s ',start_time_string);

elseif day_flag == 1    % using the day time scale
  xlabel_text = sprintf('Days past %s ',start_time_string);

end % if min_flag == 1    

% determine the total number of satellites
% start by sorting the prn numbers and looking for changes
prn_sort = sort(prn);
prn_change = find(diff(prn_sort) ~= 0);

% create a matrix that has sorted and reduced prns [1 2 4 5 6 8 ... 28]
active_sats = [prn_sort(prn_change); prn_sort(length(prn_sort))];

% compute the total number of prns 
num_prns = length(active_sats);

% Set to black background
%colordef none

% set the colors that will be used in this function
avail_colors = ['ymcrgbw'];

num_colors = length(avail_colors);

% Determine number of plots to be created
num_plots = size(plot_data,2);

% get the current figure handles
max_wins = max(get(0,'children'));
if isempty(max_wins),
  max_wins = 0;
end;

% Get the screen size in pixels to use for location of plots
set(0,'units','pixels');
screen_size = get(0,'screensize');
y_max = screen_size(2) + screen_size(4) - 60;
x_max = screen_size(1) + screen_size(3) - 70;
x_step = 80;

% Determine the location for the first azimuth plot in upper right corner
x_min = x_max/2;
y_min = 30;

% Loop over necessary number of plots to be created
for mmm=1:num_plots,

  % check to see if an individual title was provided for this plot
  if size(plot_title,1) > 1,
    title_fig = plot_title(mmm,:);
  else,
    title_fig = plot_title(1,:);
  end;

  p_handle(mmm) = figure(max_wins+mmm);
  set(p_handle(mmm),'position',[x_min-x_step*(mmm-1) y_min x_max/2 y_max/2], ...
        'NumberTitle','off','Units','Pixels', ...
        'Name',sprintf('%s',deblank(title_fig)));
  hold on

  % check to see if an individual y_label was provided for this plot
  if size(y_label,1) > 1,
    y_label_fig = y_label(mmm,:);
  else,
    y_label_fig = y_label(1,:);
  end;

  for i = 1:num_prns               % loop over the active satellites
    J = find(prn == active_sats(i));     % all satellites with the same prn number 
    % Determine how many passes for this PRN
    pass_sort = sort(pass(J));
    pass_change = find(diff(pass_sort) ~= 0);

    % create a matrix that has sorted and reduced pass numbers [1 2 4 5 6 8 ... 28]
    pass_nums = [pass_sort(pass_change); pass_sort(length(pass_sort))];
    num_passes = length(pass_nums);

    if num_prns > 1, % For more than 1 prn, change colors each PRN
      label_color = avail_colors(rem(i,num_colors) + 1);  % cycle through colors
    end;

    for k = 1:num_passes,
      if num_prns == 1, % For 1 prn, change colors each pass
        label_color = avail_colors(rem(k,num_colors) + 1);  % cycle through colors
      end;

      I = find(prn == active_sats(i) & pass == pass_nums(k));   

      if nargin == 6,
        % check for a wrap in data (i.e. azimuth)
        daz = diff(plot_data(I,mmm));
        % find all of the data changes > (y_scale(2)-y_scale(1))/1.5 (a wrap)
        if size(y_scale,1) > 1,
          I_daz = find(abs(daz) > (y_scale(mmm,2)-y_scale(mmm,1))/1.5); 
        else,
          I_daz = find(abs(daz) > (y_scale(1,2)-y_scale(1,1))/1.5); 
        end;
  
        if any(I_daz)
          % make this generic to handle multiple wraps 
          % set the starting and ending points for the first wrap plot
          first = 1;
          last = I_daz(1);
    
          % loop over the number of wraps
          for j = 1:length(I_daz) + 1                                          
            plot(t_plot(I(first:last)),plot_data(I(first:last),mmm), ...
                  'Color',label_color);
  
            % change the end points for the next pass through the loop
            first = last + 1;    
       
            if j >= length(I_daz)
              last = length(I);
            else  
              last = I_daz(j+1);
            end
          end
        else,
          % plot the az/az pairs for this satellite
          plot(t_plot(I),plot_data(I,mmm),'Color',label_color);
        end;  % if any(daz)
 
        hold on                  % turn hold on so the next satellite is added
      else,
        % plot the az/az pairs for this satellite
        plot(t_plot(I),plot_data(I,mmm),'Color',label_color);
      end;

      % add prn label to each data arc
      % put the label in the center if the satellite is visible then, otherwise put the
      % label at the beginning or the end
      if length(t_plot(I)) < 3  % not even 2 points, no labels
        I_mid = []; 
        I_end = [];  
        I_beg = []; 
  
      else
        I_mid = fix(length(I)/2);      % mid point
        I_end = length(I);             % end point
        I_beg = 1;                     % start point
  
      end % if size(t_plot(I) < 3)  % not even 2 points, no labels

      if active_sats(i) ~= 0,  
        if any(I_mid)
          % put the label on
          qts = text(t_plot(I(I_mid)),plot_data(I(I_mid),mmm), ...
                  num2str(active_sats(i)),'Color',label_color);     
          clear I_mid
        elseif any(I_end)
          % put the label on
          qts = text(t_plot(I(I_end)),plot_data(I(I_end),mmm), ...
                  num2str(active_sats(i)),'Color',label_color);         
          clear I_end
        elseif any(I_beg)
          % put the label on
          qts = text(t_plot(I(I_beg)),plot_data(I(I_beg),mmm), ...
                 num2str(active_sats(i)),'Color',label_color);         
          clear I_beg
        end % if any(I_mid)
      end %if if active_sats(i) ~= 0, 
    end;  % for k=1:num_passes,
  end % for i=1:num_prns,


  % label the plot
  ylabel(y_label_fig)

  if nargin==6,
    v = axis;
    if size(y_scale,1) > 1,
      axis([v(1) v(2) y_scale(mmm,1) y_scale(mmm,2)]);
    else,
      axis([v(1) v(2) y_scale(1,1) y_scale(1,2)]);
    end;
  end;

  % set the axis box on
  set(gca,'box','on');

  % set the plot label to be title_fig
  set(gcf,'Name',title_fig)

  % set the title
  title(title_fig)
  zoom on       

  % finally label the x axis
  xlabel(xlabel_text) 

end % for mmm=1:num_plots,

%%%%% END ALGORITHM CODE %%%%%

% end PLOTPASS