function [x,XI,YI,ZI,ZI1,ZI2]=profilex(X,Y,sp,Z,Z1,Z2)
%
% [x,XI,YI,Zout1,Zout2,Zout3]=profilex(X,Y,sp,Z,Z1,Z2)
%
% profilex extracts profile, at spacing sp, from 
% the grid X,Y,Z (same profile from Z1 Z2 if present)
% X,Y,Z should already be contoured in the current figure window
% when profilex is called

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Nov, 1997)  
%                             
%***********************************************************************************
%  
disp('EXTRACT A PROFILE FROM THE GRID')
disp('CLICK ON THE BEGINNING AND END POINTS OF THE PROFILE')
[x_x,y_y]=ginput(2);
ll=sqrt((x_x(2)-x_x(1))^2 + (y_y(2)-y_y(1))^2); % length of the profile

  if ~exist('sp') % if the sampling interval, sp, does not exist prompt    
  sp=input('ENTER THE SAMPLING INTERVAL ALONG THE PROFILE ')
  end % end if

nx=round(ll/sp); % num of samples
r=[0:nx]*sp; % sample locations along profile
flip='n'; % flag for profile orientation

  if(x_x(2))<x_x(1)  % if the profile is easternmost to westernmost
  x_x=flipud(x_x);   % flip the ends
  y_y=flipud(y_y);
  r=[nx:-1:0]*sp;    % new sample locations                                      
  flip='y';          % flag that the ends have been flipped
  end                                         

slope=(y_y(2)-y_y(1))/(x_x(2)-x_x(1)); % slope of profile in xy plane
XI=x_x(1)+[0:nx]*sp*cos(atan(slope));  % x locations
YI=y_y(1)+[0:nx]*sp*sin(atan(slope));  % y locations
size(X);size(Y);size(Z);
ZI=interp2(X,Y,Z,XI,YI,'linear');      % interpolated values
if nargin==5;ZI1=interp2(X,Y,Z1,XI,YI,'linear');elseif nargin==6;ZI2=interp2(X,Y,Z2,XI,YI,'linear');ZI1=interp2(X,Y,Z1,XI,YI,'linear');end
clf
hgrid=axes('position',[.1 .4 .8 .55]); % axes for grid
contour(X,Y,Z,25),hold on
plot([x_x(1) x_x(2)],[y_y(1) y_y(2)],'-','LineWidth',3) % plot profile on grid

  if flip=='y';      % if the profile ends were flipped, flip back
  x_x=flipud(x_x);
  y_y=flipud(y_y);
  end                                         

text(x_x(1)+sp,y_y(1)-sp,'S') % print S(tart) of profile
text(x_x(2)+sp,y_y(2)-sp,'E') % print E(nd) of profile  
hprofile=axes('position',[.1 .1 .8 .2]);  % axes for profile
plot(r,ZI,'o')  % plot profile

x=r;  % copy of positions along profile
  if flip=='y';r=fliplr(r); end % if profile ends were flipped, flip r
text(r(1)+sp,ZI(1)+(max(ZI)-min(ZI))/10,'S')          % print S(tart) of profile
text(r(length(r))-sp,ZI(1)+(max(ZI)-min(ZI))/10,'E')  % print E(nd) of profile  

set(gcf,'currentaxes',hgrid) % make current axes the grid