function draw(r, bmin, bmax, pmin, pmax)
%
% Usage:  draw(r)
%
% Plots the columns of array r, displaced by global array DISP,
% using '+' for the "key" profile (given by IKEY) and lines for the rest.
%

% 96/10/02 JF: modified to never plot bins with percent_good less than 30
% 97/05/15 JF: changed display message when showing smoothing and
%  truncation effects of set_lgb (formerly botmpas3 and last_85) to avoid
%  confusion of old and new bottom editing commands; matlab routines called
%  here retain the old names
% 97/06/19 JF: global variables BOTMPAS3 & LAST_85 removed from editglob.m
%  list, BOTMEDIT added; if statement now checks BOTMEDIT

global PMIN PMAX BMIN BMAX IKEY DISP
global PFLAG PLOTBAD PLOTBAD_ORIGINAL  NBINS NPRFS
global DAYS IBLKPRF DBNAME LON LAT U V W AMP PGOOD E CURRX U_ROTATE V_ROTATE
global BADACC BADBIN BADPRF BOTTOM BOTMPAS3 LAST_85 BOTM_EDIT
global MAB YAXIS
global GOOD_PG

%editglob

% plotting symbols:

wvar = '*m';           % high W var = magenta *
d2   = 'ow';           % d2 failure = white o
ev   = 'oy';           % high E vel = yellow o
db   = 'xw';           % already flagged in database = white x
amp  = '*c';           % max amp bin = cyan *
old_amp = 'xg';        % max amp bin already in database = blue x
user_added = '+y';     % user-added editing = yellow +
range_pg = 'xr';       % flagged via range criterion and/or pg

if nargin < 5
  pmax = PMAX;
  if nargin < 4
    pmin = PMIN;
    if nargin < 3
      bmax = BMAX;
      if nargin < 2
        bmin = BMIN;
      end
    end
  end
end

%% The following block compares the input matrix being plotted
%% to each data matrix, to figure out which one it is,
%% and therefore what it should be called.  This seems very
%% inefficient, but it does work, and does not seem to take much
%% time at all.  Since draw is normally called by a function like
%% u.m, it would be easy to eliminate this search by simply making
%% u.m and similar functions pass in the variable name (or an index
%% into a name array) as an argument to draw. Or, yet another
%% global variable could even be used.
is_velocity = 1;
imask = ~isnan(r);
if r(imask) == AMP(imask)
   label = ' AMPLITUDE';
   is_velocity = 0;
elseif r(imask) == PGOOD(imask)
   label = ' PERCENT-GOOD';
   is_velocity = 0;
elseif r(imask) == U(imask)
   label = ' U-VELOCITY in m/s';
elseif r(imask) == V(imask)
   label = ' V-VELOCITY in m/s';
elseif r(imask) == W(imask)
   label = ' W-VELOCITY in m/s';
elseif r(imask) == U_ROTATE(imask)
   label = ' Forward-VELOCITY in m/s';
elseif r(imask) == V_ROTATE(imask)
   label = ' Port-VELOCITY in m/s';
elseif r(imask) == E(imask)
   label = ' ERROR-VELOCITY in m/s';
end


x = (r + DISP) + i * YAXIS(ones(NPRFS, 1), :).';
      % complex x spares us indexing into YAXIS

if is_velocity,
   x = x .* GOOD_PG;        % Never show things below the PGOOD_THRESHOLD
end                            % EF 99/02/24

if ~PLOTBAD_ORIGINAL
   x = x .* BADACC;          % blanking out low pgood & bins > lgb
                             % and anything already flagged
                             % in the database with flags
                             % other than the glitch flag.
end

if ~PLOTBAD

  ibad = find(BADBIN ~= 0);    % blanking out bad bins
  x(ibad) = NaN * ibad;

  ibot = find(BOTTOM <= NBINS);  % blanking out bins below bottom
  if ~isempty(ibot)
    botm = abs(BOTTOM);
    if BOTM_EDIT                % old: check BOTMPAS3 & LAST_85
      disp('Simulating effect of set_lgb')
      botm = botmpas3(botm);

      j = find(botm <= NBINS);
      botm(j) = floor(0.85 * botm(j));
    end
    for j = ibot
      x(botm(j)+1:NBINS,j) = NaN * [botm(j)+1:NBINS]';
    end
  end

  ibad = find(BADPRF ~= 0);    % blanking out bad profiles
  x(:,ibad) = NaN * ones(NBINS,length(ibad));
end

% stagger plots:
plot(x(bmin:bmax,pmin:pmax))
hold on

if (~isempty(IKEY)), plot(x(bmin:bmax,IKEY), '+w'), end

x = [1:NBINS]' * (-i) * ones(1,NPRFS) + (r+DISP);     % reinitialize
if is_velocity,
   x = x .* GOOD_PG;        % Never show things below the PGOOD_THRESHOLD
end                            % EF 99/02/24

% flag markings:
jj = pmin-1+find(BADPRF(pmin:pmax) == 2);
   % Sun  99/02/21 EF to keep bottom of bad profile from
   % skewing the  whole plot axis.
%%   ibad = find(isnan(BADACC));
%%   x(ibad) = NaN;
if (~isempty(jj)), plot(x(bmin:bmax, jj), wvar), end

px = NaN * x;
ibad = find(BADBIN == 2 | BADBIN == 6);
px(ibad) = x(ibad);
plot(px(bmin:bmax,pmin:pmax), d2)

px = NaN * x;
ibad = find(BADBIN == 4 | BADBIN == 6);
px(ibad) = x(ibad);
plot(px(bmin:bmax,pmin:pmax), ev)

px = NaN * x;
ibot = find(BOTTOM <= NBINS);
ibad = NBINS * (ibot-1) + abs(BOTTOM(ibot));
px(ibad) = x(ibad);
plot(px(bmin:bmax,pmin:pmax), amp)

px = NaN * x;
ibot = find(MAB <= NBINS);
ibad = NBINS * (ibot-1) + abs(MAB(ibot));
px(ibad) = x(ibad);
plot(px(bmin:bmax,pmin:pmax), old_amp)


px = NaN * x;
ibad = find(BADBIN < 0);
px(ibad) = x(ibad);
ibad = find(BADPRF < 0);
px(:,ibad) = x(:,ibad);
plot(px(bmin:bmax,pmin:pmax), user_added)

px = NaN * x;
ibad = find(mod(BADBIN,2) == 1);  % changed rem to mod
px(ibad) = x(ibad);
ibad = find(mod(BADPRF,2) == 1);
px(:,ibad) = x(:,ibad);
plot(px(bmin:bmax,pmin:pmax), db)

if PLOTBAD_ORIGINAL
   px = NaN * x;
   ibad = find(isnan(BADACC));
   px(ibad) = x(ibad);
   plot(px(bmin:bmax,pmin:pmax), range_pg)
   plot(px(bmin:bmax,pmin:pmax)) % Default colored lines also
end

hold off

ylabel('Bin Number')
set(gca, 'ygrid', 'on', 'YLim', [-BMAX-1, -BMIN+1])
ax = axis;
xs = abs(ax(2) - ax(1));
ys = abs(ax(4) - ax(3));

key_label = '';
if ~isempty(IKEY)
  if IKEY >= pmin & IKEY <= pmax
    key_label = '; + is ';
  end
end

text(ax(1)+0.7*xs, ax(3)-0.05*ys, ['Time: ' num2str(DAYS(pmin)) ' to ' num2str(DAYS(pmax)) key_label ...
   sprintf('%10.6f',DAYS(IKEY))])
text(ax(1)+0.7*xs, ax(3)-0.10*ys, ['B.P.: ' int2str(IBLKPRF(pmin,1)) ',' int2str(IBLKPRF(pmin,2)) ' to ' ...
   int2str(IBLKPRF(pmax,1)) ',' int2str(IBLKPRF(pmax,2)) key_label ...
   sprintf('%0.0f,%0.0f', IBLKPRF(IKEY,1), IBLKPRF(IKEY,2))])
text(ax(1), ax(3)-0.05*ys, ['Lon.: ' num2str(LON(pmin)) ' to ' num2str(LON(pmax)) key_label ...
   num2str(LON(IKEY))])
text(ax(1), ax(3)-0.10*ys, ['Lat.: ' num2str(LAT(pmin)) ' to ' num2str(LAT(pmax)) key_label ...
   num2str(LAT(IKEY))])
text(-0.15,-0.1,[pctime],'units','normalized','rotation',90.0)


[md nd] = size(DISP);
if nd >= 2
  title([DBNAME label '(offset = ' num2str(DISP(1,2)) ')'])
else % only 0 or 1 profile so no offset
  title([DBNAME label])
end

CURRX = r;
