%This program 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.

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


clear all
close all

%properties of particles, mean size and standard deviation:
mean=1.992;
sd=0.022;
%mean=2;
%sd=0.08;
disp(['the mean bead size is: ',num2str(mean),'um'])
disp(['the standard deviation around the mean size is: ',num2str(sd), 'um'])

%nominal wavelength of sensor 
nominal_wl=715;
disp(['the nominal wavelength of sensor is: ',num2str(nominal_wl), 'nm'])

lambda=[380:.5:800]/1000;
figure(1);
plot(lambda*1000,W1(lambda*1000,nominal_wl),'r')
A=W1(lambda*1000,nominal_wl);
title('wavelength weighting function');
xlabel('wavelength');
ylabel('arbitrary amplitude');

%Range of wavelengths of which computation are done 3FMHW around wavelength
%of maximum response.
I1=find(A==max(A));
wl_max=lambda(I1);
I=find(A>=max(A)/2);
min_wl=wl_max-3*(wl_max-min(lambda(I)));
max_wl=wl_max+3*(wl_max-min(lambda(I)));
wl=[min_wl:(max_wl-min_wl)/99:max_wl];

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

%refractive index of medium (water)
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*max(nm);
min_x=pi*min_D/max_wl*min(nm);

nang=100*3+1;
dang=pi/2/(nang-1);
nan=2*nang-1;
ang=[0:dang:pi]*180/pi;

beta_=zeros(601,1);
Qb2=0;
for i=1:100
    for j=1:100
        D=min_D+(i-1)*(max_D-min_D)/99;   
        weight=((sd*(2*pi)^0.5)^(-1))*exp(-((D-mean)^2)/(2*sd^2))*W1(wl(j)*1000,nominal_wl);
        x=pi*D*nm(j)/wl(j);
        [S1 S2 Qb Qc Qback]=bhmie(x,n(j),nang);
	    S11=0.5*((abs(S1)).^2 + (abs(S2)).^2);
        W(i,j)=weight;
        beta_=beta_+weight*S11;
    end
    disp(i)
end
        sin_fac=sin(ang*pi/180)';
        dtheta=(ang(2)-ang(1))*pi/180;
        ntgrnd=beta_.*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;
        PF=beta_/S11_integral; %phase function

figure(2)
semilogy(ang,S11,'r');
ylabel('phase function [sr^-1]');
xlabel('angle');
DD=[min_D:(max_D-min_D)/99:max_D];
%outfile=['EcoVSF',num2str(nominal_wl),'nm_LED.mat'];
outfile=[num2str(nominal_wl),'nm_LED_1_992.mat'];
save(outfile, 'ang', 'PF', 'W', 'mean', 'sd', 'wl', 'DD')


