% input spreadsheet should contain subsection of 'plutonic mantle' worksheet
% cut/paste columns K-U (use Paste Special, Values)
% delete Junk Basket lines at top
% save file as .csv format

%--------------------- get file & run parameters
%fin=input('Lithology unit spreadsheet (name): ');
%fout=input('Save mean values in file (name):  ');
%dz=input('Interval to average (m): ');
%rstep=input('Step size- shift averaging interval down core by- (m): ');
%iave=input('Compute as percent of drilled (=1) or recovered (=2) interval: ');

%hardwire settings for repeated file/parameter runs
fin='1473lith_0126a.csv';
fout='chk.txt';
dz=20;
rstep=1;
iave=1;

%--------------------- data read
% read in subtable with depth, domain percent, and lithology
% z1 is top of subsection (mbsf)
% z2 is bottom of subsection (mbsf)
% dnum is domain number within interval (lithology unit)
% dperc is fraction of core section that is this domain
%   >> unspecified values are set to NaN, initially
% rock is lithology name

lith_in=readtable(fin);
z1=table2array(lith_in(:,1));
z2=table2array(lith_in(:,2));
dnum=table2array(lith_in(:,3));
dperc=table2array(lith_in(:,5));
rock=table2array(lith_in(:,10));

% for 1105lith_1219.csv
%z2=table2array(lith_in(:,2));
%dnum=table2array(lith_in(:,3));
%dperc=table2array(lith_in(:,4));
%rock=table2array(lith_in(:,5));

nrec=length(z1);

%-------------------- define/group categories
% number of categories=nk
nk = 13;
rtype(nrec,1:8)='        ';

for jj=1:nrec,
   rname=cell2mat(rock(jj));
    if ( strcmp(rname, 'diorite')==1 || ...
        strcmp(rname, 'trondhjemite')==1 || ...
        strcmp(rname, 'oxide-bearing trondhjemite')==1 || ...
        strcmp(rname, 'felsic vein')==1 || ...
        strcmp(rname, 'felsic material')==1 || ...
        strcmp(rname, 'felsic patch')==1 ),
        rtype(jj,1:8)='Felsic  ';  
        
    elseif (strcmp(rname, 'diabase')==1 || ...
            strcmp(rname, 'basalt')==1 ),
       rtype(jj,1:8)='Diabase ';

    elseif (strcmp(rname, 'oxide gabbro')==1 || ...
            strcmp(rname, 'amphibole oxide gabbro')==1 || ...
            strcmp(rname, 'olivine oxide gabbro')==1 || ...
            strcmp(rname, 'olivine oxide microgabbro')==1 || ...
            strcmp(rname, 'olivine-bearing oxide microgabbro')==1 || ...
            strcmp(rname, 'olivine-bearing oxide gabbro')==1 || ...
            strcmp(rname, 'oxide olivine gabbro')==1 ),
       rtype(jj,1:8)='XGabbro ';

    elseif ( (strcmp(rname, 'disseminated-oxide olivine gabbro')==1 ) || ...
     (strcmp(rname,'disseminated-oxide olivine-bearing gabbro')==1) || ...
     (strcmp(rname,'disseminated-oxide gabbro')==1) || ...
   (strcmp(rname,'disseminated-oxide orthopyroxene-bearing olivine gabbro')==1) || ...
   (strcmp(rname,'disseminated oxide- and orthopyroxene bearing olivine gabbro')==1)),
      rtype(jj,1:8)='dxGabbro';

    elseif ( (strcmp(rname, 'oxide-bearing olivine gabbro')==1 )|| ...
            (strcmp(rname, 'oxide-bearing gabbro')==1) || ...
            (strcmp(rname, 'oxide- and orthopyroxene-bearing olivine gabbro')==1) || ...
            (strcmp(rname, 'oxide- and olivine-bearing gabbro')==1) || ...
            (strcmp(rname, 'olivine- and oxide bearing olivine gabbro')==1 ) || ...
            (strcmp(rname, 'olivine- and oxide-bearing gabbro')==1 ) ),
       rtype(jj,1:8)='xbGabbro';

    elseif (strcmp(rname, 'gabbro')==1 || ...
             strcmp(rname, 'olivine- and opx-bearing gabbro')==1 || ...
           strcmp(rname, 'olivine-bearing gabbro')==1 || ...
            strcmp(rname, 'microgabbro')==1 ),
       rtype(jj,1:8)='Gabbro  ';

    elseif (strcmp(rname, 'olivine gabbro')==1 || ...
            strcmp(rname, 'amphibole-bearing olivine gabbro')==1 || ...
            strcmp(rname, 'orthopyroxene-bearing olivine gabbro')==1 || ...
            strcmp(rname, 'olivine microgabbro')==1 ),
       rtype(jj,1:8)='OGabbro ';
       
    elseif (strcmp(rname, 'oxide-gabbronorite')==1 || ...
            strcmp(rname, 'olivine-bearing oxide gabbronorite')==1 || ...
            strcmp(rname, 'oxide gabbronorite')==1),
       rtype(jj,1:8)='XNGabbro';
       
    elseif (strcmp(rname, 'gabbronorite')==1 || ...
            strcmp(rname, 'norite')==1 ),
       rtype(jj,1:8)='NGabbro ';
       
    elseif (strcmp(rname, 'troctolitic gabbro')==1 ),
       rtype(jj,1:8)='TGabbro ';

    elseif ( (strcmp(rname, 'troctolite')==1 ) || ...
            (strcmp(rname, 'clinopyroxene-rich troctolite')==1) ),      
        rtype(jj,1:8)='Trocto  ';

    elseif (strcmp(rname, 'dunite')==1 ), ...
       rtype(jj,1:8)='Dunite ';
           
    elseif (strcmp(rname, 'clinopyroxenite')==1 || ...
            strcmp(rname, 'harzburgite')==1 || ...
            strcmp(rname, 'wehrlite')==1), ...
         rtype(jj,1:8)='Perid   ';
         
    elseif (strcmp(rname, 'rubble' )),
       rtype(jj,1:8)='Rubble';
    end
end

%------------------- determine number of averages needed (nint)
zmin = z1(1);
zmax = z2(nrec);
ndz = floor( (zmax-zmin-dz)/dz );
nadd = floor( (zmax -zmin -dz - ndz*dz)/rstep);
nint=ndz*dz/rstep+nadd;
zold = zmin;
zfind = zold + rstep;
znext = zmin + dz;

%--------------- convert domain percent to fraction
%  >> convert NaN values to 100% (fraction=1.0)

for jj=1:nrec
  if ( isnan(dperc(jj)) )
    dperc(jj)=1;
  else
    dperc(jj) = dperc(jj)/100;
  end
end

%------------ loop thru averaging windows, steping downhole 

jct = 1;
jj=1;
ii=1;
rt(1,1:nk) = 0.;
zmid(nint)=0.;
zmid(1) = zmin + dz/2.;

zrock(1:nint)=0.;
rt(1:nint,1:nk) = 0.;

%------------- determine if subsections are in current window
% keep track of partial subsection at top/bottom of averaging window

 while( (jj < nrec+1) && (ii < nint+1) ),

   if ( (jj >1) && (jj < nrec) ),

     while (((z1(jj) < zfind) && (z2(jj) < zfind) )==1 || ...
         ((z1(jj) == z1(jj+1)) && (dperc(jj) == 1))), 
       jj=jj+1;
       if ((z1(jj)==z1(jj-1)) && (dperc(jj)==1)), 
         rtype(jj,1:8)=rtype(jj-1,1:8);
       end
     end

	 if  (z1(jj) ~= z1(jj-1)) ,
        if ( (z1(jj-1) <= zfind) && (z1(jj) > zfind) ),
	      if (z2(jj-1) > zfind),
  	        jct = jj-1;
	      else
	        jct = jj;
	      end
        end
     elseif (dperc(jj) == 1),
        if ( (z1(jj-2) <= zfind) && (z1(jj) > zfind) ),
	      if (z2(jj-2) > zfind),
  	        jct = jj-2;
	      else
	        jct = jj;
	      end
	    end
     else
        if ( (z1(jj-dnum(jj)) <= zfind) && (z1(jj) > zfind) ),
	      if (z2(jj-dnum(jj)) > zfind),
  	        jct = jj-dnum(jj);
	      else
	        jct = jj;
	      end
	    end
	 end

   end

    ztop = z1(jj);
    if (z1(jj) < zold), ztop = zold; end
    if (z1(jj) < znext), 
      zinc = (z2(jj)-ztop); 
      if (z2(jj) > znext), zinc = znext - ztop; end
    end
    
%--------------- add up portions of each rock type within window

    if (strcmp(rtype(jj,1:8) , 'Felsic  ')==1),
       rt(ii,13) = rt(ii,13) + zinc*dperc(jj);
     
    elseif (strcmp(rtype(jj,1:8), 'Diabase ')==1),
       rt(ii,12) = rt(ii,12) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'XGabbro ')==1),
       rt(ii,11) = rt(ii,11) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'xbGabbro')==1 ),
       rt(ii,10) = rt(ii,10) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'dxGabbro' )==1),
       rt(ii,9) = rt(ii,9) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'XNGabbro')==1),
       rt(ii,8) = rt(ii,8) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'NGabbro')==1),
       rt(ii,7) = rt(ii,7) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'Gabbro  ')==1),
       rt(ii,6) = rt(ii,6) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'OGabbro ')==1),
       rt(ii,5) = rt(ii,5) + zinc*dperc(jj);
       
    elseif (strcmp(rtype(jj,1:8), 'TGabbro ')==1),
       rt(ii,4) = rt(ii,4) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'Trocto  ')==1),
       rt(ii,3) = rt(ii,3) + zinc*dperc(jj);
       
    elseif (strcmp(rtype(jj,1:8), 'Dunite ' )==1),
       rt(ii,2) = rt(ii,2) + zinc*dperc(jj);

    elseif (strcmp(rtype(jj,1:8), 'Perid   ' )==1),
       rt(ii,1) = rt(ii,1) + zinc*dperc(jj);

	else
	  fprintf('Unclassified: %d %d %s\n', jj, z1(jj), cell2mat(rock(jj)) )      
    end
   
%    if  (strcmp(rtype(jj,1:8),'Rubble ')==0), 
	      zrock(ii)=zrock(ii)+zinc*dperc(jj);
%    end

   if (zrock(ii) > dz), fprintf('too much in window:  %f  %f  %f\n', ...
	       z1(jj),znext, zrock(ii)); 
   end

   if (jj==1),
     if ( (z2(jj) >= znext) || (z1(jj+1) > znext) ),
% set/initialize for next window down (ii)
	  zold = zold + rstep;  
	  jj=jct;
	  ii=ii+1;
      zfind = zold + rstep;
      znext = zold + dz;
      zmid(ii) = zold + dz/2;
     else
      jj=jj+1;
     end

   elseif (jj==nrec),
% set/initialize for next window down (ii)
    zold = zold + rstep;  
	jj=jct;
	ii=ii+1;
    zfind = zold + rstep;
    znext = zold + dz;
    zmid(ii) = zold + dz/2;
   
   elseif ( dperc(jj)==1),
     if ( (z2(jj) >= znext) || (z1(jj+1) > znext) ),
% set/initialize for next window down (ii)
	  zold = zold + rstep;  
	  ii=ii+1;
      zfind = zold + rstep;
      znext = zold + dz;
      zmid(ii) = zold + dz/2;
      jj=jct;
     else
% proceed summing subsections within this window
	   jj=jj+1;
     end
   else
% multi-domain interval with % specified
	 if (z1(jj)==z1(jj-1)),
	  if ( (z2(jj) >=znext) || (z1(jj+1) > znext) ),
       	zold = zold + rstep;  
	    ii=ii+1;
        zfind = zold + rstep;
        znext = zold + dz;
        zmid(ii) = zold + dz/2;
        jj=jct;
      else
        jj = jj+1;
      end
     else
      jj = jj+1;
     end

   end
end

%------------------ save results in table
c(nint,nk+1)=0;
c(:,1) = -zmid(1:nint);

if (iave == 1), div = dz;  end

for ii=1:nint,
  if (iave == 2), div = zrock(ii);  end
  for kk=1:nk,
    c(ii,kk+1)=rt(ii,kk)/div;
  end
end

save(fout,'c','-ascii');

%------------------ plot
       
l=length(c(:,1))+3;

cc(l,nk) = 0.;

% hardwired for 14 categories here...
cc(2:l-2,1)=c(:,2);
cc(2:l-2,2)=c(:,2)+c(:,3);
cc(2:l-2,3)=c(:,2)+c(:,3)+c(:,4);
cc(2:l-2,4)=c(:,2)+c(:,3)+c(:,4)+c(:,5);
cc(2:l-2,5)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6);
cc(2:l-2,6)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7);
cc(2:l-2,7)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8);
cc(2:l-2,8)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9);
cc(2:l-2,9)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9)+c(:,10);
cc(2:l-2,10)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9) ...
  +c(:,10)+c(:,11);
cc(2:l-2,11)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9) ...
  +c(:,10)+c(:,11)+c(:,12);
cc(2:l-2,12)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9) ...
  +c(:,10)+c(:,11)+c(:,12)+c(:,13);
cc(2:l-2,13)=c(:,2)+c(:,3)+c(:,4)+c(:,5)+c(:,6)+c(:,7)+c(:,8)+c(:,9) ...
  +c(:,10)+c(:,11)+c(:,12)+c(:,13)+c(:,14);


% depth vector; close it for fill region
z(2:l-2)=c(:,1);
zmn=max(z(2:l-2));
zmx=min(z);
z(1) = zmn; z(l-1)=zmx; z(l)=zmn;

% color definition for each category (RGB)
%peridotite
cm(:,1)=[25 115 75]/255;
%dunite
cm(:,2)=[75 185 0]/255;
%troctolite
cm(:,3)=[157 63 255]/255;
%troctolitic gabbro
cm(:,4)=[0 12 255]/255;
%ol gabbro
cm(:,5)=[9 162 224]/255;
%gabbro
cm(:,6)=[191 226 235]/255;
%gabbro norite
cm(:,7)=[255 255 182]/255;
%oxide gabbronorite
cm(:,8)=[255 255 0]/255;
%disseminated oxide gabbro
cm(:,9)=[238 189 8]/255;
%oxide bearing gabbro
cm(:,10)=[255 126 51]/255;
%oxide gabbro
cm(:,11)=[245 28 20]/255;
%diabase
cm(:,12)=[90 90 90]/255;
%felsic
cm(:,13)=[255 217 255]/255;

subplot(187)
for ll=nk:-1:1,
  cvec = cc(:,ll)';
  cmvec=cm(:,ll)';
  fill(cvec,z,cmvec);
  if (ll ==nk), hold on; end
end

axis([0 1 zmx  zmn])
%title('1473A 26Jan2016')
%ylabel('Depth  (mbsf), 20m running average')
xlabel(' Drilled Fraction')
%xlabel(' Recovered Fraction')
