%%      Mie Theory For A Single Sphere: Calculates Efficiencies, 
%%    Asymmetry Factors, And Scattering Phase Function & Plots Efficiencies

% Rationale:    Computation of Mie Efficiencies for given complex 
%   refractive-index ratio, (m = n + ki"), and size parameter, (x = k0*a), 
%   where k0 is the wave number in ambient medium (k0 = 2pi/wavelength), 
%   a = sphere radius, using complex Mie Coefficients a_n and b_n for 
%   n = 1 to nmax.  Computation of angle functions pi_n and tau_n for
%   -1 <= cos(theta) <= 1, n1 integer from 1 to nmax.
%               Computation of a spectrum of Mie Efficiencies with the 
%   intent to plot them over a range of size parameter values.

% Sources:  Efficiencies: S. Bohren and Huffman (1983) BEWI:TDD122, p.
%               103,119-122,477.
%           Mie coefficients: Eq. (4.88) of Bohren and Huffman (1983), 
%               BEWI:TDD122 using the recurrence relation (4.89) for Dn on 
%               p. 127 and starting conditions as described in Appendix A.
%               C. Mtzler, July 2002
%           Angle Functions: Bohren and Huffman (1983), p. 94 - 95

% Results:  n, k, x, efficiencies for extinction (qext), scattering (qsca), 
%   absorption (qabs), backscattering (qb), asymmetry parameter (asy),
%   qratio (qb/qsca), and angle functions pi_n and tau_n.

clc;
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Start User-Defined Inputs: %%%%%%%%%%%%%%%%%%%%%%%%%


r = 1.0;              % Radius of sphere in MICROmeters

lamda = 2.0;         % Wavelength of radiation in MICROmeters

m1 = 5 + 0.4i;       % Complex refractive index of particle, m1 = n1 + ik1, 
                        % n1: refractive index of particle
                        % k1: absorption index of particle
                        
m2 = 1.0;            % Surrounding material's refractive index, m2 = n2            

theta = 0;       % Scattering Angle: specify in DEGREES
                        % If there is no angle to specify, 
                        % input: "theta = 500i"

                        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% The following section is to plot over a range of size parameters only
% If not plotting, input: "xi = 500i"

xi = 500i;      % Initial size parameter

xf = 10;         % Final size parameter

dstep = 0.01;    % Step size, the smaller it is, the higher the resolution,
                        % but it takes longer to compute.

m1_p = 5 + 0.4i;    % Complex refractive index of particle, m1 = n1 + ik1, 
                        % _p denotes plotting 
                        % n: refravtive index 
                        % k: absorption index of particle
                        
m2_p = 1.0;         % Surrounding material's refractive index, m2 = n2   
                        % _p denotes plotting
                        
                        
%%%%%%%%%% End User-Defined Input: DO NOT Modify Below This Line %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


m1 = real(m1)/m2+imag(m1)*1i;        % compute for relative refractive index
m2 = real(m1_p)/m2_p+imag(m1_p)*1i;  % compute for relative refractive index 
                                        % (plotting)
x = 2*pi*r/lamda;       % Size parameter (nondimensional)
nmax = round(x+4*x^(1/3)+2);   % Truncation parameter for infinite sums

if isreal(xi)   % Sets relative refractive index according to plotting
    m = m2;
    display('Complex Refractive Index')
    disp(m);
else
    disp('Size Parameter')
    disp(x)
    m = m1;
    display('Complex Refractive Index')
    disp(m)
end

if isreal(xi) == 0  % Differentiates plotting and non-plotting situations
    xsize = 1;
else
    xsize = xi:dstep:xf;     % Array of size parameters (for plotting)
end

for k = 1: length(xsize)    % Will loop if plotting is occuring
    if isreal(xi) == 0  % Sets size parameter dependent on plot condition
        xv = x(k);          % No plotting
    else
        xv = xi:dstep:xf;   % Plotting
    end
    x = xv(k);              % Sets size paramter for each iteration
    
    %%% Calculates Mie coefficients, asymmetry factors and efficiencies %%%
if x == 0       % Done to avoid singularity for x
    if isreal(xi) == 0 % Differentiates plotting & non-plotting situations
        disp(['qext ' 'qsca ' 'qabs ' 'qb ' 'asy ' 'qratio'])
        disp(['  0   ' ' 0   ' ' 0   ' '0  ' '0   ' '1.5  ']);
    end
elseif x > 0    % Normal running situation
    n1 = nmax-1;
    n =(1:nmax);
    cn = 2*n+1; 
    c1n = n.*(n+2)./(n+1); 
    c2n = cn./n./(n+1);
    x2 = x.*x;
    
    
    %%% Calculates a_n and b_n mie coefficients %%%
    z = m.*x;     
    nmx = round(max(nmax,abs(z))+16); 
    n =(1:nmax); 
    nu = (n+0.5); 
    sx = sqrt(0.5*pi*x);
    px = sx.*besselj(nu,x);         % First kind Bessel function
    p1x = [sin(x), px(1:nmax-1)];   % Spherical Bessel function
    chx = -sx.*bessely(nu,x);       % Second kind Bessel function
    ch1x = [cos(x), chx(1:nmax-1)]; % Spherical Bessel function
    gsx = px-1i*chx; 
    gs1x=p1x-1i*ch1x;
    dnx(nmx) = 0+0i;
    for j = nmx:-1:2    % Computation of Dn(z) from (4.89) of B+H (1983)
        dnx(j-1) = j./z-1/(dnx(j)+j./z);
    end;
    dn = dnx(n);        % Dn(z), n=1 to nmax
    da = dn./m+n./x; 
    db = m.*dn+n./x;
    an =(da.*px-p1x)./(da.*gsx-gs1x);   % Solve for a_n (Equation 4.53)
    bn =(db.*px-p1x)./(db.*gsx-gs1x);   % Solve for b_n (Equation 4.53)
    f = [an; bn];   % Matrix of Mie coefficients, a_n and b_n
    
    
    %%% Calculates Efficiencies and Asymmetry Factors %%%
    anp = (real(f(1,:)));   % Real component of a_n
    bnp = (real(f(2,:)));   % Real component of b_n
    anpp = (imag(f(1,:)));  % Imaginary component of a_n
    bnpp = (imag(f(2,:)));  % Imaginary component of a_n
    g1(1:4,nmax) = [0; 0; 0; 0]; % g1 = matrix for asymmetry parameters
    g1(1,1:n1) = anp(2:nmax);    
    g1(2,1:n1) = anpp(2:nmax);
    g1(3,1:n1) = bnp(2:nmax);
    g1(4,1:n1) = bnpp(2:nmax);   
    
    % Calculate: q extinction (Equation 4.62)
    dn = cn.*(anp+bnp);
    q = sum(dn);
    qext = 2*q/x2;
    
    % Calculate: q scatter (Equation 4.61)
    en = cn.*(anp.*anp+anpp.*anpp+bnp.*bnp+bnpp.*bnpp);
    q = sum(en);
    qsca = 2*q/x2;      
    
    % Calculate: q absorb (Equation 3.25)
    qabs = qext-qsca;   
    
    % Calculate: q backscatter (From pp. 122)
    fn =(f(1,:)-f(2,:)).*cn;
    gn =(-1).^n;
    f(3,:)= fn.*gn;
    q = sum(f(3,:));
    qb = q*q'/x2;       
    
    % Calcuate: q ratio
    qratio = qb/qsca;   
    
    % Calculate: asymmetry factor 
    asy1 = c1n.*(anp.*g1(1,:)+anpp.*g1(2,:)+bnp.*g1(3,:)+bnpp.*g1(4,:));
    asy2 = c2n.*(anp.*bnp+anpp.*bnpp);
    asy = 4/x2*sum(asy1+asy2)/qsca;
    
    % Outputs results / Stores Values
    if isreal(xi) == 0 % Differentiates plotting & non-plotting situations
    disp(['     qext  ' '    qsca ' '     qabs ' '      qb  ' '     asy '...
        '     qratio'])
    disp([qext, qsca, qabs, qb, asy, qratio])
    end
    res(k,:) = [qext, qsca, qabs, qb, asy, qratio]; % Stores values
end
end

      
    %%% Plotting Efficiencies And Asymmetry Factors %%%
if isreal(xi) == 0 % Differentiates plotting & non-plotting situations
else
    plot(xsize(:),res(:,1:6)) % Plots efficiencies against size parameter
    xlabel('x'); grid on;
    title(sprintf('Mie Efficiencies, m = %g + %gi', real(m), imag(m)));
    legend('Q_{ext}','Q_{sca}','Q_{abs}','Q_{b}','cos(\theta)',...
        'Q_{b}/Q_{sca}');        
end


    %%% Calculates angle functions pi_n and tau_n %%%
if isreal(xi)
    break;
else
    if isreal(theta) == 0   % For when theta is not specified 
        break;
    else
        u = cos(theta*pi/180);  % -1 <= u <= 1: uses scattering angle
        p(1) = 1;       % Initial pi angle function term
        t(1) = u;       % Initial tau angle fucntion term
        p(2) = 3*u;     % Secondary pi angle function term
        t(2) = 3*cos(2*acos(u));    % Secondary tau angle function term
        for n1 = 3:nmax,    % Calculates angle functions from n = 3 to nmax
            p1 = (2*n1-1)./(n1-1).*p(n1-1).*u;
            p2 = n1./(n1-1).*p(n1-2);
            p(n1) = p1-p2;      % pi_n: array of pi_n angle functions
            t1 = n1*u.*p(n1);
            t2 = (n1+1).*p(n1-1);
            t(n1) = t1-t2;      % tau_n: array of tau_n angle functions
        end
        % Outputs results
        disp(['   pi_n' '  tau_n'])
        disp([p', t'])
    end
end