function x =cslice_interactive(cdf, var, time, res, scale, zfig)
%  CSLICE_INTERACTIVE produces a plot of a slice from an ECOMSI.CDF file.  The
%                user specifies the start and end locations using the mouse
%                the current graphics window.
%         
%                The amount of data returned depends upon the parameter 
%                'res', for resolution, which is the number of data points
%                from the starting point to ending point.  res must be >= 2.
%
%     The variable must be a function of x, y, z and time: (e.g. salinity)
%     USAGE cslice_interactive(cdf,var,time,res,scale,fig)
%        cdf = name of netcdf file
%        var = name of netcdf variable
%        time = time step, 0 by default
%        res = resolution = # of data points at which to interpolate, 10 by
%              default 
%        scale = how the data differs between the screen and memory
%        zfig = graphics handle for figure to take cross slice of

if (nargin < 2)
  help cslice_button; return
end
 
if (nargin < 3)
  time = 0;
  res = 10;
  scale = 1;
  zfig = gcf;
end
if (nargin < 4)
  res = 10;
  scale = 1;
  zfig = gcf;
end
if (nargin < 5)
  scale = 1;
  zfig = gcf;
end

% define the line handle to be initially zero
line_handle = 0;

mexcdf('setopts',0);
ncid=mexcdf('open',cdf,'nowrite');

% get x, y, and z dimensions
[name, nx]=mexcdf('diminq',ncid,'xpos');
[name, ny]=mexcdf('diminq',ncid,'ypos');
[name, nz]=mexcdf('diminq',ncid,'zpos');


% get the x and y grids
[xgrid,status] = mexcdf('varget', ncid, 'x', [0 0], [ny nx], 1);
[ygrid,status] = mexcdf('varget', ncid, 'y', [0 0], [ny nx], 1);


% Need to make get the bounding polygon defined by the grid.  This is needed
% to check whether or not the points chosen by the user are within the grid.
xp = zeros(1,2*nx + 2*ny - 4);
yp = zeros(1,2*nx + 2*ny - 4);

xp(1:nx-1) = xgrid(1,1:ny-1);
yp(1:nx-1) = ygrid(1,1:ny-1);

xp(ny:nx+ny-2) = xgrid(1:nx-1,ny);
yp(ny:nx+ny-2) = ygrid(1:nx-1,ny);

xp(nx+ny-1:2*ny+nx-3) = xgrid(nx,ny:-1:2);
yp(nx+ny-1:2*ny+nx-3) = ygrid(nx,ny:-1:2);

xp(2*ny+nx-2:2*nx+2*ny-4) = xgrid(nx:-1:2,1);
yp(2*ny+nx-2:2*nx+2*ny-4) = ygrid(nx:-1:2,1);



% open a new figure window for the cross slice plot
nfig = figure;

% loop until the sun dies or the user gets bored, which ever comes first
while 1

  % get the endpoints from the user after making the old figure current
  figure(zfig);
  [x,y] = ginput(2);

  
  % This looks like a good place to delete the new plot window
  if (nfig)
    figure(nfig);
    cla;
    figure(zfig);
  end
      
  p1 = [x(1) y(1)]*scale;  p2 = [x(2) y(2)]*scale;
  vec = p2-p1;
  
  
  % check the first point
  x = p1(1); y = p1(2);
  c = 0;
  j = 2*nx+2*ny-4;
  for i = 1:2*nx+2*ny-4
  
    if ( ( (yp(i)<=y) & (y<yp(j)) ) | ( (yp(j)<=y) & (y<yp(i)) ) ) 
      if (x < (xp(j) - xp(i)) * (y - yp(i)) / (yp(j) - yp(i)) + xp(i) )
        c = ~c;
      end 
    end

    j = i;
  end

  if (~c)
    fprintf(1,'point 1 was not within boundary....\n');
    fprintf(1,'If you are sure you placed it correctly, check your scale.\n');
  end

  % check 2nd point
  x = p2(1); y = p2(2);
  c = 0;
  j = 2*nx+2*ny-4;
  for i = 1:2*nx+2*ny-4

    if ( ( (yp(i)<=y) & (y<yp(j)) ) | ( (yp(j)<=y) & (y<yp(i)) ) ) 
      if (x < (xp(j) - xp(i)) * (y - yp(i)) / (yp(j) - yp(i)) + xp(i) )
        c = ~c;
      end
    end

    j = i;
  end

  if (~c)
    fprintf(1,'point 2 was not within boundary....\n');
    fprintf(1,'If you are sure you placed it correctly, check your scale.\n');
  end


  % Now locate the starting point.
  dist = abs(xgrid - p1(1)) + abs(ygrid - p1(2));
  index = find(dist==min(min(dist)));

  j1 = floor((index-1)/ny) + 1;
  i1 = index - (j1-1)*ny;

  % Now locate the ending point.
  dist = abs(xgrid - p2(1)) + abs(ygrid - p2(2));
  index = find(dist==min(min(dist)));

  j2 = floor((index-1)/ny) + 1;
  i2 = index - (j2-1)*ny;


%  hold on
  xxx = [xgrid(i1,j1) xgrid(i2,j2)];
  yyy = [ygrid(i1,j1) ygrid(i2,j2)];

  % everytime but the first time, delete the line
  if (line_handle)
    delete(line_handle);
  end
  line_handle = line(xxx/scale,yyy/scale,'color','white','erasemode','xor');
  hold off

  figure(nfig);
  [u,dd,xx,yy,zz] = cslice(cdf,var,i1,j1,i2,j2,time,res);
  shading interp;
  pslice(dd,zz,u);set(gca,'aspect',[nan nan]);
  %extcontour(dd,zz,u,'label','fill');
  mexcdf('close',ncid);

  x = 1;

end % while

return;

