% demoMie.m
%
% Demonstrate use of Maetzler's MATLAB Mie calculation.
% Plots scattering, anisotropy, and reduced scattering 
% for a range of wavelengths and range of particle diameters.
% - Steven Jacques 2009
%
% Christian Maetzler, 
% MATLAB Functions for Mie Scattering and Absorption
% University of Bern, Institut fuer Angewandte Physik,
% Research Report No. 2002-08,June, 2002 
% http://arrc.ou.edu/~rockee/NRA_2007_website/Mie-scattering-Matlab.pdf

% Corrected some errors, added table printing function in getMieScatter.m
% -Rob Brown, Rockwell Collins, June 2017
home
disp('=============')
disp('lambda dia rho mus g ')  %header for printing table of interesting parameters for use elsewhere...
%disp('lambda dia rho npar nmed real(m) imag(m) x qext qsca qabs qb asy qratio mus g musp')

%%%%
% USER CHOICES
%%%
nmed = 1; % medium refractive index (air)
npar = 1.33; % particle refractive index (water)
vf   = 0.05; % volume fraction of spheres in medium (PAY ATTENTION TO THIS. Rho (concentration) changes with particle size!)

lambdalambda = .5:.1:14; % [um] (starting wavelength, interval, ending wavelength)
diadia = [5]'; % [um]  list them, [a b c], up to three particle diameters

musgp = zeros(length(lambdalambda),length(diadia),3);

for j=1:length(diadia)
    dia = diadia(j);
    
    for i=1:length(lambdalambda)
        lambda = lambdalambda(i);

        musgp(i,j,:) = getMieScatter(lambda, dia, vf, npar,nmed);

    end % i
    
end % j

%% plot results
figure(100);clf
clr = 'rgb';
sz  = 14;

yname(1).s = '\mu_s [cm^{-1}]';
yname(2).s = 'g [-]';
yname(3).s = '\mu_s'' [cm^{-1}]';
ymax = [3000 1 400];

for k=1:3
    subplot(3,1,k)
    for j=1:length(diadia)
        plot(lambdalambda, musgp(:,j,k),[clr(j) '-'])
        hold on
        y = ymax(k); dy = y/10; x = mean(lambdalambda);
        text(x,y-dy*j,sprintf('dia = %0.3f um', diadia(j)), ...
            'fontsize',sz, 'color',clr(j))
    end
    set(gca,'fontsize',sz)
    xlabel('wavelength \lambda [\mum]')
    ylabel(yname(k).s)
end

disp('done')
