% Scattering Coefficients for Coated Particle
Qsca=[];
Qext=[];
Qabs=[];
wr=[0.1:0.1:10];
for k=2*pi./wr
tempQsca=0;
tempQext=0;
tempQabs=0;

for n=1:5;

a=1;
b=6.265;
m1=1.33+1e-8i;
m2=1;

x=a*k;
y=b*k;

phi_y=y*besselj(n+1/2,y);
phi_m1x=m1*x*besselj(n+1/2,m1*x);
phi_m2x=m2*x*besselj(n+1/2,m2*x);
phi_m2y=m2*y*besselj(n+1/2,m2*y);
phiD_y=besselj(n+1/2,y)+n/(2*n+1)*y*besselj(n-1+1/2,y)-(n+1)/(2*n+1)*y*besselj(n+1+1/2,y);
phiD_m2y=besselj(n+1/2,m2*y)+n/(2*n+1)*m2*y*besselj(n-1+1/2,m2*y)-(n+1)/(2*n+1)*m2*y*besselj(n+1+1/2,m2*y);
phiD_m1x=besselj(n+1/2,m1*x)+n/(2*n+1)*m1*x*besselj(n-1+1/2,m1*x)-(n+1)/(2*n+1)*m1*x*besselj(n+1+1/2,m1*x);
phiD_m2x=besselj(n+1/2,m2*x)+n/(2*n+1)*m2*x*besselj(n-1+1/2,m2*x)-(n+1)/(2*n+1)*m2*x*besselj(n+1+1/2,m2*x);

ksi_y=y*besselh(n+1/2,1,y);
ksiD_y=besselh(n+1/2,1,y)+n/(2*n+1)*y*besselh(n-1+1/2,1,y)-(n+1)/(2*n+1)*y*besselh(n+1+1/2,1,y);

kai_m2y=-m2*y*bessely(n+1/2,m2*y);
kai_m2x=-m2*x*bessely(n+1/2,m2*x);
kaiD_m2y=-bessely(n+1/2,m2*y)-n/(2*n+1)*m2*y*bessely(n+1/2,m2*y)+(n+1)/(2*n+1)*m2*y*bessely(n+1+1/2,m2*y);
kaiD_m2x=-bessely(n+1/2,m2*x)-n/(2*n+1)*m2*x*bessely(n+1/2,m2*x)+(n+1)/(2*n+1)*m2*x*bessely(n+1+1/2,m2*x);

An=(m2*phi_m2x*phiD_m1x-m1*phiD_m2x*phi_m1x)/(m2*kai_m2x*phiD_m1x-m1*kaiD_m2x*phi_m1x);
Bn=(m2*phi_m1x*phiD_m2x-m1*phi_m2x*phiD_m1x)/(m2*kaiD_m2x*phi_m1x-m1*phiD_m1x*kai_m2x);
an=(phi_y*(phiD_m2y-An*kaiD_m2y)-m2*phiD_y*(phi_m2y-An*kai_m2y))/(ksi_y*(phiD_m2y-An*kaiD_m2y)-m2*ksiD_y*(phi_m2y-An*kai_m2y));
bn=(m2*phi_y*(phiD_m2y-Bn*kaiD_m2y)-phiD_y*(phi_m2y-Bn*kai_m2y))/(m2*ksi_y*(phiD_m2y-Bn*kaiD_m2y)-ksiD_y*(phi_m2y-Bn*kai_m2y));

tempQsca=tempQsca+(2*n+1)*(an*conj(an)+bn*conj(bn))*2/(k*b)^2;
tempQext=tempQext+(2*n+1)*real(an+bn)*2/(k*b)^2;
tempQabs=tempQext-tempQsca;
end

Qsca=[Qsca tempQsca];
Qext=[Qext tempQext];
Qabs=[Qabs tempQabs];
end

plot(wr,Qabs,'r')