% RUN_MIEMONTE Calculate Mie VSF and attenuation for bead suspensions
%   This program determines theoretical (concentration normalized) VSF and
%   beam attenuation for Gaussian polydispersions of beads, using Mie
%   theory. The VSFs incorporate uncertainty in bead diameter using Monte
%   Carlo analysis.
%
%   See Slade and Boss, "Calibrated near-forward volume scattering function 
%   obtained from the LISST particle sizer," Optics Express (2006).
%
%   http://misclab.umeoce.maine.edu/code/lisstvsf.html
%
% Wayne H. Slade
% University of Maine, School of Marine Sciences
% Contact wayne.slade@gmail.com



% saved = ['miebead_' datestr(now,30) '/']
saved = ['miebeads/']
mkdir(saved)
diary([saved 'diary.txt'])


name =   {'4K-02(25322)', '4202A(26544)', '4220A(26484)', '6050(25307)', '4K100(25342)'}
d_mean = [1.992           2.001           20.00           50.4           99.2] % bead suspension mean diameter
d_unc =  [0.025           0.025           0.10            1.0            1.4] % uncertainty in mean diameter (+/-)
d_std =  [0.022           0.022           0.20            1.6            1.7] % standard deviation of size distribution


lambda = 0.670  % micron
i = sqrt(-1)
m_p = 1.5663 + 0.00785./(lambda.^2) + 0.000334./(lambda.^4)
m_w = 1.3481216.*exp(-lambda.*1000.*0.0000193465) + ...
                                 2.627716.*exp(-lambda.*1000.*0.0153948)
							 
comp = 100 % number of PSD bins to use in Mie calcs
nang = 2001 % number of VSF angles to use in Mie calcs
nmc = 200  % Monte Carlo sample size (number of random PSD)

VSF0 = zeros(nmc, nang*2-1);
BEAMC0 = zeros(nmc, 1);

% Calculate set of Mie theory VSFs for each bead size
for k = 1:length(name)
    
	disp(name{k})
    d = random('uniform', d_mean(k)-d_unc(k), d_mean(k)+d_unc(k), 1, nmc);
    for kd = 1:length(d)
        disp(['Monte Carlo realization ' num2str(kd) '/' num2str(length(d))])
        [theta, vsf0, c0] = mievsf_normgauss(lambda, m_p, m_w, d(kd)*1e-6, d_std(k)*1e-6, comp, nang, 'fastmie');
        VSF0(kd,:) = vsf0';
        BEAMC0(kd) = c0;
        disp('--------------------------------------------------')
    end

    theta = theta';
    save([saved 'miemonte_' name{k} '_vsf0.txt'], 'VSF0', '-ascii')
    save([saved 'miemonte_' name{k} '_psi.txt'], 'theta', '-ascii')
    save([saved 'miemonte_' name{k} '_beamc0.txt'], 'BEAMC0', '-ascii')
    
    figure, loglog(theta, VSF0)
    xlabel('\psi [rad]'), ylabel('\beta(\psi) / N_0 [sr^{-1}]')
    saveas(gcf, [saved 'miemonte_' name{k}], 'fig')
    saveas(gcf, [saved 'miemonte_' name{k}], 'jpg')
    
end

diary off