% reads lithologic unit file as .csv format
% >>> need to set correct columns to depths
% z1=top of interval, z2=bottom of interval (mbsf)


%--------------------- 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='735BlithAllopx.csv';
fout='735BlithAllopx_RecAve10m.txt';
dz=10;
rstep=1;
iave=2;

%--------------------- data read
% read in subtable with depth, domain percent, and lithology
% z1 is top of subsection (mbsf)
% z2 is bottom of subsection (mbsf)
% rock is lithology name

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

nrec=length(z1);

%-------------------- define/group categories
% number of categories=nk
nk = 13;

%------------------- 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+nadd;
zold = zmin;
zfind = zold + rstep;
znext = zmin + dz;

%------------ 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) ),
%      if (((z1(jj) < zfind) && (z2(jj) < zfind))==1 ),
      while (((z1(jj) < zold) && (z2(jj) < zold))==1 ),
         jj=jj+1; 
      end

        if ( (z1(jj-1) <= zfind) && (z1(jj) > zfind) ),
	      if (z2(jj-1) > zfind),
  	        jct = jj-1;
	      else
	        jct = jj;
	      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(rock(jj) , 'Felsic')==1),
       rt(ii,13) = rt(ii,13) + zinc;
     
    elseif ( (strcmp(rock(jj), 'Diabase')==1) || ...
            (strcmp(rock(jj),'Basalt')==1 ) ),
       rt(ii,12) = rt(ii,12) + zinc;

    elseif (  ( strcmp(rock(jj), 'XGabbro')==1 ) || ...
            ( strcmp(rock(jj),'aXGabbro')==1 )  ),
       rt(ii,11) = rt(ii,11) + zinc;

    elseif (strcmp(rock(jj), 'xbGabbro')==1 ),
       rt(ii,10) = rt(ii,10) + zinc;

    elseif (strcmp(rock(jj), 'dxGabbro' )==1),
       rt(ii,9) = rt(ii,9) + zinc;

    elseif (strcmp(rock(jj), 'XNGabbro')==1),
       rt(ii,8) = rt(ii,8) + zinc;

    elseif (strcmp(rock(jj), 'NGabbro')==1),
       rt(ii,7) = rt(ii,7) + zinc;

    elseif ( (strcmp(rock(jj), 'Gabbro')==1) || ...
            (strcmp(rock(jj),'OBGabbro')==1) || ...
             (strcmp(rock(jj),'aGabbro')==1) || ...
         (strcmp(rock(jj),'PGabbro')==1) || ...
           (strcmp(rock(jj),'MGabbro')==1) ),
       rt(ii,6) = rt(ii,6) + zinc;

    elseif (strcmp(rock(jj), 'OGabbro')==1),
       rt(ii,5) = rt(ii,5) + zinc;
       
    elseif (strcmp(rock(jj), 'TGabbro')==1),
       rt(ii,4) = rt(ii,4) + zinc;

    elseif (strcmp(rock(jj), 'Troctolite')==1),
       rt(ii,3) = rt(ii,3) + zinc;
       
    elseif (strcmp(rock(jj), 'Dunite' )==1),
       rt(ii,2) = rt(ii,2) + zinc;

    elseif (strcmp(rock(jj), 'Perid' )==1),
       rt(ii,1) = rt(ii,1) + zinc;

	else
	  fprintf('Unclassified: %d %d %s\n', jj, z1(jj), cell2mat(rock(jj)) )      
    end
   
    if  (strcmp(rock(jj),'Rubble ')==0), 
	      zrock(ii)=zrock(ii)+zinc;
    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;
   
   else,
     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
   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(144)
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('735B')
ylabel('Depth  (mbsf), 10 m running average')
%xlabel(' Drilled Fraction')
xlabel(' Recovered Fraction')
