Plot bbp final

function [eta, etaf] = plot_bbp_final(cruise)

[filename,path] = uigetfile(['home/data65/pb/HYDROSCAT/',cruise,'/final/'],'Choose HS-6 file to plot');

fid = fopen([path filename],'r');
C = textscan(fid,'%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f','delimiter',',','HeaderLines',120);
% C = textscan(fid,'%s%s%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f','delimiter',',');
fclose(fid);

%lambda = [ 442 470 510 589 671];
%newlambda = [420 442 470 510 589 700] as of 4/2010
lambda = [420 442 470 510 589 700];
% lambda = [420 700 442 510 470 589];
bbs = [C{69} C{70} C{71} C{72} C{73} C{74}]; 
%%bbs = [C{69} C{71} C{73} C{72} C{74} C{70}];%re-arrange to match uncorrected 3/7/2011
bbz = [C{75} C{76} C{77} C{78} C{79} C{80}]; 
%%bbz = [C{75} C{77} C{79} C{78} C{80} C{76}];%re-arrange to match uncorrected 3/7/2011
bbs(bbs==-9.9e+35) = NaN;
bbz(bbz==-9.9e+35) = NaN;


%after Mobley 3.31
bbw_525 = (16.06/2)*1.46e-4;
bbw = bbw_525*((525./lambda).^4.32);

bbps = bbs - repmat(bbw,size(bbs,1),1);
%plot(1:length(bbps),bbps);
depths = C{81};

% s = input('Choose which scattering correction to apply, method 1, 2, or 3; Enter 0 for none?');
%
% switch s
%     case 0
%         q = ac9.out.a_corr;
%     case 1
%         q = ac9.out.a_scorr1;
%     case 2
%         q = ac9.out.a_scorr2;
%     case 3
%         q = ac9.out.a_scorr3;
%     otherwise
% end

%
% qc = ac9.out.c_corr;
% qc(qc==-9.9e+35) = NaN;
% b = qc-q;

plot_spectrum(bbps,'b_b_p',depths,filename,bbw);
plot_spectrum(bbz,'b_b\_Zhang',depths,filename,bbw);
%plot_spectrum(b,'b',filename);

function S = plot_spectrum(q,type,d,filename,bbw)
%wl = [ 442 470 510 589 671];
wl = [420 442 470 510 589 700];
% wl = [420 700 442 510 470 589];
% d = 1:220';
% d = d';

%Plot all spectra at once
%%q = q(:,[1 6 2 4 3 5]);
qd = denan(q);

c = jet(size(qd,1));


%%%Calculate eta, using 442, 510 and 589 nm wavelengths
eta = [];
opts = optimset('fminsearch');
opts  = optimset('MaxFunEvals', 5000,'MaxIter',5000, 'TolFun',1e-6,'TolX',1e-4);

lambda_for_eta = wl([1 3 4]);
%%%%%lambda_for_eta = wl([1 5 4]);
lambda_for_eta = lambda_for_eta';

for i = 1:size(q,1)
    bbp = q(i,[1 3 4])';
  %%%%  bbp = q(i,[1 5 4])';
    temp = ([[ 1 1 1]',log(lambda_for_eta)])\log(bbp); %this regression allows for an intercept -
    eta(i) = -temp(2);

    r = corrcoef(bbp,bbp(1)*((lambda_for_eta/lambda_for_eta(1)).^(temp(2))));
    eta_R2(i) = r(1,2)^2;

    if ~any(isnan(q(i,:)))
        [X,cost] = fminsearch(@cost_eta,[1 0.005],opts,q(i,:),wl);
    else
        X = [NaN NaN];
    end
    etaf(i) = X(1);
end

figure
plot(eta,-d,'rx-')
hold on
plot(eta_R2, -d, 'ko-');
legend({'\eta fit linearly, using 442, 510 , 589 nm channels','R^2 of \eta fit'});
ylabel('Depth, m');
hold on
%plot([0.009 0.009],[0 -220], 'k-.');
%plot([0.025 0.025],[0 -220], 'k-.');
title(filename)


%%%%%%%%%%%%End of fitting for eta

figure
hold on
set(gca, 'ColorOrder',c)
plot(wl,q,'o-')
plot(wl,mean(qd),'ko-', 'LineWidth',3)
plot(wl, bbw, 'bo-', 'LineWidth', 3)
zeroline
xlabel('Wavelength, nm')
ylabel(type)
title('Mean spectrum in black, Morel saltwater bb in blue')
%title(filename)

%Plot the profiles
c = jet(6);

figure
subplot(131)
plot(q(:,1),-d,'Color',c(1,:))
hold on
plot(q(:,2),-d,'Color',c(2,:))
zeroline_horizontal
axis([-Inf Inf -220 0])
legend({[type,'442'],[type,'470']},'Location','SouthOutside');
legend boxoff

subplot(132)
plot(q(:,3),-d,'Color',c(3,:))
hold on
plot(q(:,4),-d,'Color',c(4,:))
zeroline_horizontal
axis([-Inf Inf -220 0])
legend({[type,'510'],[type,'589']},'Location','SouthOutside');
legend boxoff
title(filename)

subplot(133)
%plot(q(:,5),-d,'Color',c(5,:))
plot(q(:,5),-d,'k');
hold on
% plot(q(:,6),-d,'Color',c(6,:))
plot(q(:,6),-d,'b');
zeroline_horizontal
axis([-Inf Inf -220 0])
legend({[type,'671']},'Location','SouthOutside');
legend boxoff