function fig=sectionplot(section_number,record_number,varname,varargin)
%Author: Francois Counillon 30/08/2011 based on a routine from Knut Arild Lisaeter
%function fig=sectionplot(section_number,record_number,varname,varargin)
%A function which plots time record <record_number>
%from section <section_number> - Reads files generated by pak2sec
%
%
% Optional parameter/value pairs:
% 'showgrid'  - turn on(1)/off(0) hycom grid display, default is 0
% 'crange'    - Color range - vector with two values, max and min
% 'contours'  - vector - the values will be shown as contours on the plot
% 'staircase' - turn on(1)/off(0) staircase plot of variable
% 'xaxis'     - Choose x-axis variable - options: lon lat dist index - dist is default
% 'depthscale'- Factor for depth scaling
%
%
% Examples:
%    sectionplot(1,1,'temp') -- Plots temperature for section 1, time index 1
%    sectionplot(4,1,'temp') -- Plots temperature for section 4, time index 1
%
% Example with optional values, showing grid overlay and contour lines:
%    sectionplot(1,1,'temp','showgrid',1,'contours',[ 0 2 4 8 12]) 
%
% Example with optional values, using longitude along xaxis
%    sectionplot(1,1,'temp','showgrid',1,'xaxis','lon')
% Example with staircase plotting
%    sectionplot(1,1,'temp','staircase',1)
%
%

grid=0;
prange=[];
staircasebottom=0 ; 
Vcontours=[];
xaxis=[];
DIST=[];
depthscale=1.;
if(nargin>3 ) 
   for i=1:2:size(varargin,2)
      if (strcmp(varargin{i},'showgrid'))
         grid=varargin{i+1};
      end
      if (strcmp(varargin{i},'crange'))
         prange=varargin{i+1};
      end
      if (strcmp(varargin{i},'contours'))
         Vcontours=varargin{i+1};
      end
      if (strcmp(varargin{i},'staircase'))
         staircasebottom=varargin{i+1};
      end
      if (strcmp(varargin{i},'xaxis'))
         xaxis=varargin{i+1};
      end
      if (strcmp(varargin{i},'depthscale'))
         depthscale=varargin{i+1};
      end
   end
end
   

% Check for section file
csec=num2str(section_number,'%3.3d');
fname=['section' csec '.nc'];
nc=netcdf(fname);
if (isempty(nc))
   disp([' No file ' fname '. I Quit!']);
   return
end 
record=recdim(nc);
rec_dim=record(:);
if (rec_dim<record_number)
   disp([' Record exceeds number of records in netcdf file']);
   return
end 


% We should be set
isempty(xaxis);
if (strcmp(xaxis,'index'))
   DIST=1:numel(nc{'distance',1}(:));
   xtext='index';
   fac=1.;
elseif (strcmp(xaxis,'lon'))
   DIST=nc{'lon',1}(:);
   DIST=DIST';
   xtext='Longitude';
   fac=1.;
elseif (strcmp(xaxis,'lat'))
   DIST=nc{'lat',1}(:);
   DIST=DIST';
   xtext='Latitude';
   fac=1.;
elseif (isempty(xaxis) || strcmp(xaxis,'dist'))
   DIST=nc{'distance',1}(:)'; % Distance along section
   xtext='Distance [km]';
   fac=1e-3;
else
   disp(['Unknown x axis variable ' xaxis ]);
   disp('Options : lon, lat, index, dist and index. Default is dist')
   return ;
end

if (isempty(DIST)) 
   disp(['Could not get cumulative distance ']);
   return
end
% Depthc is variable for corner plotting
depthc=nc{'intf_lower',1}(record_number,:,:);
kdm=size(depthc,1);
dist2=repmat(DIST,kdm,1);
depthc2=zeros(kdm,size(depthc,2));
depthc2(2:kdm,:)=depthc(1:kdm-1,:);

depthc=depthc2;
if (isempty(depthc)) 
   disp(['Could not get depthc ']);
   return
end

% Get variable:
tmp=nc{varname,1}(record_number,:,:); 
if (isempty(tmp)) 
   disp(['Could not get variable ' varname]);

   % Show list of variables
   tst=var(nc);
   disp('Names of variables in netcdf file:')
   for i=1:numel(tst)
      disp(ncnames(tst(i)));
   end
   return
else
   %Display variable description if available
   tmpatt=nc{varname}.description(:);
   if (~isempty(tmpatt)) ;
      disp(['Short description of variable: ' tmpatt])
   end
end

% Treat empty layer:
for k=2:kdm
   I=find(tmp(k,:)<-20000.);
   tmpvar=tmp(k,:);
   tmpvar(I)=tmp(k-1,I);
   tmp(k,:)=tmpvar;
end

% Fix vertical arrangement so that lowest interface coincides with 
% lowest index into matrix -- this makes pcolor work better
mlim=0.01;
mask=find(depthc<mlim);
mask2=find(depthc>=mlim);

med=median(tmp(mask2));
mn=mean(tmp(mask2));
stdev=std(tmp(mask2));

%V holds min/max values
gcf; clf; F=gcf;
if (isempty(prange))
   V=[ mn-3*stdev    mn+3*stdev ];
else
   V=prange;
end
%V



%KAL -- always use pcolor
I=find(tmp>max(V)); 
tmp(I)=max(V);
I=find(tmp<min(V)); 
tmp(I)=min(V);
tmp(1,1)=min(V);
tmp(end,end)=max(V);

% This is if you want a "staircase" - type plot
if (staircasebottom==1)
   % Create new staircase arrays
   %size(dist)
   %size(depthf)
   stdepthc=zeros(size(depthc,1),2*numel(DIST)-1);
   for k=1:numel(DIST)
      for l=1:size(depthc,1);
         if (k<numel(DIST))
            stdist2 (l,2*k-1)    =dist2 (l,k);
            stdist2 (l,2*k  )    =dist2 (l,k+1);
            stdepthc(l,2*k-1:2*k)=depthc(l,k);
            sttmp   (l,2*k-1:2*k)=tmp   (l,k);
         else
            stdist2 (l,2*k-1)    =dist2 (l,k);
            stdepthc(l,2*k-1)    =depthc(l,k);
            sttmp   (l,2*k-1)    =tmp   (l,k);
         end
      end
   end 
end 


if (staircasebottom==1)
   P=pcolor(stdist2.*fac,stdepthc,sttmp);
else
   P=pcolor(dist2*fac,depthc2,tmp);
   shading flat;
end
set(P,'LineStyle','none')
A=gca;
set(A,'FontSize',11,'FontName','verdana','Color',[0.5 0.5 0.5]);
 xlabel(xtext);
ylabel('Depth (m)','FontSize',12);
title(['Variable ' varname ' record #: ' num2str(record_number)])
colorbar;
view([0 0 -1]);
hold on;
% showgrid
if (grid==1)
   for k=1:kdm
      if (staircasebottom)
         P3=plot(stdist2(kdm,:)*fac,stdepthc(k,:));
      else
         P3=plot(DIST*fac,depthc(k,:));
      end
      set(P3,'Color',[.7 .7 .7]);
      set(P3,'LineWidth',1);
   end
end

if (staircasebottom==1)
      P=stairs(stdist2(kdm,:)*fac,stdepthc(1,:));
      %P=stairs(stdist2(kdm,:)*fac,stdepthc(kdm,:));
   set(P,'Color','k');
   set(P,'LineWidth',2);
else
      P2=plot(DIST*fac,squeeze(depthc(kdm,:)));
   set(P2,'Color','k');
   set(P2,'LineWidth',2);
end


if (~isempty(Vcontours))
   [C,H]=contour(dist2*fac,depthc,tmp,Vcontours);
   set(H,'EdgeColor','k')
   set(H,'LineWidth',1)
   clabel(C,H,'FontSize',14,'FontName','verdana','Color','k','Rotation',0,'LabelSpacing',250)
end
% Return figure handle
close(nc);
print('-dpng',['Sec_' num2str(section_number) '_rec_' num2str(record_number) '_var_' varname '.png']);
