% RUN_LISSTCALRINGS  Calibration for LISST to VSF
%   This program reads the results of LISST experiments with known PSD from
%   standard beads and compares them with the expected Mie solutions 
%   for the same PSD and index of refraction. Calibration factors are then 
%   determined that can be used to later convert other LISST scat data into
%   the absolute VSF at the LISST-measured angles.
%
%   See Slade and Boss, "Calibrated near-forward volume scattering function 
%   obtained from the LISST particle sizer," Optics Express (2006).
%
%   http://misclab.umeoce.maine.edu/code/lisstvsf.html
%
% Wayne H. Slade
% University of Maine, School of Marine Sciences
% Contact wayne.slade@gmail.com

[runlist_num,runlist_txt] = xlsread('runlist.xls', 'calibration');
runlist = runlist_txt;
runlist(:,1) = num2cell(runlist_num(:,1));
clear runlist_num runlist_txt
runlist(1:runlist{1,1},:) = [];
for jj = 1:size(runlist,1)
    runlist{jj,2} = eval(runlist{jj,2});
    runlist{jj,8} = eval(runlist{jj,8});
    runlist{jj,9} = eval(runlist{jj,9});
end

CSCAT = cell(size(runlist,1),1);
CSCAT_ANG1 = cell(size(runlist,1),1);
CSCAT_ANG2 = cell(size(runlist,1),1);
CSCAT_RINGS = cell(size(runlist,1),1);
CMEAS = cell(size(runlist,1),1);
VSF0 = cell(size(runlist,1),1);
VSFBAR0 = cell(size(runlist,1),1);
VSF0_ANG = cell(size(runlist,1),1);
CMIE0 = cell(size(runlist,1),1);

% (see run_snrtest.m) 
rejectsnr = 10;

tt = char(9);
for kf = 1:size(runlist,1)

	disp('----------------------------------------------------------------------')
	
	[cscat_angm,cscat,cmeas,snr, cscat_ang1,cscat_ang2,cscat_rings] = readpscat( ...
		runlist{kf,5}, runlist{kf,6}, runlist{kf,7}, ...
		    runlist{kf,8}, runlist{kf,9}, rejectsnr);
	
	% Discard any remaining scans with negative in cscat or cmeas
	%     badscan = any(cscat'<=0)' | cmeas<=0;
	% 	if any(badscan==logical(1))
	% 		cscat(badscan,:) = [];
	% 		cmeas(badscan) = [];
	% 		disp(['** Discarded due to negative in scan:  ' num2str(find(badscan)')])
	% 	end
	
	
	% read Monte Carlo set for this run
    if kf>1 & strcmp(runlist{kf,10}, runlist{kf-1,10})
        vsf0 = VSF0{kf-1};
        vsf0_ang = VSF0_ANG{kf-1};
        cmie0 = CMIE0{kf-1};
    else
        vsf0 = load([runlist{kf,10} '_vsf0.txt']);
        vsf0_ang = load([runlist{kf,10} '_psi.txt']);  %in radians
        cmie0 = load([runlist{kf,10} '_beamc0.txt']);
        %vsf0_ang = 1.33*vsf0_ang;
	end


	% Calculate the expected average beta for each LISST detector 
	% (integrate over the angle range of the ring)
	vsfbar0 = ringavgbeta(cscat_ang1,cscat_ang2, vsf0_ang,vsf0);
	
	% Could include correction coefficient here to account for finite
	% acceptance angle of transmission sensor? (less than a percent or so)
	%cmie0 = cmie0.*delta;
	
	% Hold on to pertinent data for each run
	CSCAT{kf} = cscat;
    CSCAT_ANG1{kf} = cscat_ang1;
    CSCAT_ANG2{kf} = cscat_ang2;
    CSCAT_RINGS{kf} = cscat_rings;
    CMEAS{kf} = cmeas;
    VSF0{kf} = vsf0;
    VSF0_ANG{kf} = vsf0_ang;
    VSFBAR0{kf} = vsfbar0;
    CMIE0{kf} = cmie0;
	
	
end             

CMEAS_m = cell(size(runlist,1),1);
CMEAS_s = cell(size(runlist,1),1);
CSCAT_m = cell(size(runlist,1),1);
CSCAT_p84 = cell(size(runlist,1),1);
CSCAT_p16 = cell(size(runlist,1),1);
VSFBAR_m = cell(size(runlist,1),1);
VSFBAR_s = cell(size(runlist,1),1);

% Loop through the runs to calculate error bars and regression coefs
for kf = 1:size(runlist,1)

    CMEAS_m{kf} = median(CMEAS{kf});
    CMEAS_s{kf} = (prctile(CMEAS{kf},84)-prctile(CMEAS{kf},16))/2;
    CSCAT_m{kf} = median(CSCAT{kf});
    CSCAT_p84{kf} = prctile(CSCAT{kf},84);
    CSCAT_p16{kf} = prctile(CSCAT{kf},16);

    vsfbar0 = VSFBAR0{kf};
    vsfbar0_m = median(VSFBAR0{kf});
    vsfbar0_s = (prctile(VSFBAR0{kf},84)-prctile(VSFBAR0{kf},16))/2;

	vsf0 = VSF0{kf};
    vsf0_m = median(VSF0{kf});
    vsf0_s = (prctile(VSF0{kf},84)-prctile(VSF0{kf},16))/2;

    cmie0_m = median(CMIE0{kf});
    cmie0_s = (prctile(CMIE0{kf},84)-prctile(CMIE0{kf},16))/2;
    
    VSFBAR_m{kf} = CMEAS_m{kf} * vsfbar0_m./cmie0_m;
    VSFBAR_s{kf} = (abs(CMEAS_s{kf})/CMEAS_m{kf} + abs(vsfbar0_s)/vsfbar0_m + abs(cmie0_s)/cmie0_m)*VSFBAR_m{kf};
    
	VSF_m{kf} = CMEAS_m{kf} * vsf0_m./cmie0_m;
    VSF_s{kf} = (abs(CMEAS_s{kf})/CMEAS_m{kf} + abs(vsf0_s)/vsf0_m + abs(cmie0_s)/cmie0_m)*VSF_m{kf};	

end



tt = char(9);
disp(['Size' tt 'Run' tt 'Stock' tt  '$n_{LISST}$' tt '$c_{LISST}$' tt '$\delta c_{LISST}/c_{LISST}$' tt 'Derived $N_0$' ' \\'])
for kf = 1:size(runlist,1)
    
    cmeas_n = num2str(length(CMEAS{kf}));
    cmeas_m = num2str(CMEAS_m{kf}, '%0.2f');
    cmeas_s = num2str(CMEAS_s{kf}./CMEAS_m{kf}, '%0.2f');
    n0 = num2str(median(CMEAS{kf})/median(CMIE0{kf}), '%0.2e');
	disp([tt runlist{kf,4} tt tt (cmeas_n) tt (cmeas_m) tt (cmeas_s) ...
			tt (n0) ' \\'])
end


% Do regression for each ring
chi1_rf = nan(1,32);
chi1_se_rf = nan(1,32);
chi1_rmse_rf = nan(1,32);
chi1_mad_rf = nan(1,32);

for kr = [1:32]

	disp(['DETECTOR RING ' num2str(kr)])

	x = [];
	dx = [];
	y = [];
	dy = [];
	
	for kf = 1:size(runlist,1)
		x = [x; CSCAT_m{kf}(CSCAT_RINGS{kf}==kr)];
		dx = [dx; 0.5.*(CSCAT_p84{kf}(CSCAT_RINGS{kf}==kr)-CSCAT_p16{kf}(CSCAT_RINGS{kf}==kr))];
		y = [y; VSFBAR_m{kf}(CSCAT_RINGS{kf}==kr)];
		dy = [dy; VSFBAR_s{kf}(CSCAT_RINGS{kf}==kr)];
	end
	
	[b,stats] = robustfit(x,y, 'welsch',5, 'off');
	%[b,stats] = robustfit(x,y, 'bisquare',4.685, 'off')
	chi1_rf(kr) = b(1);
	chi1_se_rf(kr) = stats.se;
	chi1_rmse_rf(kr) = sqrt(sum((y-x.*chi1_rf(kr)).^2)/length(y));
	chi1_mad_rf(kr) = mean(abs((y-x.*chi1_rf(kr))./y));
	
	xhat = logspace(log10(min(x)),log10(max(x)), 100);
	yhat_rf = chi1_rf(kr).*xhat;

end
	
save(['lisstcalrings_' datestr(now,30)])



	

