% RINGAVGBETA Average VSF over angle range
%   Calculates average VSF over the range ang1 to ang2 (radians), given the
%   VSF over a set of angles that (hopefully) span ang1 and ang2.
%
%   VSFBAR = RINGAVGBETA(ANG1,ANG2,ANG,VSF)
%
%   ANG1 and ANG2 are length k vectors of ring detector edge angles. VSFBAR
%   is the calculated average betas (normalized or not), averaged over the 
%   (radian) angle ranges ANG1 to ANG2. VSF is the angularly-resolved VSF 
%   over the (radian) angles ANG (VSF and ANG must be same length vectors).
%   VSF is interpolated using pchips over each angle range extending beyond
%   ANG1(k) and ANG2(k), and then integrated numerically. 
%
%   This function is not dummy-proof. It will likely break when given
%   improper parameters.
%
% Wayne H. Slade
% University of Maine, School of Marine Sciences
% Contact wayne.slade@gmail.com

function vsfbar = ringavgbeta(ang1,ang2, vsf_ang,vsf)

% Ensure that ang1 and ang2 are same length vectors
% Ensure that ang1 and ang2 are row vectors
% Ensure that all ang2>ang1 

% Ensure that vsf_ang and vsf have same number of columns
% Ensure that vsf_ang is row vector, while vsf can have many rows


% set(0, 'DefaultLineLineWidth', 0.75)
% set(0, 'DefaultLineMarkerSize', 8)
% set(0, 'DefaultAxesLineWidth', 0.75)
% set(0, 'DefaultTextFontSize', 12)
% set(0, 'DefaultAxesFontSize', 12)
% set(0, 'DefaultAxesBox', 'on')
% hf = figure;
% loglog(vsf_ang, vsf, 'k.'), hold on
% ylim = get(gca,'ylim');

nring = size(ang1,2);
vsfbar = zeros(size(vsf,1), nring);

for kr = 1:nring
	
	% Select subset of vsf_ang spanning and extending beyond [ang1,ang2]
	% over which we'll generate a higher resolution angle vector (ang_hat)
	% to be used for interpolation --
	% vsf_ang point lower than and nearest to ang1
	nearest1 = find(vsf_ang<=ang1(kr), 1,'last');
	% vsf_ang point greater than and nearest to ang2
	nearest2 = find(vsf_ang>=ang2(kr), 1,'first');
	
	% idx1 and idx2 are the end indices into vsf_ang defining the extended
	% angle range for integration
	np = nearest2-nearest1;
	% idx1 = nearest1 - 2*np
	% idx2 = nearest2 + 2*np
	idx1 = nearest1 - 5;
	idx2 = nearest2 + 5;
	if idx1<1
		idx1 = 1;
	end
	if idx2>length(vsf_ang) % (should never happen, but just in case)
		idx2 = length(vsf_ang);
	end
	
	% Do the interpolation and integration for each row of vsf
	ang_hat = linspace(ang1(kr),ang2(kr), np*10);
	
	for kv = 1:size(vsf,1)
		vsf_hat = interp1(vsf_ang(idx1:idx2),vsf(kv,idx1:idx2), ang_hat, 'pchip');
		vsfbar(kv,kr) = integral(ang_hat, vsf_hat.*sin(ang_hat)) ...
					   ./ (cos(ang1(kr)) - cos(ang2(kr)));

		%plot([ang1(kr) ang1(kr) nan ang2(kr) ang2(kr)], [ylim nan ylim])
		%plot(ang_hat, vsf_hat, 'b-')
		%plot([ang1(kr) ang2(kr)], [vsfbar(kv,kr) vsfbar(kv,kr)], 'r')
	end
	
end	
