function [PF]=PF_beads(mean,sd)

%calculates the wavelengths and size averaged phase function for a gaussian polydispersion of polystyrene beads
%with a mean size and st. dev. with a user-defined number of component
%vsfs, a given mean and stdev for the particles and a wavelength weighting
%function (W1).

%Emmanuel Boss
%207-581-4378
%emmanuel.boss@maine.edu

%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

lambda=400:1:700;
plot(lambda,W1(lambda))
A=W1(lambda);
I1=find(A>max(A));
wl_max=lambda(I1);
I=find(A>max(A)/2);
min_wl=wl_max-3*(wl_max-min(lamda(I)));
max_wl=wl_max+3*(wl_max-min(lamda(I)));
wl=[min_wl:(max_wl-min_wl)/199:max_wl];
%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);
n=np/nm;
min_D=mean-3*sd;
max_D=mean+3*sd;
max_x=pi*max_D/min_wl/min(n);
min_x=pi*min_D/max_wl/max(n);

nang=100*3+1;
dang=pi/2/(nang-1);
nan=2*nang-1;
for i=1:300
	%calculate Mie size parameter
	x=min_x+(i-1)*(max_x-min_x)/299
	[S1 S2 Qb Qc Qback]=bhmie(x,n,nang);
	S11(i,:)=0.5*((abs(S1)).^2 + (abs(S2)).^2);
    sin_fac=sin(ang*pi/180);
    dtheta=(ang(2)-ang(1))*pi/180;
    ntgrnd=S11(i,:).*sin_fac;
    S11_integral=3/8*(ntgrnd(1)+ntgrnd(nan))+7/6*(ntgrnd(2)+ntgrnd(nan-1))+23/24*(ntgrnd(3)+ntgrnd(nan-2));
    S11_integral=(S11_integral+sum(ntgrnd(4:nan-3)))*2*pi*dtheta;
    S11(i,:)=S11(k,:)/S11_integral*Qb; %VSF
    Qb_(i)=Qb;
    xx(i)=x;
	disp(i);
end

beta_=zeros(601,1);
Qb_=0;
for i=1:200
    for j=1:200
        D=min_D+(i-1)*(max(D)-min(D))/199;   
        weight=((sd*(2*pi)^0.5)^(-1))*exp(-((D-min)^2)/(2*sd^2))*W1(wl(j);
        x=pi*D/wl(j)/n(j);
        beta_=beta_+weight*interp1(xx,S11,x,'linear);
        Qb_=Qb_+interp1(xx,Qb,x,'linear);
    end
end


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





