   function axesHandle = underwayGrid(sbeFile, column)
% underwayGrid - Create a contour of data in Monterey Bay
%
% Use as:  underwayGrid(seabirdFile, column)
%       or underwayGrid(sbeFile, dataFile)
%
% Inputs: seabirdFile = name of the seaBird underway data file
%         column      = Column of data in the seaBird data file to grid
%         dataFile = An underway data file containing two columns of data. The first
%                    column is time (in Excel's format - serial days since 1900). The second
%                    column is the data points.
%
% Example: underwayGrid('Example_seabird.prc', 6)
%          underwayGrid('Example_seabird.prc', 'Example_pco2.prn')
%
% Note: All lines of the headers in any file used must be removed or commented out.
%       To Comment a header line add '%' at the start of the line.
%       Example files, Example_seabird.prc and Example_pco2.prn, are provided
%       for reference
% 
%       SeaBird files a fromatted with the first column is Serial Day of the Year,
%       the second column contains latitude in decimal degrees, and the third
%       column is longitude  in decimal degrees.

% Brian Schlining
% 13 Apr 2000

switch nargin
case 0
   fprintf(1, 'No input arguments were specified\n Use as underwayGrid(sbeFile, <column or dataFile>)\n')
   return
case 1
   fprintf(1, 'Use as underwayGrid(sbeFile, <column or dataFile>)\n')
   return
end

%*********************
% Load the data files
%*********************
x   = load(sbeFile);              % Load underway data from seabird
t   = jul2date(x(:,1), 2000*ones(size(x(:,1))));     % (1st Column contains JDay)
lat = x(:,2);                     % Grab the latitude from the 2nd column
lon = x(:,3);                     % Grab the longitude from the 3rd column

if isstr(column)                  % If the second input argument is a string it specifes a file of underway data (like pCO2)
   x   = load(column);            % Load the data, (a delimited file saved from an excel sheet)
   tx  = xl2sdn(x(:,1));          % Convert the xl times in the first column to matlab dates
   dat = x(:,2);                  % Second column contains data
   lat = interp1(t, lat, tx);     % Navigate the data base on the seabird data file
   lon = interp1(t, lon, tx);
else
   dat = x(:, column);
end


if ~isempty(x)

   % Marcelino!! Change the collowing line to adust color limits
   CLim = [prctile_(dat, 3) prctile_(dat, 97)];   % Specifies the prctile of data ato use for setting colorLimits

   %***************
   % Grid the data
   %***************
   good = find(~isnan(lat) & ~isnan(lon));
   lat  = lat(good);
   lon  = lon(good);
   dat  = dat(good);
   XI   = [min(lon):(max(lon) - min(lon))./100:max(lon)];
   YI   = [min(lat):(max(lat) - min(lat))./100:max(lat)]';
   ZI   = griddata(lon,lat,dat,XI,YI, 'linear');  % Type options are linear, cubic or nearest

   %***************
   % Create a plot
   %***************
   surf(XI,YI,ZI)                % Data Contour
   set(gca, 'CLim', CLim)
   shading interp
   view(2)
   hold on
   ZLim = get(gca, 'ZLim');
   plot3(lon, lat, ones(size(dat))*(max(ZLim)+1), 'k.', 'MarkerSize', .1)  % Plot sample positions
   load CST                      % Load the coastline data
   XLim = [-122.422358651448  -121.766379333516];
   YLim = [36.56187756710288  37.13674116292442];
   bad = find(CST(:,1) < XLim(1) | CST(:,1) > XLim(2) | CST(:,2) < YLim(1) | CST(:,2) > YLim(2));
   CST(bad) = NaN;
   plot3(CST(:,1), CST(:,2), ones(size(CST(:,1)))*(max(ZLim)+1), 'k') % Plot the coastline
   hold off
   colorbar2
   set(gca, 'XLim', XLim, 'YLim', YLim)

   %axis equal
   xlabel('Longitude')
   ylabel('Latitude')
   box on
   
   [r c] = size(ZI);
   XI = repmat(XI, r, 1);
   YI = repmat(YI, 1, c);
   dat = [XI(:), YI(:) ZI(:)];
   save temp.txt dat -ascii
end

if nargout
   axesHandle = gca;
end

function y = prctile_(x,p);
%PRCTILE gives the percentiles of the sample in X.
%   Y = PRCTILE(X,P) returns a value that is greater than P percent
%   of the values in X. For example, if P = 50  Y is the median of X.
%
%   P may be either a scalar or a vector. For scalar P, Y is a row
%   vector containing Pth percentile of each column of X. For vector P,
%   the ith row of Y is the P(i) percentile of each column of X.

[prows pcols] = size(p);
if prows ~= 1 & pcols ~= 1
    error('P must be a scalar or a vector.');
end

if any(p > 100) | any(p < 0)
    error('P must take values between 0 and 100');
end

xx = sort(x);
[m,n] = size(x);

if m==1 | n==1
    m = max(m,n);
   if m == 1,
      y = x*ones(length(p),1);
      return;
   end
    n = 1;
    q = 100*(0.5:m - 0.5)./m;
    xx = [min(x); xx(:); max(x)];
else
    q = 100*(0.5:m - 0.5)./m;
    xx = [min(x); xx; max(x)];
end

q = [0 q 100];
y = interp1(q,xx,p);






