%Mike Twardowski
%401-783-1787
%mtwardo@wetlabs2.com

%calculates the phase function for a gaussian polydispersion of polystyrene beads
%  with a mean size and st. dev. with a user-defined number of component vsfs

wl=input('Enter wavelength (um, in air): ');

%refractive index of beads relative to air
np=1.5663 + 0.00785/(wl^2) + 0.000334/(wl^4);

%refractive index of medium
nm=1.3481216.*exp(-wl.*1000.*0.0000193465)+2.627716.*exp(-wl.*1000.*0.0153948);

disp(' ');
disp('   ***This program assumes a polydispersion of normally distributed particles***');
d=input('Enter modal particle diameter (um): ');
disp(' ');
disp('   ***Below, must enter a non-zero standard deviation***');
disp('   ***[D - 3*SD] must be positive***');
sd=input('Enter standard deviation of the particle distribution (um): ');
disp(' ');
disp('   ***Below, at least 101 component vsfs are recommended in the calculation (spans');
disp('   ***    plus or minus 3 standard deviations)***');
disp('   ***For a single particle size, enter 1 for the number of component vsfs***');
comp=input('Enter the number of component vsfs to use in the calculation (must be odd): ');

n=np/nm;
angle=[0:1:180]';
S11_w=zeros(181,1);
Qbb_w=0;Qb_w=0;Qc_w=0;

disp(' ');
disp('Number of components calculated:');

for i=1:comp
	%calculate component size (compd) and 
	%gaussian weighting function (weight) over 3 st. dev.
	if comp==1
		compd=d;
		weight=1;
	else
		compd=d-3*sd+(sd/((comp-1)/(2*3))*(i-1));
		weight=((sd*(2*pi)^0.5)^(-1))*exp(-((compd-d)^2)/(2*sd^2))*(sd/((comp-1)/(2*3)));
	end
	
	%calculate Mie size parameter
	rho=pi*compd/wl*nm;

	%calculate scattering vector of component
	[S1 S2 Qb Qc Qback]=bhmie(rho,n,91);
	S11=0.5*((abs(S1)).^2 + (abs(S2)).^2);

    %apply Gaussian weighting to scattering vector
    S11_w=S11_w+S11.*weight;
   
    %calculating the weighted efficiencies
    S11_rad=S11.*sin(angle/180*pi)*2*pi;
    S11_int=sum(S11_rad)*pi/180;
    S11_rad=S11.*sin(angle/180*pi)*2*pi*dang;
    S11_int=3/8*(S11_rad(1)+S11_rad(comp))+7/6*(S11_rad(2)+S11_rad(comp-1))+23/24*(S11_rad(3)+S11_rad(comp-2));
    S11_int=(S11_int+sum(S11_rad(4:comp-3)));
    back_integral=sum(S11_rad(91:181))*pi/180;
    Qbb=back_integral/S11_int*Qb;
    Qbb_w=Qbb_w+Qbb.*weight;
    Qb_w=Qb_w+Qb.*weight;
    Qc_w=Qc_w+Qc.*weight;
    
	disp(i);
end

%S11 is the scattering vector, but it must be normalized so the integrated total is 1
%to obtain the phase function

%calculate sum and multiply by the radian equivalent of 1 degree
    dang=pi/2/(nang-1);
   	nang=180+1;
    S11_rad=S11_w.*sin(angle/180*pi)*2*pi*dang;
    S11_int=3/8*(S11_rad(1)+S11_rad(comp))+7/6*(S11_rad(2)+S11_rad(comp-1))+23/24*(S11_rad(3)+S11_rad(comp-2));
    S11_int=(S11_int+sum(S11_rad(4:comp-3)));
%S11_w_norm is the phase function

phase_func=S11_w_norm;
Qb=Qb_w;Qc=Qc_w;Qbb=Qbb_w;

figure; semilogy(angle,phase_func);axis([0 180 0.0001 1000]);

clear Qbb_w Qc_w Qb_w Qback S11_w_norm S1 S11_int S11 S11_w S11_w_int S11_w_norm S11_w_rad S11_rad S2 back_integral comp  compd  d   i    n   nm   np rho sd weight  wl    




