%bracketing_crosscheck.m
%Edwin Kite

%Purpose: Compare meander widths and wavelengths when these have been
%measured for the same channel. Based on results reported for Earth, the
%width/wavelength ratio shouldn't change dramatically with elevation and
%the ratio should be about 15 (depends on slope, grainsize, e.t.c.)

%Desired changes: Make everything part of a 'bc' (bracketing crosscheck) structure. 
savedebugfigures = 1;

WhichPass = 3;
load_channel_widths

%Identify the wavelength associated with width measurements that "bracketed" a meander
load parseTransect1ChannelWidthMeasurementsSecondPassForLPSC2014_v2_output_terse.mat
load parseTransect1MeanderWavelengthMeasurements_output_terse.mat

%Align coordinate systems
assert(sum(d(1).origin < origin) == 2)
bc.px = px - (d(1).origin(1) - origin(1)); %Points along banks (for width measurements)
bc.py = py - (d(1).origin(2) - origin(2));
%cx = cx - (d(1).origin(2) - origin(2)); %Bank-pair average positions
%cy = cy - (d(1).origin(2) - origin(2)); 


%For each channel width, determine if there is a meander-wavelength point
%inside the convex hull. Only bother with first pick:

    d(1).bracketed = zeros(size(d(1).vx));
    bc.widthlambda    = zeros(size(d(1).vx));
    bc.widthlambdaz   = zeros(size(d(1).vx));
    bc.lambda(1:(max(bid/2))) = 0;
   bc.zlambda(1:(max(bid/2))) = 0;
 bc.lambdaerr(1:(max(bid/2))) = 0;
      bc.acid(1:(max(bid/2))) = 0;
for nn = 1:(max(bid)/2) %Channel consisting of two banks
    thesebanks = find(bid==((nn*2) ) | bid==((nn*2)-1));
        bc.w(nn) = s(11).cw(nn); %Using the "uniform dotting" width measurement strategy
       bc.cx(nn) = nanmedian(px(bid==(nn*2) | bid==((nn*2)-1)));
       bc.cy(nn) = nanmedian(py(bid==(nn*2) | bid==((nn*2)-1)));
       
    if isempty(thesebanks)==0
        for m = 1:max(size(d(1).p))

        this_convex_hull = convhull(bc.px(thesebanks),bc.py(thesebanks));
        IN = inpolygon(d(1).p(m).x,d(1).p(m).y, ...
             px(thesebanks(this_convex_hull)), ...
             py(thesebanks(this_convex_hull)));
               d(1).p(m).bracketed = IN;
        if sum(IN)>0
            if median(d(1).p(m).arhm(IN)) > -1
                  bc.ahm(nn) = median(d(1).p(m).arhm(IN)); %associated half-meander
                 bc.acid(nn) = m; %associated channel ID
               bc.lambda(nn) = rhm(   median(d(1).p(m).arhm(IN))).lambda;
            bc.lambdaerr(nn) = rhm(   median(d(1).p(m).arhm(IN))).err;
              bc.zlambda(nn) = rhm(   median(d(1).p(m).arhm(IN))).z;
               bc.widthlambda(thesebanks) = rhm(median(d(1).p(m).arhm(IN))).lambda;
              bc.widthlambdaz(thesebanks) = rhm(median(d(1).p(m).arhm(IN))).z;
            end
        end
        
        end
    end
end

ww(142) = NaN; 
disp('Manually excluding bank-pair ID 142, from channel-ID 29! See notes from 21FEB2014, Rivers notebook')

%Check plot

hold on
for m = 1:max(size(d(1).p))
%figure
    scatter(d(1).p(m).x,d(1).p(m).y,'b.')
    axis equal; set(gca,'FontSize',18); hold on
text(mean(d(1).p(m).x),mean(d(1).p(m).y),num2str(m),'FontWeight','bold','FontSize',18)
xlimits = get(gca,'xlim');
ylimits = get(gca,'ylim');
scatter(px,py,'k.');
%set(gca,'xlim',xlimits)
%set(gca,'ylim',ylimits)
box on
scatter(bc.px(bc.widthlambdaz~=0),bc.py(bc.widthlambdaz~=0),30,'k','Filled'),%widthlambda(widthlambdaz~=0)
scatter(bc.cx(bc.lambda>0),bc.cy(bc.lambda>0),100,bc.lambda(bc.lambda>0).*2./bc.w(bc.lambda>0),'Filled')
    colorbar
    for n = 1:length(bc.cx)
        if bc.acid(n) == m
     text(bc.cx(n),bc.cy(n),{strcat('bankpair ID:',num2str(n)); ...
                       strcat('width(m):',num2str(bc.w(n))); ...
                       strcat('widtherr(m):',num2str(s(11).cwdev(n))); ...
                       strcat('lambda(m):',num2str(bc.lambda(n))); ...
                       strcat('lambdaerr(m):',  num2str(bc.lambdaerr(n)))}, ...
                       'FontAngle','oblique','FontSize',12,'Color','r');
        end
    end
    xlabel('distance x (m)'); ylabel('distance y (m)');
    title(strcat('Bracketing Crosscheck: Channel Centerline ID:',num2str(m)))
    if savedebugfigures == 1 
            %save to eps color file
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheck',num2str(m),'.eps'), ...
            'epsc');
    end
end




[b bint] = regress(bc.lambda(bc.lambda>0)',bc.w(bc.lambda>0)')

figure
scatter(bc.w(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)',100,bc.zlambda(bc.lambda>0)','Filled')
set(gca,'FontSize',18)
grid on; colorbar
text(2.5+bc.w(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
xlabel('width (m)');ylabel('half-meander lambda * 2'); title('color is elevation (m)'); box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckA.eps'), ...
            'epsc');

figure;
errorbar(bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',2.*bc.lambdaerr(bc.lambda>0)'./bc.w(bc.lambda>0)','.')
hold on; colorbar; set(gca,'FontSize',18)
scatter(bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',100,bc.w(bc.lambda>0)','Filled')
text(10+bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
set(gca,'FontSize',18)
xlim([-2400 -2050])
grid on
xlabel('elevation (m)');ylabel('half-meander lambda * 2 / width'); title('color is width (m)'); box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckB.eps'), ...
            'epsc');

figure;scatter(bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',100, ...
    log(cpstyle(bc.lambda>0)')./log(2),'Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
xlim([-2400 -2050])
xlabel('elevation (m)');ylabel('half-meander lambda * 2 / width'); title('color is log2(preservation style)');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckC.eps'), ...
            'epsc');

figure;scatter(bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',100,cqual(bc.lambda>0)','Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(10+bc.zlambda(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
xlim([-2400 -2050])
xlabel('elevation (m)');ylabel('half-meander lambda * 2 / width'); title('color is cpqual');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckD.eps'), ...
            'epsc');


figure;scatter(cqual(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',100,bc.w(bc.lambda>0)','Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(0.2+cqual(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
xlim([1 9])
xlabel('preservation quality');ylabel('half-meander lambda * 2 / width'); title('color is width (m)');box on
            saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckE.eps'), ...
            'epsc');

figure;scatter(cqual(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',100,log(cpstyle(bc.lambda>0))'./log(2),'Filled')
grid on; colorbar; set(gca,'FontSize',18)
text(0.2+cqual(bc.lambda>0)',2.*bc.lambda(bc.lambda>0)'./bc.w(bc.lambda>0)',num2str(bc.acid(bc.lambda>0)'))
xlabel('preservation quality');ylabel('half-meander lambda * 2 / width'); title('color is log2(preservation style)');box on
xlim([1 9])
saveas(gca, ...
            strcat('DebugParseTransect1BracketingCrosscheckF.eps'), ...
            'epsc');
