clear all, clf

%set parameters
% nwr: # whole round records
% nsh: # shmsl records
% gnul: value (threshold?) of profile record when no core
% stpsize: value used for step detection, gradient of profile measure 
% short: minimum piece length to keep MS values
% pshort: minimum piece length to keep pt MS values
% eoff: offset from piece edge where data becomes good (cm)
% peoff: offset from piece edge for valid pt MS & GRA measurement
% dmin: minimum acceptable density (g/cm3)

gnul=2.1;
stpsize=0.05;
short=8;
pshort=2;
eoff=4;
peoff=1;
dmin=2.1;

%load data files (temporary for tst)
prof=csvread('prof.csv');
msl=load('msl.txt');
msp=load('msp.txt');
din=load('d.txt');
zp=prof(:,1); p=prof(:,2);
zm=msl(:,1); ml=msl(:,3); czm=msl(:,2);
zmp=msp(:,1); mp=msp(:,3);czp=msp(:,2);
zd=din(:,1); d=din(:,3); czd=din(:,2);

fin=input('section: ','s');

% length of each data series
nwr=length(ml);
nsh=length(p);
nshp=length(mp);
nd=length(d);

%find laser-detected gaps
k=0;
igap=0;

for i=1:nsh,
 if (p(i) <= gnul) & (igap==0),
  k=k+1;
  gs(k)=zp(i);
  igap=1;
 end
 if (p(i) > gnul) & (igap==1)
  ge(k)=zp(i-1);
  igap=0;
 end
end

if (igap==1), ge(k)=zp(nsh); end
ng1=k;

% replace short pieces with gap & move in from most piece edges

iflag=0;
for j=2:ng1,
  dzP = gs(j) - ge(j-1);
  if (dzP < short),
    ge(j-1)= gs(j);
    if (dzP < pshort),
      gpe(j-1)= gs(j);
      gps(j)=gs(j);
    else
      gpe(j-1)= gs(j)+peoff;
      gps(j)=gs(j)-peoff;
    end
    iflag=1;
  else
    if (iflag==0),
      gpe(j-1)=ge(j-1)+peoff;
      gps(j) = gs(j)-peoff;
      ge(j-1)=ge(j-1)+eoff;
      gs(j) = gs(j)-eoff;
    else
      gpe(j-1)=ge(j-1)+peoff;
      gps(j) = gs(j);
      ge(j-1)=ge(j-1)+eoff;
      gs(j) = gs(j);
      iflag=0;
    end
  end
end
if (gs(1)>1), 
    gps(1)=gs(1)-peoff; 
    gs(1)=gs(1)-eoff; 
end
if (ge(ng1) ~= nsh), 
    gpe(ng1)=ge(ng1)+peoff; 
    ge(ng1)=ge(ng1)+eoff; 
end

% insert dummy value into MSL for gap/short piece
% move past top edge effect

j=1;
for i=1:nwr,
  if ( (gs(j) <= zm(i) ) & ( zm(i) <= ge(j)) ),
    sm(i) = NaN;
  elseif (zm(i)<=eoff),
    sm(i) = NaN;
  else
    sm(i) = ml(i);
  end
  while ( (zm(i) >= ge(j)) & (j<ng1) ),
      j=j+1;
  end
end

% insert dummy value into Point MS for gap/short piece
% move past top edge effect

j=1;
for i=1:nshp,
  if ( (gps(j) <= zmp(i) ) & ( zmp(i) <= gpe(j)) ),
    smp(i) = NaN;
  elseif (zmp(i)<=peoff),
    smp(i) = NaN;
  else
    smp(i) = mp(i);
  end
  while ( (zmp(i) >= gpe(j)) & (j<ng1) ),
      j=j+1;
  end
end

% insert dummy value into GRA for gap/short piece
% move past top edge effect

j=1;
for i=1:nd,
  if ( (gps(j) <= zd(i) ) & ( zd(i) <= gpe(j)) ),
    sd(i) = NaN;
  elseif ((zd(i)<=peoff) || (d(i)<dmin) ),
    sd(i) = NaN;
  else
    sd(i) = d(i);
  end
  while ( (zd(i) >= gpe(j)) & (j<ng1) ),
      j=j+1;
  end
end


% findpiece-edge/cracks not detected, set pts to each side to NaN

k=0;
dp=diff(p);

for i=1:nsh-1,
 if (abs(dp(i)) > stpsize ),
  k=k+1;
  ss(k)=zp(i);
 end
end
nstp=k;

j=1;
% search within smaller offset for WRMS then NaN next adjacent pt
% since larger offset gets lost in counting steps for WRMS
for i=1:nwr,
  if ( ( (ss(j)-peoff) <= zm(i) ) & ( zm(i) <= (ss(j)+peoff) ) ),
    sm(i) = NaN;
    if (i< nwr-2) 
      sm(i+1) = NaN; 
      s(i+2) = NaN;
    end
  end
   while ( (zm(i) > ss(j)+ peoff ) & (j<nstp) ),
      j=j+1;
   end
end

j=1;
for i=1:nshp,
  if ( ( (ss(j)-peoff) <= zmp(i) ) & ( zmp(i) <= (ss(j)+peoff) ) ),
    smp(i) = NaN;
  end
   while ( (zmp(i) > ss(j)+ peoff ) & (j<nstp) ),
      j=j+1;
   end
end

j=1;
for i=1:nd,
  if ( ( (ss(j)-peoff) <= zd(i) ) & ( zd(i) <= (ss(j)+peoff) ) ),
    sd(i) = NaN;
      elseif ((zd(i)<=peoff) || (d(i)<dmin) ),
    sd(i) = NaN;

  end
   while ( (zd(i) > ss(j)+ peoff ) & (j<nstp) ),
      j=j+1;
   end
end

mmax=max(abs(ml));
mpmax=max(abs(mp));

plot(zp,p,'k'), hold on; title(fin(1:17))
plot(zm,2*ml/mmax,'m'), plot(zmp,2*mp/mpmax,'c'),plot(zd,d,'k:')
hold on, plot(zm,2*sm/mmax,'r+'),plot(zmp,2*smp/mpmax,'b+')
grid minor, plot(zp(2:nsh),10*dp+2,'g'),plot(zd,sd,'ko')
axis([0 150 -0.5 4.5])
%plot(ss,2,'mo') 
%plot(gs,2.2,'cd'), plot(ge,2.2,'c+')
ax=gca;
ax.XTick = [10 20 30 40 50 60 70 80 90 100 110 120 130 140 150];
xlabel('section depth (cm)'), ylabel('profile (cm), gradient (offset) & scaled MS')
text(130,4,'MSL: red'), text(130,3.7,'pt MS: blue')

filtms(:,1)=zm; filtms(:,2)=sm; filtms(:,3)=czm;
filtmsp(:,1)=zmp; filtmsp(:,2)=smp;filtmsp(:,3)=czp;
filtd(:,1)=zd; filtd(:,2)=sd; filtd(:,3)=czd;

%fout1(1:5)='fdat/';fout1(6:8)='fMS'; fout1(9:23)=fin(1:15);fout1(24:27)='.txt';
%fout2(1:5)='fdat/';fout2(6:9)='fpMS'; fout2(10:24)=fin(1:15);fout2(25:28)='.txt';
pfout(1:5)='fplt/';pfout(6:21)=fin(1:16); pfout(22:25)='.png';
%save(fout1,'filtms','-ascii')
%save(fout2,'filtmsp','-ascii')
save('fMS_1473n.txt','filtms','-ascii','-append')
save('fMSp_1473n.txt','filtmsp','-ascii','-append')
save('fd_1473n.txt','filtd','-ascii','-append')
print(pfout,'-dpng')

