% function [bb_lambda_data]=get_bb_from_ecovsf(lambda,apg_lambda_data,beta_data)
%
% calculate bb from the ecovsf data in accordance with Twardowski et al
% add bbw to data
%
% INPUTS: 
%		lambda: sensor wavelength (ex: 530)
%     apg_lambda: apg data at sensor wavelength (n x 1 array)
%     beta_data: beta values at 100 125 and 150 deg (n x 3 array)
%      
% OUTPUTS:
%		bb_lambda_data: bb values at sensor wavelength (n x 1 array)    
%     
% Written by Scott Pegau, Jan 2001, OSU
% Modified by Francois Baratange May 2001, OSU

function [bb_lambda_data]=get_bb_from_ecovsf(lambda,apg_lambda_data,beta_data)

% transpose (might not be needed) to (n_sample,n_variable)

[n,r]=size(beta_data);

bb100=(beta_data(:,1)-128.6)*3.091e-5;
bb125=(beta_data(:,2)-127.716)*2.78e-5;
bb150=(beta_data(:,3)-132.9)*1.159e-5;

% correct measurements for attenuation
bb100=bb100.*exp(.0314*apg_lambda_data);
bb125=bb125.*exp(.0441*apg_lambda_data);
bb150=bb150.*exp(.0804*apg_lambda_data);

% convert bb data to flux data
flux(:,1)=2.*pi.*sin(pi*100/180).*bb100;
flux(:,2)=2.*pi.*sin(pi*125/180).*bb125;
flux(:,3)=2.*pi.*sin(pi*150/180).*bb150;
flux(:,4)=0;

%fit flux data with a 3rd order polynomial and integrate from 90-180 deg
rad=[pi*100/180 pi*125/180 pi*150/180 pi];
for i=1:n
   [dwrpolycoeff(i,1:4) s]=polyfit(rad,flux(i,:),3);
   coeff=dwrpolycoeff(i,1:4);
   bb_lambda_data(i)=coeff(1)*(pi^4)/4+coeff(2)*(pi^3)/3+...
          coeff(3)*(pi^2)/2+coeff(4)*pi-(coeff(1)*((pi/2)^4)/4+...
          coeff(2)*((pi/2)^3)/3+coeff(3)*((pi/2)^2)/2+coeff(4)*pi/2);
end

bb_lambda_data=bb_lambda_data+1.38*((lambda/500)^(-4.32))*2*pi*(1+(1-0.09)/(3*(1+0.09)))*10^(-4);
bb_lambda_data=bb_lambda_data';
return