circle_mesh.m

In this program we calculate the scattering field of a infinite dielectric cylinder by two methods.

%               Advanced Engineering Electromagnetic
%                   Pouya Dastmalchi  8600844

%                Triangular mesh plot for a Circule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xx,yy,S] = circle_mesh(a,dp,dr)  
% close all
% clear all
% clc
% c = 3e8;
% f = 1e9;
% lambda = c/f;
% a = lambda/4;
% a=1;
% dp=pi/6;
% dr=a/3;
z=1;
X(1)=0;
Y(1)=0;

for k=1:(a/dr)
    Np(k)=2*k-1;
end

N=(sum(Np(1:end)))*2*pi/dp;


for h=1:2*pi/dp
  T(h,1,1)=1;
end

for n=1:2*pi/dp
    for m=2:a/dr+1
        
        for k=1:m
        z=z+1;    
        x(n,m,k)=(m-1)*dr*cos(((n-1)*dp)+((k-1)*dp/(m-1)));
        y(n,m,k)=(m-1)*dr*sin(((n-1)*dp)+((k-1)*dp/(m-1)));
        X(z)=x(n,m,k);
        Y(z)=y(n,m,k);
        T(n,m,k)=z;
        end
    end
end

   for n=1:2*pi/dp
    for m=1:(a/dr)
        for k=1:m
            l=m*(m-1)/2+k+(n-1)*sum(Np(1:end));
            tri(l,1)=T(n,m+1,k);
            tri(l,2)=T(n,m+1,k+1);
            tri(l,3)=T(n,m,k);            
        end
    end
    
    for m=1:(a/dr)-1
        for k=1:m
            l=l+1;
            tri(l,1)=T(n,m+1,k);
            tri(l,2)=T(n,m+1,k+1);
            tri(l,3)=T(n,m+2,k+1);
        end
    end    
   end 

for k=1:N
    xx(k)=(X(tri(k,1))+X(tri(k,2))+X(tri(k,3)))/3;
    yy(k)=(Y(tri(k,1))+Y(tri(k,2))+Y(tri(k,3)))/3;
    
    L1 = sqrt((X(tri(k,1))-X(tri(k,2)))^2+(Y(tri(k,1))-Y(tri(k,2)))^2);
    L2 = sqrt((X(tri(k,3))-X(tri(k,2)))^2+(Y(tri(k,3))-Y(tri(k,2)))^2);
    L3 = sqrt((X(tri(k,3))-X(tri(k,1)))^2+(Y(tri(k,3))-Y(tri(k,1)))^2);
    P = (L1+L2+L3)/2;
    S(k) = sqrt(P*(P-L1)*(P-L2)*(P-L3));
%     O(k,1)=xx(k);
%     O(k,2)=yy(k);    
end


% Z=zeros(1,length(X));
% figure(1)
% trimesh(tri,X,Y,Z)

% figure(2)
% plot(X,Y,'k.')
% hold on
% triplot(tri,X,Y,'b')
% hold on
% plot(xx,yy,'r.')
% axis equal

%___________________________________________________________________
% c = 3e8;
% f = 1e9;
% lambda = c/f;
% w = 2*pi*f;
% a = lambda;
% phi0 = pi/3;
% G = 1.781072418;
% I0 = 1;
% eps0 = (1/(36*pi))*10^(-9);
% eps = 3;
% eps1 = 3*eps0;
% neta = 377;
% B = 2*pi/lambda;
% u = 4*pi*10^(-7);
% r = 4*lambda;
% phi = 0:0.05:2*pi;
% X = r*cos(phi);
% Y = r*sin(phi);
% La = length(xx);
% 
% for q=1:La 
%     for p=1:La
%         xm = xx(p);
%         ym = yy(p);
%         xn = xx(q);
%         yn = yy(q);
%         if p~=q
%             A(p,q) = (j*pi*B*a/2)*(eps-1)*besselj(1,B*a)*besselh(0,2,B*sqrt((xm-xn)^2+(ym-yn)^2));
%         elseif p==q
%             A(p,q) = 1+(eps-1)*(j/2)*(pi*B*a*besselh(1,2,B*a)-2*j);
%         end
%     end
%             B1(q,1) = exp(j*B*(xn*cos(phi0)+yn*sin(phi0)));
% % 
% end
% % 
% ALFA = [];
% ALFA = A\B1;
% J = j*w*(eps-eps0)*ALFA;
% 
% for tt = 1:length(phi)
%     for gg = 1:La
%         R = sqrt((X(tt)-xx(gg))^2+(Y(tt)-yy(gg))^2);
%         Ez(gg) = -J(gg)*a^2*(w*u/4)*besselh(0,2,B*R);
%     end
%     EZ(tt) = sum(Ez(1:end));
% end
% figure(2)
% RCS =2*pi*r*(abs(EZ).^2)/lambda;
% hold on
% plot(phi*180/pi,RCS,'-')