%Edwin Kite
%Purpose: Take as input hand-picked bank pairs and output paleo-channel
%width information.

%************************SET CONTROL PARAMETERS *****************************************************************
WhichTransect = 1; %Which geologic transect number
distancecutoff = 500; %Distance, in m, above which line segments are no longer evaluated (for channel widths)
dottingstrategy = [1]; %if = 1, dot uniformly along channel banks; if = 0 (deprecated), find width from vertices.
dbd = 2.5; %distance between dots (m)
%************************SET GRAPHICS AND HOUSEKEEPING PARAMETERS ***********************************************
picturesque = 1;
picture_all_measurements = 0; %Time-intensive, used for debugging.
output_type = 'terse'; %'none', 'terse', or 'complete'
%************************IMPORT DATA FROM CSV FILES ***************************
load_channel_widths
%************************END OF IMPORT DATA FROM CSV FILES ***************************

assert(min(px)==0); assert(min(py)==0); assert(minx~=0); assert(miny~=0);

%check image:
figure
scatter(px(zc(:,2)~=0),py(zc(:,2)~=0),10,z(zc(:,2)~=0)-zc(zc(:,2)~=0,2));
    
%********* Pepper Points Uniformly Along Channel Centerlines **************
if sum(dottingstrategy==1)>=1
        ctl(1:max(bid)) = 0; %channel (actually bank) total length
        %Interpolate uniformly in s, dir coordinates
        for l = 1:(size(px,1)-1)
            if bid(l) == bid(l+1)%Is it from the same bank?
                lseg(l)     =   sqrt((px(l+1) - px(l)).^2 + ...
                                     (py(l+1) - py(l)).^2); %length of segment
                ctl(bid(l)) = ctl(bid(l)) + lseg(l); %bank total length
                lac(l+1)    = ctl(bid(l)); %length along bank for this vertex
            else %last point on channel trace
                %lac(l) = ctl(cid(l));
            end
        end
        
        assert(min(ctl(ctl>0))>dbd,'Error. Decrease distance between dots (dbd)')

        for m = 1:max(bid) 
            if sum(bid==m) >= 1
            %position vector of interpolated points
            p(m).x = interp1(lac(bid==m),px(bid==m),[0:dbd:ctl(m)]);
            p(m).y = interp1(lac(bid==m),py(bid==m),[0:dbd:ctl(m)]);
            p(m).z = interp1(lac(bid==m), z(bid==m),[0:dbd:ctl(m)]);
            p(m).bid(1:max(size(p(m).x)))=m;
            end
        end
end

%*********************************FIND CHANNEL WIDTHS ***************************
dcls_temp = NaN;

for strategy_index = 1:length(dottingstrategy)
    dotuniformly = dottingstrategy(strategy_index)
    
    if dotuniformly == 0
        nstartpoints = size(px,1);
    elseif dotuniformly == 1
        nstartpoints = max(size([p(:).bid]));
        udx = [p(:).x]; udy = [p(:).y]; udz = [p(:).z]; udbid = [p(:).bid]; %uniformly-dotted x, y, z, and bank-ID
    end
    
%if dotuniformly == 0
for l = 1:nstartpoints %for each vertex
    if rem(l,100) == 0;disp(strcat('dotting strategy:',num2str(dotuniformly), ...
                                   ',l:',num2str(l)));end
    %if bid(min(size(px,1),l+1)) == bid(l) && bid(max(1,l-1)) == bid(l) ...
    %   && l > 1 && l < size(px,1)   %omit endpoints of banks
    
    if l > 1 && l < nstartpoints   %include endpoints of banks
    %find distance and coordinates of closest line segment
        if dotuniformly == 0
            x3 = px(l);y3 = py(l);
        elseif dotuniformly ==1
            x3 = udx(l);y3 = udy(l);
        end
        
        pd = sqrt((x3 - px).^2 + (y3 - py).^2);
        %pd(1) = distancecutoff*9e9; 
        pd(end) = distancecutoff*9e9; %prevent evaluation of first and last m-points
        xc = NaN*zeros(size(px));
        yc = xc;
        
        okmlist = find(pd<distancecutoff);
        
    for mindex = 1:length(okmlist) %2:(size(px,1)-1) 
           m = okmlist(mindex);
           %if pd(m) < distancecutoff
           %per http://stackoverflow.com/questions/6176227/for-a-point-in-an-irregular-polygon-what-is-the-most-efficient-way-to-select-th
           %but with sqrt added to denominator of equation for u(m) 
            x1 = px(m)  ; y1 = py(m);
            x2 = px(m+1); y2 = py(m+1);
        
        %segment is finite
            dQ1(m) = sqrt((x1 - x3)^2 + (y1 - y3)^2); %<-- Bugs result if these lines are moved inside the "filter out line segment between different bank ID's" loop
            
            dQ2(m) = sqrt((x2 - x3)^2 + (y2 - y3)^2);
            
        %DO consider the first line segment (for which bid(m-1)~=bid(m))
        if dotuniformly == 0
            acceptable_endpoint = (bid(m+1) == bid(m) && l~=m && bid(m)~=bid(l));
        elseif dotuniformly ==1
            acceptable_endpoint = (bid(m+1) == bid(m) && bid(m)~=udbid(l));
        end
            
        if acceptable_endpoint == 1 
      
%Following http://paulbourke.net/geometry/pointlineplane/
%using lambda instead of u and PC instead of P
lambda(m) = ( (x3-x1).*(x2-x1) + (y3-y1).*(y2-y1) ) ./ ...
             ((x2-x1).^2 + (y2-y1).^2);


if lambda(m)<0
    xc(m) = x1; yc(m) = y1;
    lambda(m) = 0;
elseif lambda(m)>1
    xc(m) = x2; yc(m) = y2;
    lambda(m) = 1;
else
    xc(m) = x1 + lambda(m)*(x2-x1);
    yc(m) = y1 + lambda(m)*(y2-y1);
end
          
        else
            %skip. Zero values of xc will be used to set dcls_temp(m) to NaN
        end
            %skip. Zero values of xc will be used to set dcls_temp(m) to NaN
    end
  
 dcls_temp = sqrt((x3 - xc).^2 + (y3 - yc).^2);        
 dcls_temp(isnan(xc)) = NaN; % Zero values of xc are used to set dcls_temp(m) to NaN
    
    [C,I] = nanmin(dcls_temp);
    if dotuniformly == 0
        lambdaforl(l) = lambda(I);
        clx(l) = xc(I);
        cly(l) = yc(I);
        dcls(l) = C;
        cls(l) = I; %line segment number, not point number
        clsbank(l) = bid(I);
    elseif dotuniformly == 1
        udlambdaforl(l) = lambda(I);
        udclx(l) = xc(I);
        udcly(l) = yc(I);
        uddcls(l) = C;
        udcls(l) = I; %line segment number, not point number
        udclsbank(l) = bid(I);
    end
    
     %We want to suppress the case where the closest *point* (not line segment) is a
            %vertex at the END of a bank (banks are indexed by 'bid')
                dQdebug(l) = 0;
            if I == 1; %The closest point is the first point in the list of all points on all banks
                dQdebug(l) = 9; %Suppress the width measurement
            elseif dQ1(I) == dcls_temp(I) && bid(I-1)~=bid(I)
                dQdebug(l) = 1;
            elseif dQ2(I) == dcls_temp(I) && bid(min(length(bid),I+2))~=bid(I)
                dQdebug(l) = 2;
            end
            
            if dQdebug(l) > 0
                if dotuniformly == 0
                    dcls(l) = NaN; cls(l) = NaN; clsbank(l) = NaN;
                elseif dotuniformly == 1
                    udcls(l) = NaN; uddcls(l) = NaN; udclsbank(l) = NaN;
                end
            end
    
    else
        if dotuniformly == 0
            cls(l) = NaN; dcls(l) = NaN; clsbank(l) = NaN;
        else
            udcls(l) = NaN; uddcls(l) = NaN; udclsbank(l) = NaN;
        end
    end
end
end

%***********END OF FIND CHANNEL WIDTHS ***************************

if picturesque == 1
    figure %check plot
    scatter(px,py,'k*')
    hold on
    title('plot for debugging')
    EdwinColor = {'k','b','r','g','y',[.5 .6 .7],[.8 .2 .6]};
    if picture_all_measurements == 1
            for l = 2:(size(px,1)-1)
                if isnan(cls(l)) == 0
                    line([px(l) clx(l)],[py(l) cly(l)],'Color', EdwinColor{1+rem(l,6)})
                    %line([px(l) px(cls(l))],[py(l) py(cls(l))],'Color', EdwinColor{1+rem(l,6)})
                    scatter(px(l),py(l),20,EdwinColor{1+rem(l,6)},'LineWidth',2)
                end
            end
            %    scatter(px(1:length(dQdebug)),py(1:length(dQdebug)),1500,dQdebug);

    for l = 2:max(size(udclx))
        if isnan(udcls(l)) == 0
            line([udx(l) udclx(l)],[udy(l) udcly(l)],'Color', EdwinColor{1+rem(l,6)})
            %line([px(l) px(cls(l))],[py(l) py(cls(l))],'Color', EdwinColor{1+rem(l,6)})
            scatter(udx(l),udy(l),20,EdwinColor{1+rem(l,6)},'LineWidth',2)
        end
    end

          %     scatter(px(1:length(dQdebug)),py(1:length(dQdebug)),50,dQdebug);
    end
end 
    
%**********AGGREGATE POINT-BY-POINT MEASUREMENTS TO FIND CHANNEL FIRST MOMENTS ***************************

%now find median width, dz, and rootsemi for each channel
    %chann = floor((bid+2)/2); %channel number
    
    %Note 21FEB2014: This section should be cleaned up so that all
    %measurements go into the 's' (strategy) structure. Check lower down in
    %the script where some of the 'cx', 'cy', e.t.c. are used for plotting.
for nn = 1:(max(bid)/2) %each channel has two banks; they are assumed to be in order of this pairing!
    
    css(nn)  = sum(bid==((nn*2) ) | bid==((nn*2)-1));%channel sample size (for calculating pooled variance, later)
    if css(nn)>0
        cx(nn)     = nanmean(px(bid==(nn*2) | bid==((nn*2)-1)));
        cy(nn)     = nanmean(py(bid==(nn*2) | bid==((nn*2)-1)));
        cqual(nn)  = nanmean(pqual((bid==(nn*2) | bid==(nn*2)-1))); %note that FID == a(:,2)
        cz(nn)     = nanmean(z((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        if (WhichTransect == 1 || WhichTransect == 2);
        for iss = 1:5
            czs(nn,iss)= nanmean(z((bid==(nn*2) | bid==(nn*2)-1),iss) - zc((bid==(nn*2) | bid==(nn*2)-1),iss)); 
        end
        end 
        czdev(nn)      = nanstd(z((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        cpstyle(nn)    = nanmean(pstyle((bid==(nn*2) | bid==(nn*2)-1)));
        ctienumber(nn) = nanmean(tienumber((bid==(nn*2) | bid==(nn*2)-1)));
        
        
        if sum(dottingstrategy==1)>=1
            banktobankdistances = uddcls(udbid==(nn*2) | udbid==((nn*2)-1));
            s(11).cw(nn)   = nanmean(banktobankdistances);
            s(11).cwrange(nn) = max(banktobankdistances)./min(banktobankdistances);
            s(11).cwmin(nn) = nanmin(banktobankdistances);
            s(11).cwmax(nn) = nanmax(banktobankdistances);
            s(11).cwdev(nn) = nanstd(uddcls(udbid==(nn*2) | udbid==(nn*2)-1));
        end
                
        if exist('dz','var') == 1
            cdz(nn) = nanmean(dz((bid==(nn*2) | bid==(nn*2)-1) & isnan(dcls')==0));
        end
    else
        cz(nn) = NaN;
    end
end

%**********END OF AGGREGATE POINT-BY-POINT MEASUREMENTS TO FIND CHANNEL FIRST MOMENTS ***************************


%**************TO BE IMPLEMENTED: AGGREGATE RESULTS WITH THE SAME CHANNEL ID *******


%******************** Aggregate Replicable Half-Meanders on A Single Centerline Trace ******************************

%Combine the mean and standard deviation of bank-pairs with the same
%tienumber to obtain a characteristic width for each tienumber.
%A common tienumber indicates that the bank-pairs probably sit on the same
%channel thread.

clear  bpti t
ti = 0;

for m = 1:length(s(11).cw) %For each channel
    if cqual(m)<8 %Do not include candidates
        
    ctie = find(ctienumber == ctienumber(m));%bank-pairs with the same tienumber as this bank-pair
        
    if (length(ctie)==1 || ctienumber(m) ==0) %This bank-pair has a unique tienumber.
        ti = ti + 1; %tie index
        t(ti).cw    = s(11).cw(m);
        t(ti).cwdev = s(11).cwdev(m);
        t(ti).cx     = cx(m);
        t(ti).cy     = cy(m);
        t(ti).cz     = cz(m);
        t(ti).czdev  = czdev(m);
        bpti(m) = ti; %bank-pair tie index
            if (WhichTransect == 1 || WhichTransect == 2);
                t(ti).czc    = czc(m,:);
            end
        t(ti).intie = ctie;
    else %Combine the data for the bank-pairs with the same tienumber.
        if m == min(ctie); %only one t (tie) entry per tienumber
            
        ti = ti + 1; %tie index
        t(ti).cw    = nanmean(s(11).cw(ctie));
        t(ti).cx    = nanmean(cx(ctie));
        t(ti).cy    = nanmean(cy(ctie));
        t(ti).cz    = nanmean(cz(ctie));
        t(ti).intie = ctie;    
            if (WhichTransect == 1 || WhichTransect == 2);
                        t(ti).czc   = nanmean(czc(ctie,:));
            end
    %Based on Wikipedia article
    %'Standard_deviation#Combining_standard_deviations' as of 24FEB2014
    %Non-overlapping populations
    t(ti).cwdev = sqrt( ( sum(( [s(11).cwdev(ctie)].^2 + [s(11).cw(ctie)].^2 )) ...
                    / length(ctie))...
                    - t(ti).cw.^2);
    t(ti).czdev = sqrt( ( sum(( [czdev(ctie)].^2 + [cz(ctie)].^2 )) ...
                    / length(ctie))...
                    - t(ti).cz.^2);            
                
        bpti(ctie) = ti;
        end
    end
    else
        bpti(m) = NaN;
    end

end

%Debug plots
figure
    errorbar([1:length(t)],[t.cw],[t.cwdev],'r.')
    hold on
    errorbar(0.25+bpti,[s(11).cw],[s(11).cwdev],'*');
figure
    errorbar(0.25+bpti,[cz],[czdev],'m*');
        hold on
    errorbar([1:length(t)],[t.cz],[t.czdev],'g.')


%**************DISPLAY RESULTS*******

%Kludge fix for mixed preservation styles
cpstyle(cpstyle==9) = 8;
cpstyle(cpstyle==10) = 8;

figure('units','normalized','outerposition',[0 0 1 1])
    h = errorbar(cz,s(11).cw,s(11).cwdev,'xk','LineWidth',2)
    hold on;grid on; colorbar; errorbar_tick(h,0)
    scatter(cz,s(11).cw,30,cx,'LineWidth',2)
        scatter(cz(cqual>=8),s(11).cw(cqual>=8),300,'ro','LineWidth',1)
    title('Raw elevation, candidates circled in red'); xlabel('Raw elevation (m)';ylabel('Channel width (m)');
                saveas(gca,'DebugPlanarFit30MAR2014_1.eps','epsc');
    

figure('units','normalized','outerposition',[0 0 1 1])
    h = errorbar(cz(cqual<=4) - czc(cqual<=4,2)' ,s(11).cw(cqual<=4),s(11).cwdev(cqual<=4),'xk','LineWidth',1)
    hold on;colorbar;grid on; errorbar_tick(h,0)
    set(gca,'FontSize',18)
    title('Quadratic correction, color is log2(pstyle), no candidates')
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(cz(cqual<=4) - czc(cqual<=4,2)' ,s(11).cw(cqual<=4),30,log(cpstyle(cqual<=4))./log(2),'LineWidth',2)    
                    saveas(gca,'DebugPlanarFit30MAR2014_3.eps','epsc');

    
    
figure('units','normalized','outerposition',[0 0 1 1])
        for irhm = 1:length(t)
            tzadjusted(irhm) = [t(irhm).cz] - [t(irhm).czc(2)];
        end
    h = errorbar(tzadjusted,[t(:).cw],[t(:).cwdev],'xk','LineWidth',1)
    hold on;grid on;colorbar; errorbar_tick(h,0);
    set(gca,'FontSize',18)
    title({'Quadratic correction';'color is x position'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(tzadjusted,[t(:).cw],30,[t(:).cx],'LineWidth',2)
                    saveas(gca,'DebugPlanarFit30MAR2014_4.eps','epsc');                    
                    
                   
figure('units','normalized','outerposition',[0 0 1 1])
    h = errorbar(cz(cqual<=4) - czc(cqual<=4,2)', ...
             s(11).cw(cqual<=4),s(11).cwdev(cqual<=4),'xk','LineWidth',1)
    hold on;grid on;colorbar; errorbar_tick(h,0);
    czcorr = cz(cqual<=4) - czc(cqual<=4,2)';
    set(gca,'FontSize',18)
    title({'Quadratic correction';'color is log2(quality), no candidates'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not kridged stratigraphic position!) is usually negligible
    scatter(cz(cqual<=4) - czc(cqual<=4,2)', ...
            s(11).cw(cqual<=4),30,log(cqual(cqual<=4))./log(2),'LineWidth',2);%log(cqual(cqual<=4))./log(2),'LineWidth',2)
                    saveas(gca,'May2014DebugQuadraticFit30MAR2014_4_PQualExcludeCandidates.eps','epsc');
                    
                    
figure('units','normalized','outerposition',[0 0 1 1])
    h = errorbar(cz(cqual<=4) - czc(cqual<=4,2)', ...
             s(11).cw(cqual<=4),s(11).cwdev(cqual<=4),'xk','LineWidth',2)
    hold on;grid on;colorbar; errorbar_tick(h,0)
    czcorr = cz(cqual<=4) - czc(cqual<=4,2)';
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality), no candidates'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not kridged stratigraphic position!) is usually negligible
    scatter(cz(cqual<=4) - czc(cqual<=4,2)', ...
            s(11).cw(cqual<=4),30,rem(ctienumber(cqual<=4),7)','LineWidth',2);%log(cqual(cqual<=4))./log(2),'LineWidth',2)

                    saveas(gca,'May2014DebugQuadraticFit30MAR2014_5_PQualExcludeCandidates.eps','epsc');

                    
                                 
%******************************************For 8th Mars 2014*********************************                  
figure('units','normalized','outerposition',[0 0 0.35 1])
    errorbar(cz(cqual<=4) - (cetl + cx(cqual<=4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4).*pfyc(WhichTransect)), ...
             s(11).cw(cqual<=4),s(11).cwdev(cqual<=4),'xr','LineWidth',2)
    hold on;grid on;
    czcorr = cz(cqual<=4) - (cetl + cx(cqual<=4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4).*pfyc(WhichTransect));
    set(gca,'FontSize',18)
    %title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality), no candidates'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
  %  scatter(cz(cqual<=4) - cetl + cx(cqual<=4).*pfxc(WhichTransect) ...
  %                               + cy(cqual<=4).*pfyc(WhichTransect)), ...
  %          s(11).cw(cqual<=4),30,log(cqual(cqual<=4))./log(2),'LineWidth',2)
    wunit0(1,1) =    nanmean(s(11).cw(czstrat <(0))); %Mean channel width above and below paleohydrologic break
    wunit0(2,1) =    nanmean(s(11).cw(czstrat>=(0))); %Mean channel width above and below paleohydrologic break
    wunit0(1,2) =  nanmedian(s(11).cw(czstrat <(0))); %Median channel width above and below paleohydrologic break
    wunit0(2,2) =  nanmedian(s(11).cw(czstrat>=(0))); %Median channel width above and below paleohydrologic break
    wunit0(1,3) =     nanstd(s(11).cw(czstrat <(0))); %s.d. of channel width above and below paleohydrologic break
    wunit0(2,3) =     nanstd(s(11).cw(czstrat>=(0))); %s.d. of channel width above and below paleohydrologic break
    
  line([0 0],[0.01 300],'Color','k')
  line([-150 0],[wunit0(1,1) wunit0(1,1)],'Color','k')
  line([0 150],[wunit0(2,1) wunit0(2,1)],'Color','k')
view([-90, 90])
set(gca,'ydir','reverse')
grid off;% set(gca,'yminortick','on')
set(gca,'yscale','log')
xlim([-70 130])
ylim([3 130]) 
set(gca,'ytick',[5 10 15 20 30 40 60 80 120])

                    saveas(gca,'For8thMars2014_DebugPlanarFit30MAR2014_4_PQualExcludeCandidates.fig','fig')
                    saveas(gca,'For8thMars2014_DebugPlanarFit30MAR2014_4_PQualExcludeCandidates.eps','epsc');

 %Split by x = 1e4 (for dtm seperation on transect 1:)
 
 figure('units','normalized','outerposition',[0 0 0.35 1])
    errorbar(cz(cqual<=4 & cx<1e4) - (cetl + cx(cqual<=4 & cx<1e4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4 & cx<1e4).*pfyc(WhichTransect)), ...
             s(11).cw(cqual<=4 & cx<1e4),s(11).cwdev(cqual<=4 & cx<1e4),'xr','LineWidth',2)
    hold on;grid on;
    czcorr = cz(cqual<=4 & cx<1e4) - (cetl + cx(cqual<=4& cx<1e4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4 & cx<1e4).*pfyc(WhichTransect));
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality), no candidates'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
  %  scatter(cz(cqual<=4) - (cetl + cx(cqual<=4).*pfxc(WhichTransect) ...
  %                               + cy(cqual<=4).*pfyc(WhichTransect)), ...
  %          s(11).cw(cqual<=4),30,log(cqual(cqual<=4))./log(2),'LineWidth',2)
    wunit0(1,1) =    nanmean(s(11).cw(czstrat <(0) & cx<1e4)); %Mean channel width above and below paleohydrologic break
    wunit0(2,1) =    nanmean(s(11).cw(czstrat>=(0) & cx<1e4)); %Mean channel width above and below paleohydrologic break
    wunit0(1,2) =  nanmedian(s(11).cw(czstrat <(0) & cx<1e4)); %Median channel width above and below paleohydrologic break
    wunit0(2,2) =  nanmedian(s(11).cw(czstrat>=(0) & cx<1e4)); %Median channel width above and below paleohydrologic break
    wunit0(1,3) =     nanstd(s(11).cw(czstrat <(0) & cx<1e4)); %s.d. of channel width above and below paleohydrologic break
    wunit0(2,3) =     nanstd(s(11).cw(czstrat>=(0) & cx<1e4)); %s.d. of channel width above and below paleohydrologic break
    
  line([0 0],[0.01 300],'Color','k')
  line([-150 0],[wunit0(1,1) wunit0(1,1)],'Color','k')
  line([0 150],[wunit0(2,1) wunit0(2,1)],'Color','k')
view([-90, 90])
set(gca,'ydir','reverse')
grid off;% set(gca,'yminortick','on')
set(gca,'yscale','log')
xlim([-70 130])
ylim([3 130]) 
set(gca,'ytick',[5 10 15 20 30 40 60 80 120])

                    saveas(gca,'For8thMars2014_PQualExcludeCandidates_LowX.fig','fig')
                    saveas(gca,'For8thMars2014_PQualExcludeCandidates_LowX.eps','epsc');

%High split by x = 1e4

 figure('units','normalized','outerposition',[0 0 0.35 1])
    errorbar(cz(cqual<=4 & cx>1e4) - (cetl + cx(cqual<=4 & cx>1e4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4 & cx>1e4).*pfyc(WhichTransect)), ...
             s(11).cw(cqual<=4 & cx>1e4),s(11).cwdev(cqual<=4 & cx>1e4),'xr','LineWidth',2)
    hold on;grid on;
    czcorr = cz(cqual<=4 & cx>1e4) - (cetl + cx(cqual<=4& cx>1e4).*pfxc(WhichTransect) ...
                                  + cy(cqual<=4 & cx>1e4).*pfyc(WhichTransect));
    set(gca,'FontSize',18)
    title({'Planar correction, one point per pick on F1-F2 contact,';'color is log2(quality), no candidates'})
    xlabel('Stratigraphic elevation relative to F1-F2 contact (m)'); ylabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
  %  scatter(cz(cqual<=4) - (cetl + cx(cqual<=4).*pfxc(WhichTransect) ...
  %                               + cy(cqual<=4).*pfyc(WhichTransect)), ...
  %          s(11).cw(cqual<=4),30,log(cqual(cqual<=4))./log(2),'LineWidth',2)
    wunit0(1,1) =    nanmean(s(11).cw(czstrat <(0) & cx>1e4)); %Mean channel width above and below paleohydrologic break
    wunit0(2,1) =    nanmean(s(11).cw(czstrat>=(0) & cx>1e4)); %Mean channel width above and below paleohydrologic break
    wunit0(1,2) =  nanmedian(s(11).cw(czstrat <(0) & cx>1e4)); %Median channel width above and below paleohydrologic break
    wunit0(2,2) =  nanmedian(s(11).cw(czstrat>=(0) & cx>1e4)); %Median channel width above and below paleohydrologic break
    wunit0(1,3) =     nanstd(s(11).cw(czstrat <(0) & cx>1e4)); %s.d. of channel width above and below paleohydrologic break
    wunit0(2,3) =     nanstd(s(11).cw(czstrat>=(0) & cx>1e4)); %s.d. of channel width above and below paleohydrologic break
    
  line([0 0],[0.01 300],'Color','k')
  line([-150 0],[wunit0(1,1) wunit0(1,1)],'Color','k')
  line([0 150],[wunit0(2,1) wunit0(2,1)],'Color','k')
view([-90, 90])
set(gca,'ydir','reverse')
grid off;% set(gca,'yminortick','on')
set(gca,'yscale','log')
xlim([-70 130])
ylim([3 130]) 
set(gca,'ytick',[5 10 15 20 30 40 60 80 120])

                    saveas(gca,'For8thMars2014_PQualExcludeCandidates_HighX.fig','fig')
                    saveas(gca,'For8thMars2014_PQualExcludeCandidates_HighX.eps','epsc');
                    
                    
                    
%END OF FOR 8TH MARS 2014                    
                    
                    
figure
   
    for nn = 1:(max(bid)/2)
    pointwisezstrat = (zstrat((bid==(nn*2) | bid==(nn*2)-1)));
    plot(ones(size(pointwisezstrat))*(cz(nn) - (cetl + cx(nn).*pfxc(WhichTransect) ...
                                                     + cy(nn).*pfyc(WhichTransect))), ...
            pointwisezstrat);grid on; hold on
    end
        title('Check: should be a straight 1-to-1 line')
        ylabel('median elevation with stratigraphic correction');
        xlabel('pointwise median stratigraphic elevation');
   
        
figure
    title('Check')
    scatter(cx,cy,30,(cetl + cx.*pfxc(WhichTransect) ...
                           + cy.*pfyc(WhichTransect)),'*' ...
            );grid on;colorbar    
        hold on
    %scatter(cx,cy+100,30,-(czstrat - cz),'o' ...
    %        );grid on;colorbar    
    %    hold on       
    xlabel('x');ylabel('y')
figure

    errorbar(cdz,s(11).cw,cwdev,'xk','LineWidth',2)
    hold on
    scatter(cdz,s(11).cw,50,'co','fill')
    scatter(cdz(cdz<0),         s(11).cw(cdz<0),50,'bo','fill')
    scatter(cdz(cdz>0 & cdz<50),s(11).cw(cdz>0 & cdz<50),50,'wo','fill')
    scatter(cdz(cdz>0 & cdz<50),s(11).cw(cdz>0 & cdz<50),50,'ko')

    %errorbar(cdz(cwhichdtm==2),cw(cwhichdtm==2),cwdev(cwhichdtm==2),'xr')

    %herrorbar(cdz,cw,crootsemi,'x')
xlabel('Relative Time (Height Above F1/F2 Contact (m))','FontSize',20)
ylabel('River Channel Width (m)','FontSize',20)
set(gca,'FontSize',20)
set(gca,'xminortick','on')
set(gca,'yminortick','on')
ylim([0 70])
xlim([-100 150])
%grid on
box on
line([-100,0],[29.3,29.3],'Color','b','LineWidth',2)
line([50,150],[15.8616,15.8616],'Color','c','LineWidth',2)
line([0,50],[ 26.8001, 26.8001],'Color','k','LineWidth',2)

%********** SAVE OUTPUT ***************************
switch output_type
    case 'terse'
        widthoutputname = strcat('parseWidthMeasurementsMarsRivers_output_terse',num2str(WhichTransect))        
        save(widthoutputname, ...
             's','cqual','cpstyle','cx','cy','cz','minx','miny')
    case 'complete'
        disp('Not saving complete output.')
    case 'none'
end

%********** END OF SAVE OUTPUT ***************************

%error('Not doing probability summation while debugging')

%Quick summation of probability distribution (assuming uncorrelated, Guassian error
%bars) From Wikipedia, "Multivariate normal distribution"
p = zeros(601,401);
warning('Assuming fixed sigmax for debug')
for nn = 1:(max(bid)/2) %add the probability kernel of each observation
oldp = p;
nn
            sigmax = 29.0847;%<-- Assuming fixed sigma x for debugging %crootsemi(nn);
            %The value of fixed sigma x comes from
            %f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output
            
            sigmay = s(11).cwdev(nn);
            normfactor = 1/(2*pi*sigmax*sigmay);
            mux = cz(nn) - (cetl + cx(nn).*pfxc(WhichTransect) ...
                                 + cy(nn).*pfyc(WhichTransect));
            muy = s(11).cw(nn);
            
            if sum(isnan([sigmax,sigmay,mux,muy]))==0
                for x = -300:1:300 %height above F1/F2 contact (m)
                    for y = -100:1:300 %river channel width (m)

                            p(x+301,y+101) = p(x+301,y+101) + ...
                            normfactor*exp(-1/2*(((x-mux)^2/(sigmax^2))+((y-muy)^2/(sigmay^2))));
                    end
                end
            end
                %For debug only
                %check that the additional probability sums to 1 (note that
                %the discretisation and the cutoffs will introduce errors!)
                %if abs((sum(p(:)) - sum(oldp(:)))-1) > 0.01
                %    error('MVN Dist Checksum Error!')
                %end
 end

 [X,Y] = meshgrid(-300:1:300,-100:1:300);
%contour(X,Y,p');grid on %show the sum of all kernels

%show 1 sigma and 2 sigma error windows for channel width as a function of
%stratigraphic position

pcs = cumsum(p,2);
ps = sum(p');
for l = 1:601 %cumsum normalized
pcsn(l,:) = pcs(l,:)./ps(l);
end

hold on
contour(X,Y,pcsn',[0.05,(0.31731/2),1-((0.31731/2)),0.95]);

%find pooled variance for 50m stratigaphic-position bins: (Note to self:
%should preferably use *fractional* variance?)
stratbin = 1;
for stratz = -200:50:200
    to_use(stratbin).v = find(cdz>stratz & cdz <stratz+50);
    ninbin(stratbin) = length(to_use(stratbin).v);
    pooledvar(stratbin) = sqrt(sum(css(to_use(stratbin).v).*(cwdev(to_use(stratbin).v).^2))/sum(css(to_use(stratbin).v)));
    stratbin = stratbin + 1;
end

%Overall pooled variance:
opooledvar = sqrt(sum(css.*(cwdev.^2))/sum(css))