% MIEVSF_NORMGAUSS Mie VSF for Gaussian PSD, normalized to concentration
%   Calculates VSF for Gaussian polydispersion of spheres. Size 
%   distribution is specified by mean diameter and distribution standard
%   deviation. PSD is discretized into a specified number of bins. Mie
%   theory is used to calculate angular scattering matrix and efficiencies
%   for each bin, which are then integrated. The calculated VSF and beam
%   attenuation are normalized by total concentration.
%
%   [THETA,VSF0,C0] = MIEVSF_NORMGAUSS(LAMB,MP,MW,DMEAN,DSTD,DN,NANG,MIEF)
%
%   Outputs
%     THETA is vector of angles (NANG*2-1) from 0 to pi radians;
%     VSF0 and C0 are normalized volume scattering function corresponding
%     to angles in THETA and normalized beam attenuation.
%
%   Inputs
%     LAMB is wavelength of light in air (microns);
%     MP and MW are indices of refraction for the particles and the
%     surrounding medium (MW is real, but MP can be complex);
%     DMEAN and DSTD are the PSD mean diameter and std. deviation (meters);
%     DN is the number of size bins for discretization of PSD;
%     NANG is the number of angles for the Mie code to calculate (THETA
%     will actually be of length NANG*2-1;
%     MIEF is version of Mie code to use, e.g. 'fastmie' or 'bhmie'.
%
% Wayne H. Slade, 10 Oct 2004

% 20 mar 06 In previous version Dn had to be odd, is this true now since
% we're using mmintgrl instead of Simpsons rule?

function [theta, vsf0, c0] = mievsf_normgauss(lambda, m_p, m_w, Dmean, Dstd, Dn, ntheta, miefunc)

	% calculate D bins, total of Dn bins, +/-3 stds
	D = linspace(Dmean-3*Dstd, Dmean+3*Dstd, Dn);
	
	% calculate Gaussian particle size distribution over bins
	ND0 = (Dstd*sqrt(2*pi))^(-1) .* exp(-((D-Dmean).^2)/(2*Dstd^2));
	disp(['integral_ND0 = ' num2str(integral(D, ND0))]);  % hope this is close to unity
	
	x = pi*(D/10^(-6))/lambda*m_w;    % Mie size parameter, lambda is in microns
	m = m_p/m_w;
	
	theta = linspace(0,pi, ntheta*2-1)';
	disp(['dtheta = ' num2str((theta(2)-theta(1))/pi*180)]);  % hope this is close to unity
	
	s1 = zeros(length(theta), Dn);
	s2 = zeros(length(theta), Dn);
	QB = zeros(1, Dn);
	QC = zeros(1, Dn);
	QBACK = zeros(1, Dn);

	% Calculate the (weighted) scattering function (function of theta) for each size bin
	for i = 1:Dn
		if mod(i,fix(Dn/10))==0
			disp(['  Bin ' num2str(i) '/' num2str(Dn)]);
		end
		[s1(:,i), s2(:,i), QB(i), QC(i), QBACK(i)] = feval(miefunc, x(i), m, ntheta);
	end
    
    S11 = repmat(ND0, length(theta),1) .* (0.5*(abs(s1).^2 + abs(s2).^2));
	S11_w = sum(S11, 2);
	phasef = S11_w / (2*pi*integral(theta, S11_w.*sin(theta)));
	disp(['integral_phasef = ' num2str(2*pi*integral(theta, phasef.*sin(theta)))])
	
    b0 = integral(D, ND0.*QB.*(pi*D.^2/4));
    vsf0 = b0 * phasef;

    c0 = integral(D, ND0.*QC.*(pi*D.^2/4));
    
	disp(' ')
return
