function [cs,han]=contourfP(x,y,A,ci,Fcol,xBrk,Lcol,Lthick,xlohi);
% function [cs,han]=contourfP(x,y,A,ci,Fcol,xBrk,Lcol,Lthick,xlohi);
% b klinger, 5/15/15
% Version of contourf.m with more control over color palette.
% (x,y,A,ci) as in contourf.m
% Fcol = fill color = string of letters for sequence of colors
%          D=default;A=gray,B=blue,R=red,G=green,O=orange,N=brown
%          uppercase Fcol: light-to-dark; 
%          lowercase Fcol: dark-to-light
% Lcol string = line color = 'none',single letter (b,k,etc)
%           or 3 letters for +/-/0 (such as 'rbg')
% xBrk = vector of contour breaks: vector with value of ci
%          at which palette switches
% Lthick = line thickness
% xlohi = variable for colorpal.m (see documentation for colorpal.m)
% Calls colorpal, maptoint

%          MODIFICATIONS
% 5/28/15 allow to use numeric information in for Lcol when +/-/0 
%         and correct 'w' in +/-/0 option for printing
% 12/21/15 fix bug that happens for +/-/0 option when no values in
%         range for one of the contour calls

%---------------------------------------------------------------------
% some graphics overhead and default values
%---------------------------------------------------------------------

ihold=ishold;   % remember initial hold state of figure
if ~exist('Lthick')
   Lthick=.5;
end
if ~exist('Lcol')
   Lcol='k';
end
LW='linewidth';
LC='linecolor';

%---------------------------------------------------------------------
% decide if different color lines for +/-/0; if true create params
%---------------------------------------------------------------------

nci=length(ci);
Afill=A;
LcolF=Lcol;
posneg=0;
if length(Lcol)==3
   posneg=1;                % remember to do +/-/0
   ciF=[0 ci];              % fill contour values positive only
   cin=-ci(nci:-1:1);       % negative contours
   ncol=nci;                % number of shades 
   Afill=abs(A);            % array to do fill contours
   LcolF='none';            % no lines when doing filled contours
else           % if single color, convert from white to almost-white
   if LcolF=='w'
      LcolF=[1 1 1]*.99;
   end
end      

%---------------------------------------------------------------------
% create colors and (if necessary) transformation for filled contours
%---------------------------------------------------------------------

if upper(Fcol(1))=='D'        % if using default colormap...
   Amap=Afill;
   cimap=ci;
elseif upper(Fcol(1))=='X'
   cimap=1:nci;
   Amap=maptoint(A,ci);
else                         % if not using default colormap...

   Ndir=length(Fcol);
   if ~exist('xlohi')
      xlohi=[zeros(Ndir,1) ones(Ndir,1)];
   end

   iBrko=1;
   for idir=1:Ndir-1
      iBrk=max(find(ci<=xBrk(idir)));
      ncol(idir)=iBrk-iBrko;
      iBrko=iBrk;
   end
   ncol(Ndir)=nci-iBrko;

   colpal=colorpal(Fcol,ncol,xlohi);
   colormap(colpal);
   Amap=maptoint(Afill,ci);
   cimap=1:nci;

end

%---------------------------------------------------------------------
% Do filled contour plot.
% If +/-/0 option, do contours in different colors for +/-/0 values.
% For consistency, do contours using same mapped array as for fill.
%---------------------------------------------------------------------

if upper(Fcol(1))=='X'
   [cso,han]=contour(x,y,Amap,cimap,LC,LcolF,LW,Lthick);
else
   [cso,han]=contourf(x,y,Amap,cimap,LC,LcolF,LW,Lthick);
end
cs=maptoint(cso,ci);
if posneg

% check to see if character is a number, in which case convert
% to a RGB value for a gray shade (0=black,9=almost white)

   if length(strfind('0123456789',Lcol(1)))==0
      Lcolp=Lcol(1);
      if Lcolp=='w', Lcolp=[1 1 1]*.99; end
   else
      Lcolp=[1 1 1]*str2num(Lcol(1))/10;
   end
   if length(strfind('0123456789',Lcol(2)))==0
      Lcoln=Lcol(2); 
      if Lcoln=='w', Lcoln=[1 1 1]*.99; end
   else
      Lcoln=[1 1 1]*str2num(Lcol(2))/10;
   end
   if length(strfind('0123456789',Lcol(3)))==0
      Lcol0=Lcol(3);
      if Lcol0=='w', Lcol0=[1 1 1]*.99; end
   else
      Lcol0=[1 1 1]*str2num(Lcol(3))/10;
   end

% do contour plots

   hold on
   Apos=maptoint(A,ci);
   if length(find(isnan(Apos)==0))>0
      contour(x,y,Apos,1:nci,LC,Lcolp,LW,Lthick);
   end
   cineg=-ci(nci:-1:1);
   Aneg=maptoint(A,cineg);
   if length(find(isnan(Aneg)==0))>0
      contour(x,y,Aneg,1:nci,LC,Lcoln,LW,Lthick);
   end
   contour(x,y,A,[0 0],LC,Lcol0,LW,Lthick);
end

caxis(minmax(cimap))
if ihold==0; hold off; end     % return hold state to value at start
