function [fn, gn, vn, wn] = expcoeff_mie_strat_int( an, bn, x, m )
%EXPCOEFF_MIE_STRAT_INT Calculates the expansion coefficients of the 
%   internal fields for a stratified sphere.
%
%   [FN,GN,VN,WN] = EXPCOEFF_MIE_STRAT(AN,BN,X,M) calculates the internal 
%   fields expansion coefficients FN, GN, VN and WN for a stratified sphere
%   with size parameters X and relative refractive indices M. The
%   corresponding scattered field expansion coefficients are given in AN
%   and BN.
%
%   The calculation is based on the book of Bohren and Huffman [1]:
%   [1] Bohren, C. F. and Huffman, D. R., Absorption and scattering of
%       light by small particles, Wiley-Interscience, New York, 1998.
%
%   Copyright 2012 Jan Schäfer, Institut für Lasertechnologien (ILM)
%   Author: Jan Schäfer (jan.schaefer@ilm.uni-ulm.de)
%   Organization: Institut für Lasertechnologien in der Medizin und
%       Meßtechnik an der Universität Ulm (http://www.ilm-ulm.de)

%% Calculate truncation number
M = numel(an);
N = numel(x);

%% Calculate expansion coefficients
fn = zeros(M,N);
gn = zeros(M,N);
vn = zeros(M,N);
wn = zeros(M,N);

for ilay=N:-1:1
    for n=1:M
        if ilay == N
            Sx = ricbesj(n, x(ilay));
            dSx = dricbesj(n, x(ilay));
            Smx = ricbesj(n, m(ilay)*x(ilay));
            dSmx = dricbesj(n, m(ilay)*x(ilay));
            xix = ricbesh(n, x(ilay));
            dxix = dricbesh(n, x(ilay));
            cmx = ricbesy(n, m(ilay)*x(ilay));
            dcmx = dricbesy(n, m(ilay)*x(ilay));
            denom = cmx*dSmx - dcmx*Smx;
            fn(n,ilay) = ((Sx - bn(n)*xix)*dcmx*m(ilay) + ...
                (bn(n)*dxix - dSx)*cmx)/-denom;
            gn(n,ilay) = ((dSx - an(n)*dxix)*cmx*m(ilay) + ...
                (an(n)*xix - Sx)*dcmx)/denom;
            vn(n,ilay) = ((dSx - bn(n)*dxix)*Smx + ...
                (bn(n)*xix - Sx)*dSmx*m(ilay))/denom;
            wn(n,ilay) = ((Sx - an(n)*xix)*dSmx + ...
                (an(n)*dxix - dSx)*Smx*m(ilay))/-denom;
        elseif ilay == 1
            Smp1x = ricbesj(n, m(ilay+1)*x(ilay));
            Smx = ricbesj(n, m(ilay)*x(ilay));
            cmp1x = ricbesy(n, m(ilay+1)*x(ilay));
            fn(n,ilay) = (fn(n,ilay+1)*m(ilay)*Smp1x - ...
                vn(n,ilay+1)*m(ilay)*cmp1x)/m(ilay+1)/Smx;
            gn(n,ilay) = (gn(n,ilay+1)*Smp1x - ...
                wn(n,ilay+1)*cmp1x)/Smx;
            
            % dSmx = dricbesj(n, m(ilay)*x(ilay));
            % dSmp1x = dricbesj(n, m(ilay+1)*x(ilay));
            % dcmp1x = dricbesy(n, m(ilay+1)*x(ilay));
            % fn(n,ilay) = (fn(n,ilay+1)*dSmp1x - ...
            %     vn(n,ilay+1)*dcmp1x)/dSmx;
            % gn(n,ilay) = (gn(n,ilay+1)*m(ilay)*dSmp1x - ...
            %     wn(n,ilay+1)*m(ilay)*dcmp1x)/m(ilay+1)/dSmx;
            
        else %ilay
            Smp1x = ricbesj(n, m(ilay+1)*x(ilay));
            dSmp1x = dricbesj(n, m(ilay+1)*x(ilay));
            Smx = ricbesj(n, m(ilay)*x(ilay));
            dSmx = dricbesj(n, m(ilay)*x(ilay));
            cmp1x = ricbesy(n, m(ilay+1)*x(ilay));
            dcmp1x = dricbesy(n, m(ilay+1)*x(ilay));
            cmx = ricbesy(n, m(ilay)*x(ilay));
            dcmx = dricbesy(n, m(ilay)*x(ilay));
            denom = m(ilay+1)*(Smx*dcmx - dSmx*cmx);
            fn(n,ilay) = (fn(n,ilay+1)*(m(ilay)*Smp1x*dcmx - ...
                m(ilay+1)*dSmp1x*cmx) + ...
                vn(n,ilay+1)*(m(ilay+1)*dcmp1x*cmx - ...
                m(ilay)*cmp1x*dcmx))/denom;
            gn(n,ilay) = (gn(n,ilay+1)*(m(ilay+1)*Smp1x*dcmx - ...
                m(ilay)*dSmp1x*cmx) + ...
                wn(n,ilay+1)*(m(ilay)*dcmp1x*cmx - ...
                m(ilay+1)*cmp1x*dcmx))/denom;
            vn(n,ilay) = (fn(n,ilay+1)*(m(ilay)*Smp1x*dSmx - ...
                m(ilay+1)*dSmp1x*Smx) + ...
                vn(n,ilay+1)*(m(ilay+1)*dcmp1x*Smx - ...
                m(ilay)*cmp1x*dSmx))/denom;
            wn(n,ilay) = (gn(n,ilay+1)*(m(ilay+1)*Smp1x*dSmx - ...
                m(ilay)*dSmp1x*Smx) + ...
                wn(n,ilay+1)*(m(ilay)*dcmp1x*Smx - ...
                m(ilay+1)*cmp1x*dSmx))/denom;
        end %ilay
    end %for n=1:M
end %for ilay=N:-1:1

end