clear all
close all

load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\0005\APEX\float0005_dilate.mat'
% float0005_dilate contains profile number, date, lat, long, P, theta, S, sigma,
    % B_LSS_b_corr, Chl_b_corr, O2sat, O2, POC_b, C_phy, and profile (in
    % days from start)
    
binz0 = [-1 30 75 130 185 245 315 480 750 1050]; % bins defined
nbinz = length(binz0)-1; % number of bins

pr = unique(float0005(:,1)); % get unique profile numbers 1-2XX
float0005_bin = NaN( length(pr)*nbinz , size(float0005,2)); % declare array float0005_bin with NaNs
B=NaN(length(pr)*nbinz,1); % declare array B to number of unique profiles*number of bins with NaNs

for kpr = 1:length(pr) % loop through each unique profile
    thispr = float0005(float0005(:,1)==pr(kpr),:);  % where col 1 in float0005 (profile number) == unique profile each loop
    for kbin = 1:nbinz % loop through each bin
        float0005_bin((kpr-1)*nbinz+kbin, [1:5]) = [thispr(1,1:4) binz0(kbin)]; % for each bin row 1, set first 4 cols = first 4 cols of each profile and 5th col = bin depth
        binidx = find(thispr(:,5)>=binz0(kbin) & thispr(:,5)<binz0(kbin+1)); % get all depths greater than bin1 and less than next bin
        B((kpr-1)*nbinz+kbin)=length(binidx); % get number of measurements in the bin 
        if B((kpr-1)*nbinz+kbin)==1 % if there's only one
            float0005_bin((kpr-1)*nbinz+kbin, 6:12) = thispr(binidx,6:12); % , use it not the median (otherwise MatLab will median the row and not the col)
        else % else there's more than one measurement
            float0005_bin((kpr-1)*nbinz+kbin, 6:12) = nanmedian(thispr(binidx,6:12)); % get median of all measurements found for each of the 7 cols and store in float0005_bin
        end
    end
end

% correlations
binz = binz0(1:(end-1));
cols = 6:12;
corrs = NaN(nbinz, length(cols));
robustcorr = NaN(nbinz, length(cols));

for kbin = 1:nbinz % for each bin, allcol

    surfbin = float0005_bin(float0005_bin(:,5)==binz(1), cols); % get surface binned median data
    thisbin = float0005_bin(float0005_bin(:,5)==binz(kbin), cols); % get each of the next 8 binned median data sets
       
    for kc = 1:length(cols)
        tmp = corrcoef(surfbin(:,kc),thisbin(:,kc), 'rows','complete'); % correlate surface bins with each of the next 8 bins
        corrs(kbin,kc) = tmp(1,2) % store in corrs
        
        [c,stats] = robustfit(surfbin(:,kc),thisbin(:,kc)); % correlate surface bins with each of the next 8 bins via robustfit
        robustcorr(kbin,kc) = stats.coeffcorr(1,2) % get coeffcorr from stats structure and store in robustcorr
       % clear stats
    end
end

% theta=corrs(:,1);
% sal=corrs(:,2);
% sigma=corrs(:,3);
% bb=corrs(:,4);
% chl=corrs(:,5);
% 
% thetaRobust=robustcorr(:,1);
% salRobust=robustcorr(:,2);
% sigmaRobust=robustcorr(:,3);
% bbRobust=robustcorr(:,4);
% chlRobust=robustcorr(:,5);

% figure, hold on
%     plot(corrs(:,1), binz, 'bo-')
%     plot(corrs(:,4), binz, 'rd-')
%     set(gca, 'ydir','reverse')
%     legend('TempC', 'bb')
%     
% figure, hold on
%     plot(corrs(:,2), binz, 'bo-')
%     plot(corrs(:,5), binz, 'rd-')
%     set(gca, 'ydir','reverse')
%     legend('Salinity', 'Chlorophyll')
%     
% % figure, hold on
%     plot(salRobust, binz, 'bo-')
%     plot(chlRobust, binz, 'rd-')
%     set(gca, 'ydir','reverse')
%     legend('Salinity', 'Chlorophyll')
    
% corrcoef(theta, binz) % temperature -0.6937
% [c,stats] = robustfit(theta,binz) % -0.9834
% stats.coeffcorr(1,2)
% 
% corrcoef(sal, binz) % salinity -0.6234
% [c,stats] = robustfit(sal,binz) % -0.99634
% 
% corrcoef(sigma, binz) % density -0.6192
% [c,stats] = robustfit(sigma,binz) % -0.58632
% 
% corrcoef(bb, binz) % bb -0.7375
% [c,stats] = robustfit(bb,binz) % -0.70775
% 
% corrcoef(chl, binz) % chl 0.0187
% [c,stats] = robustfit(chl,binz) % -0.14233

for kbin = 1:nbinz % for each bin, allcol
    Depth_bin_bb(:,kbin) = float0005_bin(float0005_bin(:,5)==binz(kbin), 9); 
    Depth_bin_chl(:,kbin) = float0005_bin(float0005_bin(:,5)==binz(kbin), 10); 
    Depth_bin_Sigma(:,kbin) = float0005_bin(float0005_bin(:,5)==binz(kbin), 8); 
    Depth_bin_T(:,kbin) = float0005_bin(float0005_bin(:,5)==binz(kbin), 6); 
end

time=unique(tt);

figure('Units','inches','Resize','off','position',[0 0 7.25 5.25],'paperposition',[0.5 3 7.25 5.25]);
subplot(2,2,1); % chl
semilogy(time,Depth_bin_chl(:,[1:2:9])+0.004,'LineWidth',2)
hold on
plot([732910 732910],[0.003 3],'k','LineWidth',2)
plot([732985 732985],[0.003 3],'k','LineWidth',2)
hold off
xlim([min(time),max(time)]);
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',10);
set(gca,'LineWidth',0.75);
ylabel('Chlorophyll [\mu gr/l]')
datetick(gca,'x','mmmyy','keepticks');
set(gca,'XTick',[]);
set(gca,'YTick',[0.003 0.03 0.3 3]);
axis([min(time) max(time) 0.003 3]);

subplot(2,2,2); % sigma
plot(time,Depth_bin_Sigma(:,[1:2:9]),'LineWidth',2)
hold on
plot([732910 732910],[24.5 28],'k','LineWidth',2)
plot([732985 732985],[24.5 28],'k','LineWidth',2)
hold off
xlim([min(time),max(time)]);
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',10);
set(gca,'LineWidth',0.75);
ylabel('\sigma_\theta [kg m^{-3}]')
datetick(gca,'x','mmmyy','keepticks');
set(gca,'XTick',[]);
axis([min(time) max(time) 24.5 28]);

subplot(2,2,3); % bb
semilogy(time,Depth_bin_bb(:,[1:2:9]),'LineWidth',2)
hold on
plot([732910 732910],[0.0001 0.01],'k','LineWidth',2)
plot([732985 732985],[0.0001 0.01],'k','LineWidth',2)
hold off
xlim([min(time),max(time)]);
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',10);
set(gca,'LineWidth',0.75);
set(gca,'XTick',[732115 732480 732850 733215]);
datetick(gca,'x','mmmyy','keepticks');
ylabel('Backscattering [m^{-1}]')
set(gca,'YTick',[0.0001 0.001 0.01]);
axis([min(time) max(time) 0.0001 0.01]);

subplot(2,2,4); % temp
plot(time,Depth_bin_T(:,[1:2:9]),'LineWidth',2)
hold on
plot([732910 732910],[3 20],'k','LineWidth',2)
plot([732985 732985],[3 20],'k','LineWidth',2)
hold off
xlim([min(time),max(time)]);
set(gca,'FontName','Times New Roman');
set(gca,'FontSize',10);
set(gca,'LineWidth',0.75);
set(gca,'XTick',[732115 732480 732850 733215]);
datetick(gca,'x','mmmyy','keepticks');
ylabel('Temperature [C]')
axis([min(time) max(time) 3 20]);
legend('0-30m','75-130m','185-245m','315-480m','750-1050m','location','northwest');
legend('boxon');
