```function [fnp, gnn, vnp, wnn] = expcoeff_cyl_strat_int( ann, bnp, r, k, m )
%EXPCOEFF_CYL_STRAT_INT Calculates the internal expansion coefficients for
%   the scattering by a stratified infinite cylinder at perpendicular
%   incidence.
%
%   [FNP,GNN,VNP,WNN] = EXPCOEFF_CYL_STRAT_INT(ANN,BNP,R,K,M) calculates
%   the expansion coefficients (FNP,GNN,VNP,WNN) for a stratified cylinder
%   with size parameters X and relative refractive indices M at
%   perpendicular incidence with wavenumber K.
%
%   The calculation is based on the book of Kerker [5]:
%   [5] Kerker, M., The scattering of light and other electromagnetic
%
%   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)

N = numel(r);   % Number of layers
MM = numel(ann);
M = (MM - 1)/2;

%% Calculate expansion coefficients
fnp = zeros(MM,N);
gnn = zeros(MM,N);
vnp = zeros(MM,N);
wnn = zeros(MM,N);

for ilay=N:-1:1
if ilay == N
mfrac = 1/m(ilay);
rho_lp1 = k*r(ilay);
lp1 = k;
else
mfrac =  m(ilay+1)/m(ilay);
lp1 = k*m(ilay+1);
rho_lp1 = lp1*r(ilay);
end
rho_l = k*m(ilay)*r(ilay);
for n=-M:M
if ilay == N
u_lp1 = besselj(n, rho_lp1) - bnp(n+M+1)*besselh(n, rho_lp1);
u_lp1_dr = lp1*(dbesselj(n, rho_lp1) - ...
bnp(n+M+1)*dbesselh(n, rho_lp1));
v_lp1 = besselj(n, rho_lp1) - ann(n+M+1)*besselh(n, rho_lp1);
v_lp1_dr = lp1*(dbesselj(n, rho_lp1) - ...
ann(n+M+1)*dbesselh(n, rho_lp1));
else
u_lp1 = fnp(n+M+1,ilay+1)*besselj(n, rho_lp1) - ...
vnp(n+M+1,ilay+1)*besselh(n, rho_lp1);
u_lp1_dr = lp1*(fnp(n+M+1,ilay+1)*dbesselj(n, rho_lp1) - ...
vnp(n+M+1,ilay+1)*dbesselh(n, rho_lp1));
v_lp1 = gnn(n+M+1,ilay+1)*besselj(n, rho_lp1) - ...
wnn(n+M+1,ilay+1)*besselh(n, rho_lp1);
v_lp1_dr = lp1*(gnn(n+M+1,ilay+1)*dbesselj(n, rho_lp1) - ...
wnn(n+M+1,ilay+1)*dbesselh(n, rho_lp1));
end
if ilay > 1
denom = dbesselh(n, rho_l)*besselj(n, rho_l) - ...
besselh(n, rho_l)*dbesselj(n, rho_l);

vnp(n+M+1,ilay) = mfrac/m(ilay)/k*(k*m(ilay)*...
dbesselj(n, rho_l)*u_lp1 - ...
besselj(n, rho_l)*u_lp1_dr)/denom;
fnp(n+M+1,ilay) = mfrac/m(ilay)/k*(k*m(ilay)*...
dbesselh(n, rho_l)*u_lp1 - ...
besselh(n, rho_l)*u_lp1_dr)/denom;
wnn(n+M+1,ilay) = 1/m(ilay)/k*(lp1*mfrac*...
dbesselj(n, rho_l)*v_lp1 - ...
besselj(n, rho_l)*v_lp1_dr)/denom;
gnn(n+M+1,ilay) = 1/m(ilay)/k*(lp1*mfrac*...
dbesselh(n, rho_l)*v_lp1 - ...
besselh(n, rho_l)*v_lp1_dr)/denom;
else
fnp(n+M+1,ilay) = mfrac*u_lp1/besselj(n, rho_l);
gnn(n+M+1,ilay) = mfrac.^2*v_lp1/besselj(n, rho_l);
end %ilay
end %for n=1:M
end %for ilay=N:-1:1

end

```