% computes factors needed for LSS bead callibration
clear all
%weighting functions
LSSWF=[0 0;10 1;20 4.5;30 12;40 24.5;50 36;60 49.5;70 61;80 72;90 78;...
         100 80;110 78;120 72.5;130 64;140 51;150 39;160 25;170 10.5;180 0];
AC9WF=[0 0;0.69 0;0.71 1; 180 1];
Eco125=[0 0;90 0;105 0;110 0.005;112 0.01;114 0.015;116 0.02;120 0.035;122.5 0.045;...
        125 0.055;126 0.057;127 0.055;130 0.045;133 0.035;135 0.027;137.5 0.02;...
        141 0.01;143 0.005;145 0.002;147 0;160 0;180 0];
BB4WF=[0 0;90 0; 100 0;110 0.02;115 0.04; 120 0.14;125 .28; 130 0.5;135 0.8;140 .98; 143 1; 145 .97;150 0.53;155 0.02;160 0;170 0;180 0];

load phase_func_beads_650

%data from ac-9 at 650

%calculation of 'concentration' factor
ang_new=[0.7:(180-0.7)/601:180]; %ac9 
sin_fac=sin(ang_new*pi/180);
dtheta=(ang_new(2)-ang_new(1))*pi/180;
beta=interp1(ang,phase_func,ang_new,'linear');
ntgrnd=beta.*sin_fac;
nan=length(ang_new);
beta_integral_ac9_650=3/8*(ntgrnd(1)+ntgrnd(nan))+7/6*(ntgrnd(2)+ntgrnd(nan-1))+23/24*(ntgrnd(3)+ntgrnd(nan-2));
beta_integral_ac9_650=(beta_integral_ac9_650+sum(ntgrnd(4:nan-3)))*2*pi*dtheta;


%computing the value of c880

load phase_func_beads_880
beta=interp1(ang,phase_func,ang_new,'linear');
ntgrnd=beta.*sin_fac;
beta_integral_ac9_880=3/8*(ntgrnd(1)+ntgrnd(nan))+7/6*(ntgrnd(2)+ntgrnd(nan-1))+23/24*(ntgrnd(3)+ntgrnd(nan-2));
beta_integral_ac9_880=(beta_integral_ac9_880+sum(ntgrnd(4:nan-3)))*2*pi*dtheta;

%computing the value of b_LSS880
nan=length(ang);
LSSWF2=interp1(LSSWF(:,1),LSSWF(:,2),ang,'cubic');
LSSWF2=LSSWF2/max(LSSWF2);
sin_fac=sin(ang*pi/180);
dtheta=(ang(2)-ang(1))*pi/180;
ntgrnd2=LSSWF2.*phase_func'.*sin_fac;
LSS_integral=3/8*(ntgrnd2(1)+ntgrnd2(nan))+7/6*(ntgrnd2(2)+ntgrnd2(nan-1))+23/24*(ntgrnd2(3)+ntgrnd2(nan-2));
LSS_integral=(LSS_integral+sum(ntgrnd2(4:nan-3)))*2*pi*dtheta;

c650=[0.051; 0.122; 0.298; 0.597; 0.835; 1.121; 1.71; 2.412; 2.689; 2.947]-0.013;
c650b=[0.168; 0.316; 0.989; 1.465; 1.984; 2.507]-0.007;


for i=1:length(c650)
   factor=c650(i)/beta_integral_ac9_650;
   c880(i)=factor*beta_integral_ac9_880;
   LSS(i)=LSS_integral*factor;
end

for i=1:length(c650b)
   factor=c650b(i)/beta_integral_ac9_650;
   c880b(i)=factor*beta_integral_ac9_880;
   LSSb(i)=LSS_integral*factor;
end

%uses a Fournier Forand phase function to convert bLSS to cp(880)
load FF_mu3_5n1_05
nan=length(phi);
LSSWF2=interp1(LSSWF(:,1),LSSWF(:,2),phi,'cubic');
LSSWF2=LSSWF2/max(LSSWF2);
sin_fac=sin(phi*pi/180);
dtheta=(ang(2)-ang(1))*pi/180;
ntgrnd1=beta.*sin_fac;
ntgrnd2=LSSWF2.*beta.*sin_fac;
LSS_integral=3/8*(ntgrnd2(1)+ntgrnd2(nan))+7/6*(ntgrnd2(2)+ntgrnd2(nan-1))+23/24*(ntgrnd2(3)+ntgrnd2(nan-2));
LSS_integral=(LSS_integral+sum(ntgrnd2(4:nan-3)))*2*pi*dtheta;
c_integral=3/8*(ntgrnd1(1)+ntgrnd1(nan))+7/6*(ntgrnd1(2)+ntgrnd1(nan-1))+23/24*(ntgrnd1(3)+ntgrnd1(nan-2));
c_integral=(c_integral+sum(ntgrnd1(4:nan-3)))*2*pi*dtheta;
factorLLStoc=c_integral/LSS_integral