D=readkorrel('korrel.txt');

% USER options;
% ------------

   CLRMAP   = colormap;close;
   ax       = [90 300 535 635]*1e3; % whole wad
   ax       = [90 190 530 630]*1e3; % wwad

   U.clrscale  = 4;
   U.cmax      = 10;
   U.dotsize   = 4;

   U.exportfig = 0;
   U.printfig  = 0;
   U.closefig  = 0;
   U.specials  = 1;
   U.sieves    = 0;
   U.exportxyz = 0;


   D.wwad =  INPOLYGON(D.X_COORD,...
                       D.Y_COORD,...
                       [ax(1) ax(2) ax(2) ax(1)],...
                       [ax(3) ax(3) ax(4) ax(4)]);

   UNITVEC  = (0:length(CLRMAP)-1)./(length(CLRMAP)-1);

if U.specials
   
   for pmudcases = 2
   
   tic;

      switch pmudcases

         case 1
         D.pmud = D.cumPHI(13,:)/D.cumPHI(17,:).*D.pmud16_from_63;
         titletxt = {['pMud_{<63um} = pMud_{<63um, sedatlas} \times pMud_{<16um, EVA} / pMud_{<16um, sedatlas}'],...
                     ['where pMud_{<16um, EVA} = f(pMud_{<63um, sedatlas})']};
         case 2
         D.pmud = D.pmud16_from_63;
         titletxt = {['pMud_{<16um, EVA} = f(pMud_{<63um, sedatlas})']};
         
         case 3
         D.pmud = D.cumPHI(13,:)/D.cumPHI(17,:).*D.pmud16_from_16;
         titletxt = {['pMud_{<63um} = pMud_{<63um, sedatlas} \times pMud_{<16um, EVA} / pMud_{<16um, sedatlas}'],...
                     ['where pMud_{<16um, EVA} = f(pMud_{<16um, sedatlas})']};
         case 4
         D.pmud = D.pmud16_from_16;
         titletxt = {['pMud_{<16um, EVA} = f(pMud_{<16um, sedatlas})']};

         case 5
         D.pmud = D.cumPHI(13,:)';
         titletxt = ['pMud_{<63um, sedatlas}'];

         case 6
         D.pmud = D.cumPHI(17,:)';
         titletxt = ['pMud_{<16um, sedatlas}'];
      
      otherwise
      
      end

      figmud = figure;
      %setfig2screensize(figmud)
      hold     ('on')

      % scale on max per sieve
      % -------------------
      if U.clrscale  == 1
         fullrange = UNITVEC.*(1.01*max(D.pmud)-0.995*min(D.pmud)) + 0.995*min(D.pmud);
      
      % scale per on 100 %
      % -------------------
      elseif U.clrscale  == 2
         fullrange = UNITVEC.*100;
      
      % scale on max for all sieves
      % -------------------
      elseif U.clrscale  == 3
         fullrange = UNITVEC.*maxmax(D.pmud);

      % scale on U.cmax
      % -------------------
      elseif U.clrscale  == 4
         fullrange = UNITVEC.*U.cmax;
      end
      
      % make colorscale
      % -------------------
      D.C(:,1)  = interp1(fullrange,CLRMAP(:,1),D.pmud);
      D.C(:,2)  = interp1(fullrange,CLRMAP(:,2),D.pmud);
      D.C(:,3)  = interp1(fullrange,CLRMAP(:,3),D.pmud);
      
      leave_in = or(or(~isnan(D.C(:,1)),...
                       ~isnan(D.C(:,2))),...
                       ~isnan(D.C(:,3)));
      
      surf([-1e6 -2e6;-1e6 -2e6],[-1e6 -1e6;-2e6 -2e6],[0 0;max(fullrange) max(fullrange)]);
      %colorbarwithtxt('% on sieve')
      colorbar
      H = zeros(size(D.X_COORD));
      
      for i=1:length(D.X_COORD)
         if D.wwad(i)
            if leave_in(i)
                H(i) = plot(D.X_COORD(i),D.Y_COORD(i),'o',...
                            'markerfacecolor',D.C(i,:),...
                            'markeredgecolor',D.C(i,:),...
                            'Color'          ,D.C(i,:),...
                            'markersize'     ,U.dotsize); 
            end
         end
      end
      
      axis   ('equal')
      title  (titletxt);
      axis   (ax)
      
      if U.exportxyz         
        XYZ = [D.X_COORD D.Y_COORD, D.pmud];
        save(['pmud',num2str(pmudcases),'.xyz'],'XYZ','-ascii')
        % * data from sedimentartas waddenzee
	% * fraction < 63 um used with calibration curve
	% * 0.1879.*sum(D.PHI(13:end,:),1)'; % mud if D < 63 um
	% * to obtain the classical mud content definded 
	% * as the fraction < 16 um.
      end
      
      if U.printfig
         print(figmud,'-dpng', ['figmud',num2str(pmudcases),'colorscale',num2str(U.clrscale)])
      end      
      
      if U.closefig
         close(figmud)
      end
      
      disp(['Processed case ',num2str(pmudcases)]);

   toc 
   
   end % for pmudcases = 1:6
   
   
end % if U.specials

   if U.sieves
   
      for phi=1:26
      
         fig(phi) = figure   (phi);
         %setfig2screensize(fig(phi))
         hold     ('on')
         
         % scale on max per sieve
         % -------------------
         if U.clrscale  == 1
            fullrange = UNITVEC.*(max(D.PHI(phi,:))-min(D.PHI(phi,:))) + min(D.PHI(phi,:));
         
         % scale per on 100 %
         % -------------------
         elseif U.clrscale  == 2
            fullrange = UNITVEC.*100;
         
         % scale on max for all sieves
         % -------------------
         elseif U.clrscale  == 3
            fullrange = UNITVEC.*maxmax(D.PHI);
         end
         
         % make colorscale
         % -------------------
         D.C(:,1)  = interp1(fullrange,CLRMAP(:,1),D.PHI(phi,:));
         D.C(:,2)  = interp1(fullrange,CLRMAP(:,2),D.PHI(phi,:));
         D.C(:,3)  = interp1(fullrange,CLRMAP(:,3),D.PHI(phi,:));
      
         % plot colored dots
         % -------------------
      
         surf([-1e6 -2e6;-1e6 -2e6],[-1e6 -1e6;-2e6 -2e6],[0 0;max(fullrange) max(fullrange)]);
         %colorbarwithtxt('% on sieve')
         colorbar
         H = zeros(size(D.X_COORD));
         
         for i=1:length(D.X_COORD)
            if D.wwad(i)
               H(i) = plot(D.X_COORD(i),D.Y_COORD(i),'o',...
                           'markerfacecolor',D.C(i,:),...
                           'markeredgecolor',D.C(i,:),...
                           'Color'          ,D.C(i,:),...
                           'markersize'     ,U.dotsize); 
            end
         end
         
      
         disp(['Time needed for PHI',num2str(D.PHIvalue(phi,:)),' : ',num2str(toc),' s'])
         
         % lay out
         % -------------------
      
         axis   ('equal')
         title  (['Sieve: PHI',num2str(D.PHIvalue(phi,:)),' - Diameter: ',num2str(D.D(phi,:)),' [mm]'])
         axis   (ax)
      
         % print - export
         % -------------------
      
         if U.exportfig
            saveas (fig(phi),['PHI',num2str(D.PHIvalue(phi,:)),'colorscale',num2str(U.clrscale),'.fig'])
         end
      
         if U.printfig
            print  (fig(phi),'-dpng', ['PHI',num2str(D.PHIvalue(phi,:)),'colorscale',num2str(U.clrscale)])
         end
         
         if U.closefig
            close(fig(phi))
            fig(phi) = nan;
         end
         
      end
      
   end

%  Caros PC                            WL PC
%  Time needed for PHI-20 : 7.872 s    13.118 s
%  Time needed for PHI-15 : 5.818 s     5.557 s
%  Time needed for PHI-10 : 6.229 s     5.528 s
%  Time needed for PHI-5 :  5.598 s     5.558 s
%  Time needed for PHI0 :   5.868 s     5.978 s
%  Time needed for PHI5 :   6.209 s     6.569 s
%  Time needed for PHI10 :  6.31 s      9.463 s
%  Time needed for PHI15 :  8.172 s    15.341 s
%  Time needed for PHI20 : 17.716 s    43.339 s
%  Time needed for PHI25 : 21.17 s     43.339 s
%  Time needed for PHI30 : 21.221 s    45.222 s
%  Time needed for PHI35 : 11.006 s    20.337 s
%  Time needed for PHI40 :  9.033 s    10.324 s
%  Time needed for PHI45 :  6.489 s     8.863 s
%  Time needed for PHI50 :  6.879 s     7.64 s
%  Time needed for PHI55 :  6.589 s    10.074 s
%  Time needed for PHI60 :  7.451 s     8.482 s
%  Time needed for PHI65 :  6.91 s     18.996 s
%  Time needed for PHI70 :  8.492 s    13.529 s
%  Time needed for PHI75 :  6.539 s     9.513 s
%  Time needed for PHI80 :  6.61 s      7.691 s
%  Time needed for PHI85 :  6.249 s     7.24 s
%  Time needed for PHI90 :  6.329 s     6.308 s
%  Time needed for PHI95 :  6.529 s     5.808 s
%  Time needed for PHI100 : 6.559 s     6.129 s
%  Time needed for PHI105 : 6.51 s      5.658 s




























