TE-Wave Scattering from a dielectric cylinder a.k.a Richmond's article

   clc;
clear all;
close all;
format long;

N=6; Nphi=100; c=3e8; deltaphi=2*pi/Nphi; f=3e8; lambda=c/f; R=0.5*lambda; % Radius of the dielectric cylinder. w=2*pi*f; %eps0=(1e-9)/(36*pi); eps0=8.854187e-12; eps=4; mu0=4*pi*10^(-7); mur=1; mu=mur*mu0; eta0=sqrt(mu0/eps0); %k0=2*pi/lambda; k0=w*sqrt(eps0*mu0); k1=k0*sqrt(eps); k=k0; deltax=R/N; deltay=R/N; an=sqrt(deltax*deltay/pi); phi1=pi:deltaphi:-pi; phi=0:deltaphi:2*pi; E0=1; H0=1/eta0; %Eix=0;

x=-R:deltax:R; y=-R:deltay:R;

q=1; for m=1:2*N+1 for n=1:2*N+1 if sqrt((x(m)^2)+(y(n)^2))<=R x1(q)=x(m); y1(q)=y(n); q=q+1;

      end
   end
end
%  scatter(x1, y1,'+')
%  figure

Eix=zeros(1,length(y1)); Eiy=E0*exp(1j*k0.*(x1));

Eiy=Eiy.'; Eix=Eix.';

for m=1:length(x1) for n=1:length(y1) rho=sqrt(((x1(m)-x1(n))^2)+((y1(m)-y1(n))^2)); Kp=(1j*pi*an*besselj(1,k*an)*(eps-1))/(2*rho^3);

        if m==n
            A(m,n)=1+(eps-1)*(0.25*1j*pi*k*an*besselh(1,2,k*an)+1);
            B(m,n)=0;
            C(m,n)=0;
            D(m,n)=A(m,n);
        else
A(m,n)=Kp*(k*rho*((y1(m)-y1(n))^2)*besselh(0,2,k*rho)+((x1(m)-x1(n))^2-(y1(m)-y1(n))^2)*besselh(1,2,k*rho));
B(m,n)=Kp*(x1(m)-x1(n))*(y1(m)-y1(n))*(2*besselh(1,2,k*rho)-k*rho*besselh(0,2,k*rho));
C(m,n)=B(m,n);
D(m,n)=Kp*(k*rho*((x1(m)-x1(n))^2)*besselh(0,2,k*rho)+((y1(m)-y1(n))^2-(x1(m)-x1(n))^2)*besselh(1,2,k*rho));

        end
    end
    m
end

coeffm=[A, B C,D]; Ei=[Eix Eiy];

E=(eye(size(coeffm))-coeffm)\(Ei);

Exn=E(1:length(E)/2); Eyn=E((length(E)/2)+1:length(E));

Jx=1j*w*eps0*(eps-1).*Exn; Jy=1j*w*eps0*(eps-1).*Eyn; Jx=Jx.'; Jy=Jy.';

rhog=15*lambda; % Observation radius from the center of cylinder. xg=rhog*cos(phi); yg=rhog*sin(phi);

for m=1:length(phi); Exss(m)=0; Eyss(m)=0; Esf(m)=0; Esx(m)=0; Esy(m)=0; for n=1:length(x1) rhogn=sqrt((xg(m)-x1(n))^2+(yg(m)-y1(n))^2); % Observation radius from center of the nth cell at mth angel. xgn=rhogn*cos(phi(m)); ygn=rhogn*sin(phi(m)); K=-(pi*an*besselj(1,k*an))/(2*w*eps0*rhogn^3); Exs(m,n)=K*((k*rhogn*(ygn^2)*besselh(0,2,k*rhogn)+(xgn^2-ygn^2)*besselh(1,2,k*rhogn))*Jx(n)+xgn*ygn*(2*besselh(1,2,k*rhogn)-k*rhogn*besselh(0,2,k*rhogn))*Jy(n)); Eys(m,n)=K*(xgn*ygn*(2*besselh(1,2,k*rhogn)-k*rhogn*besselh(0,2,k*rhogn))*Jx(n)+(k*rhogn*(xgn^2)*besselh(0,2,k*rhogn)+(ygn^2-xgn^2)*besselh(1,2,k*rhogn))*Jy(n)); Exss(m)=Exss(m)+Exs(m,n); % Exs, x component of the scattered electric field Eyss(m)=Eyss(m)+Eys(m,n); % Eys, y component of the scattered electric field Esf(m)=Esf(m)+((1j*sqrt((1j*pi*k)/(2*rhog))*exp(1j*k*rhog))*(eps-1)*an*besselj(1, k*an)*(Exn(n)*sin(phi(m))-Eyn(n)*cos(phi(m)))*exp(1j*k*(xg(m)*cos(phi(m))+yg(m)*sin(phi(m))))); % MoM far field approx.

    end
    m
end

 Es=(Exss.*sin(phi)-Eyss.*cos(phi));
% ed1=abs(Exss.*sin(phi))+abs(Eyss.*cos(phi));

plot(rad2deg(phi), abs(Es)) hold on

% plot(rad2deg(phi), abs(Esf), '--k') % hold on

%%%%%%% Analytical Solution %%%%%%%

     for m=1:length(phi);
        Esanphi(m)=0;
        for n=-30:30
%             num=besjp(n,k0*R)*besselj(n,k1*R)-sqrt(mu/eps)*besselj(n,k0*R)*besjp(n, k1*R);
%             denum=(sqrt(mu/eps)*besjp(n, k1*R)*besselh(n,2,k0*R)-besselj(n,k1*R)*besh2p(n, k0*R));
%             ac=(j^(-n))*num/denum;
            numerator=(besjp(n,k0*R)*besselj(n,k1*R))-sqrt(mur/eps)*besselj(n,k0*R)*besjp(n,k1*R);
            denumerator=sqrt(mur/eps)*besjp(n,k1*R)*besselh(n,2,k0*R)-besselj(n,k1*R)*besh2p(n,k0*R);
            ac=(1j^(-n))*numerator/denumerator;
            %Esanphi(m)=Esanphi(m)+(E0*ac*besselh(n,2,k0*rhog)*exp(1j*n*phi(m)));
            Esanphi(m)=Esanphi(m)+((-H0)*k0/(1j*w*eps0))*ac*besh2p(n,k0*rhog)*exp(1j*n*phi(m));
        end
       m
    end

    plot(rad2deg(phi), abs(Esanphi), '--+r')
    title(['N=',num2str(N),';  ','Epsilon=',num2str(eps),';  ','Radius of the dielectric cylinder=',num2str(R/lambda),'*lambda',';  ','Radius of the observation=',num2str(rhog/lambda),'*lambda'])
    xlabel('Observation Angle (/phi)')
    ylabel('Scattered Electric Field (E_s)');
    legend('MoM Solution','Analitycal Solution')
    grid

hi, this my code. I don't know where I am wrong. The MoM solution should be the same as analytical solution. If anyone understands and helps me , I'll appreciate. Thx for your time.