function [factor]=compute_factor(center_angle,nominal_wl);
% 2005.12
% Emmanuel Boss
% emmanuel.boss@maine.edu
% School of Marine Sciences, University of Maine, USA
%
% Calculates the calibration factor for Eco-type sensors. This factor
% represents the angle and wavelength weighted phase function. The
% wavelength and bead-size computation has been done by
% 'phase_beads_lambda_D'. The output takes acount of the finite acceptance
% angle of the ac-9, here assumed to be 0.7degrees.
%
%
% INPUT
% center_angle for backscattering                     (center-angle)    [degree]
% nominal wavelength for backscattering               (nominal_wl)    [nm]

% OUTPUT
% weighted phase function                            (factor)   [sr^-1]

infile=[num2str(nominal_wl),'nm_LED_1_992.mat'];
load(infile)
%load('EcoVSF650nm_LED.mat')
m=length(PF);
nang=100*3+1;

%insuring the phase function obeys normalization
sin_fac=sin(ang*pi/180)';
dtheta=(ang(2)-ang(1))*pi/180;
ntgrnd=PF.*sin_fac;
S11_integral=3/8*(ntgrnd(1)+ntgrnd(m))+7/6*(ntgrnd(2)+ntgrnd(m-1))+23/24*(ntgrnd(3)+ntgrnd(m-2));
S11_integral=(S11_integral+sum(ntgrnd(4:m-3)))*2*pi*dtheta;
PF=PF/S11_integral; %phase function

%Phase function for ac-9 acceptance angle
angb=[0.8:0.1:180];
N=length(angb);
PFb=interp1(ang,PF,angb,'cubic')';
sin_facb=sin(angb*pi/180)';
dtheta=0.1;
ntgrnd=PFb.*sin_facb;
S11_integral=3/8*(ntgrnd(1)+ntgrnd(N))+7/6*(ntgrnd(2)+ntgrnd(N-1))+23/24*(ntgrnd(3)+ntgrnd(N-2));
S11_integral=(S11_integral+sum(ntgrnd(4:N-3)))*2*pi*dtheta;
PFb=PFb/S11_integral; 

%calculate angle weighted phase function for full phase function
S11_rad=PF.*W2(ang',center_angle)*0.3*pi/180;
S11_int=3/8*(S11_rad(1)+S11_rad(m))+7/6*(S11_rad(2)+S11_rad(m-1))+23/24*(S11_rad(3)+S11_rad(m-2));
S11_int=(S11_int+sum(S11_rad(4:m-3)));

%calculate angle weighted phase function for ac-9 based phase function
S11_rad=PFb.*W2(angb',center_angle)*0.1*pi/180;
S11_intb=3/8*(S11_rad(1)+S11_rad(m))+7/6*(S11_rad(2)+S11_rad(N-1))+23/24*(S11_rad(3)+S11_rad(N-2));
S11_intb=(S11_int+sum(S11_rad(4:N-3)));

%calculate normalization integral to divide the weighted phase functions by
S_rad=W2(ang',center_angle)*0.3*pi/180;
S_int=3/8*(S_rad(1)+S_rad(m))+7/6*(S_rad(2)+S_rad(m-1))+23/24*(S_rad(3)+S_rad(m-2));
S_int=(S_int+sum(S_rad(4:m-3)));


factor_b=S11_int/S_int
factor=S11_intb/S_int

