
% Plot calibrated VSF


%--------------------------------------------------------------------------
clear all, close all

[filename, pathname] = uigetfile('lisstcalrings_*.mat', 'LISST Calibration Data')
load([pathname filename])

% ID = '2f'
% DATAFILE = 'beads_dec04\beads_2u_8.log'
% ZSCATFILE = 'beads_dec04\zsc_diw_122104_1.asc'
% DCALFILE = 'beads_dec04\Ringarea_1102.asc'
% BADROWS = [];
% BADRINGS = [1:6];
% MIEFILES = 'miebeads/miemonte_4202A(26544)'
% REJECTSNR = 1;
% ylim = [4 2e2]

% ID = '20h'
% DATAFILE = 'beads_dec04\beads_20u_8.log'
% ZSCATFILE = 'beads_dec04\zsc_diw_122004_2.LOG'
% DCALFILE = 'beads_dec04\Ringarea_1102.asc'
% BADROWS = [];
% BADRINGS = [5 6];s_20050213T195941/miemonte_4220A(26484)'
% REJECTSNR = 1;
% ylim = [1 5e3]

ID = '50a'
DATAFILE = 'beads_dec04\beads_50u_2.log'
ZSCATFILE = 'beads_dec04\zsc_diw_122104_3.asc'
DCALFILE = 'beads_dec04\Ringarea_1102.asc'
BADROWS = [];
BADRINGS = [];
MIEFILES = 'miebeads/miemonte_6050(25307)'
REJECTSNR = 1;
ylim = [1 5e4]



%--------------------------------------------------------------------------
[cscat_angm,cscat,cmeas,snr, cscat_ang1,cscat_ang2,cscat_rings] = readpscat( ...
	DATAFILE, ZSCATFILE, DCALFILE, BADROWS, [], 1);

cscat(:,BADRINGS) = nan;	
		
% read Monte Carlo set for this run
vsf0 = load([MIEFILES '_vsf0.txt']);
vsf0_ang = load([MIEFILES '_psi.txt']);  %in radians
cmie0 = load([MIEFILES '_beamc0.txt']);

% 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);

cmeas_m = median(cmeas);
cmeas_s = (prctile(cmeas,84)-prctile(cmeas,16))/2;
cscat_m = median(cscat);
cscat_p84 = prctile(cscat,84);
cscat_p16 = prctile(cscat,16);

vsfbar0_m = median(vsfbar0);
vsfbar0_s = (prctile(vsfbar0,84)-prctile(vsfbar0,16))/2;
vsf0_m = median(vsf0);
vsf0_s = (prctile(vsf0,84)-prctile(vsf0,16))/2;
cmie0_m = median(cmie0);
cmie0_s = (prctile(cmie0,84)-prctile(cmie0,16))/2;
    
vsfbar_m = cmeas_m * vsfbar0_m./cmie0_m;
vsfbar_s = (abs(cmeas_s)/cmeas_m + abs(vsfbar0_s)/vsfbar0_m + abs(cmie0_s)/cmie0_m)*vsfbar_m;
    
vsf_m = cmeas_m * vsf0_m./cmie0_m;
vsf_s = (abs(cmeas_s)/cmeas_m + abs(vsf0_s)/vsf0_m + abs(cmie0_s)/cmie0_m)*vsf_m;	



set(0, 'DefaultLineLineWidth', 0.75)
set(0, 'DefaultLineMarkerSize', 12)
set(0, 'DefaultAxesLineWidth', 0.75)
set(0, 'DefaultTextFontSize', 9)
set(0, 'DefaultAxesFontSize', 9)
set(0, 'DefaultTextHorizontalAlignment', 'left')
set(0, 'DefaultAxesBox', 'on')

hf1 = figure, hold on
set(gcf, 'paperunits', 'inches', 'paperorientation', 'portrait', 'papersize', [3.1 2.9])
set(gcf, 'paperpositionmode', 'manual', 'paperposition', [0.05 0.05 3.0 2.8])
plot(vsf0_ang./pi.*180, vsf_m, 'k:')
set(gca, 'xlim', [0 20], 'xscale','log', 'yscale','log')
%set(gca, 'ylim', [3 200])
%ylim = get(gca,'ylim')

hf2 = figure, hold on
set(gcf, 'paperunits', 'inches', 'paperorientation', 'portrait', 'papersize', [3.1 2.9])
set(gcf, 'paperpositionmode', 'manual', 'paperposition', [0.05 0.05 3.0 2.8])


for kr = 1:length(cscat_rings)

	figure(hf1) % -----------------------------------------------------------------------
	
	plot(repmat(cscat_ang1(kr),1,2)/pi*180, [ylim(1) ylim(2)]', 'k:') 
	delx = (cscat_ang2(kr)-cscat_ang1(kr));
	vxmie = [cscat_ang1(kr) cscat_ang2(kr) ...
			cscat_ang2(kr) cscat_ang1(kr)]/pi*180;
	vymie = [repmat(vsfbar_m(kr)+3.*vsfbar_s(kr),1,2) repmat(vsfbar_m(kr)-3.*vsfbar_s(kr),1,2)];
	vxlisst = [cscat_ang1(kr)+0.1*delx cscat_ang1(kr)+0.9*delx ...
			cscat_ang1(kr)+0.9*delx cscat_ang1(kr)+0.1*delx]/pi*180;
	vylisst = [repmat(cscat_p84(kr),1,2).*(chi1_rf2(kr)+3*chi1_se_rf2(kr)) ...
				repmat(cscat_p16(kr),1,2).*(chi1_rf2(kr)-3*chi1_se_rf2(kr))];
	patch(vxmie,vymie,[1 1 1])
	patch(vxlisst,vylisst,[1 0 0])
	
	
	figure(hf2)
		plot(chi1_rf2(kr).*cscat_m(kr), vsfbar_m(kr),'o', ...
			'markersize',6, 'MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0 0 1])
		plot([cscat_p16(kr).*(chi1_rf2(kr)-3*chi1_se_rf2(kr)) cscat_p84(kr).*(chi1_rf2(kr)+3*chi1_se_rf2(kr))], ...
			 repmat(vsfbar_m(kr),1,2),'k-')
		plot(repmat(chi1_rf2(kr).*cscat_m(kr),1,2), ...
			[vsfbar_m(kr)-3*vsfbar_s(kr) vsfbar_m(kr)+3*vsfbar_s(kr)], 'k-')

		
end

figure(hf1)
plot(repmat(cscat_ang2(kr),1,2)/pi*180, [ylim(1) ylim(2)]', 'k:') 
plot(vsf0_ang./pi.*180, vsf_m, 'k-', 'linewidth', 2*get(0, 'DefaultLineLineWidth'))
set(gca, 'xlim',[0.05 20])
set(gca, 'xtick', [0.1 1.0 10.0])
set(gca, 'xticklabel', {'0.1' '1.0' '10.0'})
set(gca, 'ylim', ylim)
ylabel('\langle\beta_{LISST}\rangle_i [m^{-1} sr^{-1}]')
xlabel('\psi [degrees]')


figure(hf2)
set(gca, 'xscale','log', 'yscale','log')
set(gca, 'xlim',ylim, 'ylim',ylim)
plot(ylim, ylim, 'r-')
xlabel('\langle\beta_{LISST}\rangle_i [m^{-1} sr^{-1}]')
ylabel('\langle\beta_{Mie}\rangle_i [m^{-1} sr^{-1}]')

% print(hf1, '-depsc', ['fig-calibratedbeta_vs_psi_' ID])
% print(hf2, '-depsc', ['fig-calibratedbeta_scatter_' ID])

