function tspoke(action,varargin)
% Plot 2 timeseries against each other and time, higlight corresponding bits
%
% syntax xytpoke(data,time_field)
%
% data is a structure with data fields to be plotted
%
% time_field is either the name of the time field in data or a vector of
%            times
% 
% name,value pairs
%  x           xdata field name  (eg. 'lat')
%  y           ydata field name  (eg. 'lon')
% limits       'eq'
%
% ------------
% %getting userdata:
% figH = gcf;
% ud = get(figH,'UserData');
% ud.d.[fieldname](ud.sel)



% Switchyard calling pattern.  General format is 'tspoke(action,vars)'.
% If data rather than an action is first, initialize a new instance.
if ischar(action)			% This is a callback
  try
    feval(action, varargin{:});      	% Switchyard
  catch
    disp(lasterr);
  end
elseif isstruct(action)
  if ischar(varargin{1})
    tstr = varargin{1};
    if isfield(action,tstr)
      t = getfield(action,tstr);
    else
      error(sprintf('%s is not a field of the input data!',tstr));
    end
  else
    t = varargin{1};			% assume it's a vector
    tstr = 'Time';
  end
  if length(varargin)>1;
    xytpinit(action, t, tstr, varargin{2:end});
  else
    xytpinit(action, t, tstr);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               SUBROUTINES                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytpinit(d, t, tstr, varargin)

% Some defaults
ud.d = d;
ud.t = t;
ud.equalscales = 0;
ud.xscale = 'linear';
ud.yscale = 'linear';
ud.markersize = 2;
ud.marker = 'o';
ud.asciidate = 0;
ud.xcolor = [1 .4 .4];       %peach
ud.xselcolor = .9*[1 0 0];   %dark red
ud.ycolor = .9*[0 1 1];      %cyan
ud.yselcolor = [0 0 1];      %blue
ud.xycolor = .5*[1 1 1];     %gray
ud.xyselcolor = .8*[0 1 0];  %green


names = fieldnames(d);

% Go through and eliminate from consideration any fields not the same
% length as time
tlength = length(t);
oknames = zeros(size(names));
for ii = 1:length(names)
  if length(getfield(d,names{ii}))==tlength & isnumeric(getfield(d,names{ii}))
    oknames(ii)=1;
  end
end
if sum(oknames)<2
  error(['There aren''t two data fields with the length of the time' ...
	 ' field.']);
end
ud.names=names(find(oknames));
ud.xname = ud.names{1};
ud.yname = ud.names{2};
figtitle='';

%%JH 2002/06/24
if (isfield(d, 'date'))
   ud.date = d.date;  %cellarray of date strings
end

% Process the input arguments
for ii = 1:2:length(varargin)
  switch(lower(varargin{ii}))
   case {'x','x1'}
    ud.xname = varargin{ii+1};
    if ~isfield(d,ud.xname)
      error(sprintf('%s is not a valid data field',xname));
    end
   case {'y','x2'}
    ud.yname = varargin{ii+1};
    if ~isfield(d,ud.yname)
      error(sprintf('%s is not a valid data field',xname));
    end
   case {'limits'}
    if strmatch('eq',lower(varargin{ii+1}))
      ud.equalscales = 1;
    end
   case {'scale'}
    if strmatch('lo',lower(varargin{ii+1}))
      xlog = 1;
      ylog = 1;
    end
   case {'title'}
     figtitle = varargin{ii+1};
   case {'asciidate'}
    if strmatch('lo',lower(varargin{ii+1}))
       ud.asciidate = 1;
    end
   case {'xscale'}
    if strmatch('lo',lower(varargin{ii+1}))
      xlog = 1;
    end
   case {'yscale'}
    if strmatch('lo',lower(varargin{ii+1}))
      ylog = 1;
    end
   case {'marker'}
    ud.marker = varargin{ii+1};
   case {'markersize'}
    ud.markersize = varargin{ii+1};
  end
end
    
%Set up figure (created with GUIDE, the fig2mfile)
tspokegui;
figH = gcf;
if (~isempty(figtitle))
   set(figH, 'name', figtitle)
end

%Set up context menus for marker symbol and size
ux=uicontextmenu('parent',figH);
uy=uimenu('parent',ux,...
	  'label','Zoom');

uu(1)=uimenu('parent',uy,...
	  'label','Both');
uu(2)=uimenu('parent',uy,...
	  'label','XY plot');
uu(3)=uimenu('parent',uy,...
	  'label','Time plot');
for ii = 1:3
  uw(ii,1)=uimenu('parent',uu(ii),...
		  'label','Selection',...
		  'Tag',get(uu(ii),'Label'), ...
		  'callback','tspoke(''xytzoomselection'',gcbo)');
  uw(ii,2)=uimenu('parent',uu(ii),...
		  'label','Out',...
		  'Tag',get(uu(ii),'Label'), ...
		  'callback','tspoke(''xytzoomout'',gcbo)');
  uw(ii,3)=uimenu('parent',uu(ii),...
		  'label','All data',...
		  'Tag',get(uu(ii),'Label'), ...
		  'callback','tspoke(''xytzoomall'',gcbo)');
  uw(ii,4)=uimenu('parent',uu(ii),...
		  'label','Rectangle',...
		  'Tag',get(uu(ii),'Label'), ...
		  'callback','tspoke(''xytzoomrect'',gcbo)');
end

uu=uimenu('parent',ux,...
	  'label','Plot Symbol',...
	  'callback','tspoke(''xytsetsymbol'',gcbf)');

uu=uimenu('parent',ux,...
	  'label','Symbol Size',...
	  'callback','tspoke(''xytsetsymbolsize'',gcbf)');

uu=uimenu('parent',ux,...
	  'label','Linear Fit',...
	  'callback','tspoke(''fitselection'',gcbf)');

% fill data selection menus
xbH = findobj(figH,'tag','xbutton');
ybH = findobj(figH,'tag','ybutton');
set(xbH,'String',ud.xname,'Callback','tspoke(''tspchoosevar'',gcbo)');
set(ybH,'String',ud.yname,'Callback','tspoke(''tspchoosevar'',gcbo)');

% Top plot is x vs y
ud.x = getfield(d,ud.xname);
ud.y = getfield(d,ud.yname);
ud.ax(1) = findobj(figH,'Tag','xyaxes');
ud.xylineH = plot(ud.x,ud.y, ...
                  'color',ud.xycolor,'Marker',ud.marker, ...
                  'MarkerSize', ud.markersize', ...
                  'LineStyle', 'none', 'HitTest', 'off');
set(ud.ax(1), 'ButtonDownFcn', ...
	      'tspoke(''xytselectpoly'',gcbo,1)', ...
	      'XScale',ud.xscale, 'YScale', ud.yscale);

if ud.equalscales; axis equal; end

% Bottom plot is timeseries.  Separate y axes in case the variables are
% very different.  This requires 2 axes overlaid.
ud.ax(2) = findobj(figH,'Tag','tsaxes1');
axes(ud.ax(2));
ud.xtlineH = plot(ud.t,ud.x, ...
                  'color',ud.xcolor,'Marker',ud.marker, ...
                  'MarkerSize', ud.markersize', ...
                  'LineStyle', 'none', 'HitTest', 'off');
set(gca,'YColor',ud.xselcolor,'YScale',ud.xscale);
if (ud.asciidate), datetick('x','keeplimits'); end   %JH 2002/06/24
ylabel(ud.xname, 'Interpreter','none');
xlabel(tstr, 'Interpreter','none');

ud.ax(3) = axes('Units',get(ud.ax(2),'Units'), ...
		'Position',get(ud.ax(2),'Position'), ...
		'Color','none', ...
		'XTick',[], ...
		'YAxisLocation', 'right', ...
		'YColor',ud.yselcolor, ...
		'Tag','tsaxes2', ...
		'YScale',ud.yscale);
axes(ud.ax(3))
ud.ytlineH = line(ud.t,ud.y, ...
		  'color', ud.ycolor,'Marker',ud.marker, ...
		  'MarkerSize', ud.markersize', ...
		  'LineStyle', 'none', 'HitTest', 'off');
ylabel(ud.yname, 'Interpreter','none');
set(ud.ax(3), 'ButtonDownFcn','tspoke(''xytselectrect'',gcbo,3)');

if ud.equalscales			% Look for extremes
  ylims = get(ud.ax(2:3),'YLim');	% returns a cell array, which is
  ylims = cat(1,ylims{:});		% convertes to double.  Seems like
  newylim(1) = min(ylims(:,1));		% there ought to be a more
  newylim(2) = max(ylims(:,2));		% graceful way to do this
  set(ud.ax(2:3),'YLim', newylim);
end

ud.sel = [];
set(figH,'UserData',ud,'Pointer','fullcrosshair', 'UIContextMenu',ux);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytselectrect(axH,axnum)

point1 = get(gca,'CurrentPoint');    % button down detected
finalRect = rbbox;                   % return figure units
point2 = get(gca,'CurrentPoint');    % button up detected

minpt = min([point1; point2]);
maxpt = max([point1; point2]);

figH = gcbf;
ud = get(figH,'UserData');

%Find which points lie in the new selection region
if axnum == 1
  sel = find(ud.x>=minpt(1) & ud.x<=maxpt(1) & ...
	     ud.y>=minpt(2) & ud.y<=maxpt(2));
elseif axnum == 3
  sel = find(ud.t>=minpt(1) & ud.t<=maxpt(1));
end

xythighlight(figH,sel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xythighlight(figH,newsel)

ud = get(figH,'UserData');

%Delete any earlier highlighted selection
lineH = findobj(figH,'Tag','highlight');
if ~isempty(lineH)
  delete(lineH);
end

%Determine the selection mode
modeH = findobj(figH,'Tag','selectmode');
strs = get(modeH,'String');
switch(strs{get(modeH,'Value')})
 case {'New'}
  ud.sel = newsel;
 case {'New not'}
  ud.sel = setdiff(1:length(ud.x),newsel);
 case {'Add'}
  ud.sel = union(ud.sel,newsel);
 case {'Subset'}
  ud.sel = intersect(ud.sel,newsel);
 case {'Subset not'}
  ud.sel = setdiff(ud.sel,newsel);
end

%plot the new highlights
erasefit(figH);
axes(ud.ax(1))
ud.xyselH = line(ud.x(ud.sel), ud.y(ud.sel), 'Linestyle','none', ...
		 'marker', ud.marker, ...
		 'MarkerSize', ud.markersize, ...
		 'Color',ud.xyselcolor, ...
		 'Tag', 'highlight', 'HitTest', 'off');
axes(ud.ax(2))
ud.xtselH = line(ud.t(ud.sel), ud.x(ud.sel), 'Linestyle','none', ...
		 'marker',ud.marker, ...
		 'MarkerSize', ud.markersize, ...
		 'Color',ud.xselcolor, ...
		 'Tag', 'highlight', 'HitTest', 'off');
axes(ud.ax(3))
ud.ytselH = line(ud.t(ud.sel), ud.y(ud.sel), 'Linestyle','none', ...
		 'marker',ud.marker, ...
		 'MarkerSize', ud.markersize, ...
		 'Color',ud.yselcolor, ...
		 'Tag', 'highlight', 'HitTest', 'off');

set(figH,'UserData',ud);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y] = figpos2data(figH,axH,figpt);

% For the moment, assume that figure position is in pixels and axes
% position is normalized.
figpos = get(figH,'position');
fracfigpt = figpt./figpos(3:4);
axespos = get(axH,'position');
fracaxespt = (fracfigpt - axespos(1:2))./axespos(3:4);
xlim = get(axH,'XLim');
ylim = get(axH,'YLim');

if(strcmp('log',get(axH,'XScale')))
  xlog = 1;
  xlim = log10(xlim);
else
  xlog=0;
end
if(strcmp('log',get(axH,'YScale')))
  ylog = 1;
  ylim = log10(ylim);
else
  ylog=0;
end

x = xlim(1) + fracaxespt(1).*diff(xlim);
y = ylim(1) + fracaxespt(2).*diff(ylim);

if xlog; x = 10.^x; end
if ylog; y = 10.^y; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytsetsymbolsize(figH)

ud = get(figH,'UserData');
sizestr = inputdlg('Enter a new symbol size', ...
		   'Dialog: symbol size', 1, ...
		   {num2str(ud.markersize)});
if ~isempty(sizestr)			% not cancelled or empty
  newsize = str2double(sizestr);
  if newsize>0 & newsize<50		% basic sanity check.  Should
    ud.markersize = newsize;		% probably prompt for new value,
    h = findobj(figH,'Type','line');	% but I'm too lazy.
    set(h,'MarkerSize',newsize);
    set(figH,'UserData',ud);
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytsetsymbol(figH)

ud = get(figH,'UserData');
symbols = {'+', 'o', '*', '.', 'x', 'square', 'diamond', ...
	   'v', '^', '>', '<', 'pentagram', 'hexagram'};
currsymnum = strmatch(ud.marker, symbols, 'exact');

[newsymnum, ok] = listdlg('ListString', symbols, ...
			  'SelectionMode', 'single', ...
			  'InitialValue', currsymnum, ...
			  'PromptString', 'Choose plotting symbol');
if ok
  ud.marker = symbols{newsymnum};
  h = findobj(figH,'Type','line');
  set(h,'Marker',ud.marker);
  set(figH,'UserData',ud);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytselectpoly(axH,axnum)

figH = gcbf;

%Figure out which selection routine to use
toolH = findobj(figH,'Tag','selecttool');
switch(get(toolH,'Value'))
 case 1
  xytselectrect(axH,axnum);
  return
 case 2
  polyH = drawpolygon(axH);
 case 3
  polyH = drawfreehand(axH);
end

%polyH = feval(func,axH);
%if isempty(polyH)
%  return
%end

ud = get(figH,'UserData');

%Get the polygon data
poly.x = get(polyH,'XData');
poly.y = get(polyH,'YData');

% inpolygon is ungodly slow if the polygon has many vertices.  Try speeding
% things up by only giving it the data points within the extremes of the
% polygon. 
xmin = min(poly.x);
xmax = max(poly.x);
ymin = min(poly.y);
ymax = max(poly.y);

rect = find(ud.x>=xmin & ud.x<=xmax & ud.y>=ymin & ud.y<=ymax);
inside = find(inpolygon(ud.x(rect),ud.y(rect),poly.x,poly.y));
sel = rect(inside);

delete(polyH);
xythighlight(figH,sel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cursor = xcursor(mode)
if nargin<1; mode = 'data'; end;

if ~isempty(strmatch(lower(mode),{'hotspot','center'}))
  cursor = [8 8];
else
  cursor = ...
      [  2   2  NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN  2   2  NaN;
	 2   1   2  NaN NaN NaN NaN NaN NaN NaN NaN NaN  2   1   2  NaN;
	NaN  2   1   2  NaN NaN NaN NaN NaN NaN NaN  2   1   2  NaN NaN;
	NaN NaN  2   1   2  NaN NaN NaN NaN NaN  2   1   2  NaN NaN NaN;
	NaN NaN NaN  2   1   2  NaN NaN NaN  2   1   2  NaN NaN NaN NaN;
	NaN NaN NaN NaN  2   2  NaN NaN NaN  2   2  NaN NaN NaN NaN NaN;
	NaN NaN NaN NaN NaN NaN NaN  2  NaN NaN NaN NaN NaN NaN NaN NaN;
	NaN NaN NaN NaN NaN NaN  2   1   2  NaN NaN NaN NaN NaN NaN NaN;
	NaN NaN NaN NaN NaN NaN NaN  2  NaN NaN NaN NaN NaN NaN NaN NaN;
	NaN NaN NaN NaN  2   2  NaN NaN NaN  2   2  NaN NaN NaN NaN NaN;
	NaN NaN NaN  2   1   2  NaN NaN NaN  2   1   2  NaN NaN NaN NaN;
	NaN NaN  2   1   2  NaN NaN NaN NaN NaN  2   1   2  NaN NaN NaN;
	NaN  2   1   2  NaN NaN NaN NaN NaN NaN NaN  2   1   2  NaN NaN;
	 2   1   2  NaN NaN NaN NaN NaN NaN NaN NaN NaN  2   1   2  NaN;
 	 2   2  NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN  2   2  NaN;
	NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytzoomselection(menuH)
figH = gcbf;
%Check which axes to operate on
switch(get(menuH,'Tag'))
 case 'Both'
  doxy = 1;
  doxyt = 1;
 case 'XY plot'
  doxy = 1;
  doxyt = 0;
 case 'Time plot'
  doxy = 0;
  doxyt = 1;
end
ud = get(figH,'UserData');
xlim = [min(ud.x(ud.sel)) max(ud.x(ud.sel))];
ylim = [min(ud.y(ud.sel)) max(ud.y(ud.sel))];
tlim = [min(ud.t(ud.sel)) max(ud.t(ud.sel))];

if ud.equalscales
  dlim = [min([xlim(1) ylim(1)]) max([xlim(2) ylim(2)])];
  xlim = dlim;
  ylim = dlim;
end

erasefit(figH);

if doxy
  set(ud.ax(1),'XLim',xlim, 'YLim',ylim);
end
if doxyt
  set(ud.ax(2),'XLim',tlim, 'YLim',xlim);
  set(ud.ax(3),'XLim',tlim, 'YLim',ylim);
  axes(ud.ax(2));
  if (ud.asciidate), datetick('x','keeplimits'); end   %JH 2002/06/24
  axes(ud.ax(3));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytzoomout(menuH)
%zoom out a fixed percentage
zoomconst = 0.5;			% Double the span of each axis
figH = gcbf;
%Check which axes to operate on
switch(get(menuH,'Tag'))
 case 'Both'
  doxy = 1;
  doxyt = 1;
 case 'XY plot'
  doxy = 1;
  doxyt = 0;
 case 'Time plot'
  doxy = 0;
  doxyt = 1;
end
ud = get(figH,'UserData');
xlim = get(ud.ax(1),'XLim');
ylim = get(ud.ax(1),'YLim');
tlim = get(ud.ax(2),'XLim');
xlim = xlim+[-1 1].*diff(xlim).*zoomconst;
ylim = ylim+[-1 1].*diff(ylim).*zoomconst;
tlim = tlim+[-1 1].*diff(tlim).*zoomconst;
if ud.equalscales
  dlim = [min([xlim(1) ylim(1)]) max([xlim(2) ylim(2)])];
  xlim = dlim;
  ylim = dlim;
end

erasefit(figH);

if doxy
  set(ud.ax(1),'XLim',xlim, 'YLim',ylim);
end
if doxyt
  set(ud.ax(2),'XLim',tlim, 'YLim',xlim);
  set(ud.ax(3),'XLim',tlim, 'YLim',ylim);
  axes(ud.ax(2));
  if (ud.asciidate), datetick('x','keeplimits'); end   %JH 2002/06/24
  axes(ud.ax(3));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xytzoomall(menuH)
figH = gcbf;
ud = get(figH,'UserData');

%Check which axes to operate on
switch(get(menuH,'Tag'))
 case 'Both'
  doxy = 1;
  doxyt = 1;
 case 'XY plot'
  doxy = 1;
  doxyt = 0;
 case 'Time plot'
  doxy = 0;
  doxyt = 1;
end

erasefit(figH);

if doxy
  set(ud.ax(1),'XLimMode','auto','YLimMode','auto');
  if ud.equalscales
    axes(ud.ax(1))
    axis equal
  end
end

if doxyt
  set(ud.ax(2),'XLimMode','auto','YLimMode','auto');
  set(ud.ax(3),'XLimMode','auto','YLimMode','auto');
  axes(ud.ax(2));
  if (ud.asciidate), datetick('x','keeplimits'); end   %JH 2002/06/24
  axes(ud.ax(3));

  if ud.equalscales			% Look for extremes
    ylims = get(ud.ax(2:3),'YLim');	% returns a cell array, which is
    ylims = cat(1,ylims{:});		% convertes to double.  Seems like
    newylim(1) = min(ylims(:,1));		% there ought to be a more
    newylim(2) = max(ylims(:,2));		% graceful way to do this
    set(ud.ax(2:3),'YLim', newylim);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function tspchoosevar(buttonH)
% Change which variables are shown
figH = gcbf;
ud = get(figH,'UserData');
if strcmp(get(buttonH,'Tag'),'xbutton')
  tstr = 'X data';
elseif strcmp(get(buttonH,'Tag'),'ybutton')
  tstr = 'Y data';
end

curstr = get(buttonH,'String');
curidx = strmatch(curstr,ud.names,'exact');

listheight = min(500, 15*floor(length(ud.names)));

[choice,ok] = listdlg('ListString',ud.names, ...
                      'ListSize', [160 listheight],...
		      'SelectionMode','single', ...
		      'InitialValue', curidx, ...
		      'Name', 'tstr', ...
		      'PromptString', 'Choose data vector');

if ok==0
  return				% User hit 'Cancel'
end

name = ud.names{choice};
set(buttonH,'String',name)

erasefit(figH);

%Now replot the data
if strcmp(tstr,'X data')
  ud.xname = name;
  ud.x = getfield(ud.d,name);
  set(ud.xylineH,'XData',ud.x);
  set(ud.xtlineH,'YData',ud.x);
  if ~isempty(ud.sel)
    set(ud.xyselH,'XData',ud.x(ud.sel));
    set(ud.xtselH,'YData',ud.x(ud.sel));
  end
  axes(ud.ax(2));
  ylabel(name,'Interpreter','none')
else
  ud.yname = name;
  ud.y = getfield(ud.d,name);
  set(ud.xylineH,'YData',ud.y);
  set(ud.ytlineH,'YData',ud.y);
  if ~isempty(ud.sel)
    set(ud.xyselH,'YData',ud.y(ud.sel));
    set(ud.ytselH,'YData',ud.y(ud.sel));
  end
  axes(ud.ax(3));
  ylabel(name,'Interpreter','none')
end

%Reset scales to include all data
set(ud.ax(1), 'YLimMode','auto', 'XLimMode','auto');
set(ud.ax(2), 'YLimMode','auto', 'XLimMode','auto');
set(ud.ax(3), 'YLimMode','auto', 'XLimMode','auto');
axes(ud.ax(2));
if (ud.asciidate), datetick('x','keeplimits'); end   %JH 2002/06/24
axes(ud.ax(3));

set(figH,'UserData',ud);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fitselection(figH)
% Quick linear fit to the selected data.
ud = get(figH,'UserData');
if isempty(ud.sel)			% If nothing highlighted, quit
  return
end

%do the fit
ok = ~(isnan(ud.x(ud.sel)) | isnan(ud.y(ud.sel)));
if (length(ok) == 1)
   fprintf('\n\ndate is %s\n', ud.date{ud.sel(ok)});
   fprintf('dday is %10.7f\n', ud.t(ud.sel(ok)))
   fit.p = ud.t(ok);
   fit.R = 1;
else   
   fit.p = polyfit(ud.x(ud.sel(ok)),ud.y(ud.sel(ok)),1);
   %and get the correlation coefficient
   fit.R = corrcoef(ud.x(ud.sel(ok)),ud.y(ud.sel(ok)));
end

erasefit(figH);

if (length(fit.p) > 1)
   %Draw the line.  First get the axes limits to figure out line coefficients
   axes(ud.ax(1));
   xlim = get(ud.ax(1),'XLim');
   ylim = get(ud.ax(1),'YLim');
   textx = xlim(1)+0.05*diff(xlim);
   texty = ylim(1)+0.95*diff(ylim);
   textstr = sprintf('y = %fx%+f; R^2 = %4.2f',fit.p(1),fit.p(2),fit.R(1,2).^2);

   fit.line = line(xlim,polyval(fit.p,xlim), 'color', 'r');
   fit.text = text(textx,texty,textstr,'color','r');
end

setappdata(figH,'fitdata',fit);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function erasefit(figH)
%Get rid of fit line and text, if it exists.
if ~isappdata(figH,'fitdata')
  return
end
fit = getappdata(figH,'fitdata');
delete(fit.line)
delete(fit.text)
rmappdata(figH,'fitdata');
