Code covered by the BSD License

# Scattering from hexagonal PEC cylinder of infinite length

### Saira Bashir (view profile)

This code calculates scattered field from a an infinite length hexagonal PEC cylinder

Sairanemtassignment.m
```%%%%%%%%%%%%%%%%%%%%%%Hexagonal cylinder code where patches of same edge areapproached using straight line equation%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
L=input('Enter Length per unit wavelength(electrical length) of 1 edge=')
N=input('Enter number of segments on 1 edge such that((10*L)<N)=')
k=2*3.14159;
delL=L/N ;%length of one patch
delphi=60*3.14159/(180*N); %Angular length of 1 patch in radians
delx=delL*cos(60*3.14159/180);
c=2*sqrt((L^2)-(L/2)^2);
for p=1:1:6
for a=1:1:N
mat(1,(a+(p-1)*N))=((a-1)+(p-1)*N)*delphi;%phi
if a==1 &&p==1
mat(2,(a+(p-1)*N))=L;
mat(3,(a+(p-1)*N))=0;
elseif a==1 &&p==2||a==1 &&p==6
mat(2,(a+(p-1)*N))=L/2;
elseif a==1 &&p==3||a==1 &&p==5
mat(2,(a+(p-1)*N))=-L/2;
elseif a==1 &&p==4
mat(2,(a+(p-1)*N))=-L;
else
if p==1||p==3
mat(2,a+(p-1)*N)= mat(2,(1+(p-1)*N))-delx*(a-1);
end
if p==4||p==6
mat(2,a+(p-1)*N)= mat(2,(1+(p-1)*N))+delx*(a-1);
end
if p==2
mat(2,a+(p-1)*N)= mat(2,(1+(p-1)*N))-delL*(a-1);
end
if p==5
mat(2,a+(p-1)*N)= mat(2,(1+(p-1)*N))+delL*(a-1);
end
end     %%%%%%x-components calculated

if a==1&&p==1
mat(3,a)=0;
elseif a==1&&p==2||a==1&&p==3
mat(3,(a+(p-1)*N))=sqrt(L^2-(L/2)^2);
elseif a==1&&p==4
mat(3,(a+(p-1)*N))=0;
elseif a==1&&p==5||a==1&&p==6
mat(3,(a+(p-1)*N))=-sqrt(L^2-(L/2)^2);
else
if p==1
m1=(sqrt(L^2-(L/2)^2)/-(L/2));
mat(3,a)=m1*mat(2,a)+c;
end
if p==2

m2=0;
mat(3,(a+(p-1)*N))=sqrt(L^2-(L/2)^2);
end
if p==3
m3=sqrt(L^2-(L/2)^2)/(L/2);
mat(3,(a+((p-1)*N)))=m3*mat(2,(a+(p-1)*N))+c;
end
if p==4
m4=-sqrt(L^2-(L/2)^2)/(L/2);
mat(3,(a+((p-1)*N)))=m4*mat(2,(a+(p-1)*N))-c;
end
if p==5
m5=0;
mat(3,(a+((p-1)*N)))=-sqrt(L^2-(L/2)^2);
end
if p==6
m6=sqrt(L^2-(L/2)^2)/(L/2);
mat(3,(a+((p-1)*N)))=(m6*mat(2,(a+(p-1)*N)))-c;
end
end

end
end
% % % %       %%%%%%%% to check geometry of cylinder(Uncomment)%%%%%%%%%
% for k=1:1:(N*6)
%    xaxis(1,k)=mat(2,k);
%     yaxis(1,k)=mat(3,k);
% end
%    plot(xaxis,yaxis)
%%%%%%%%%%%%%%%%%Rmn formation%%%%%%%
for m=1:1:N*6
for n=1:1:N*6
if m==n
z(m,n)=delL*(1-i*(2/3.14159)*log(1.781*k*delL/(4*2.719)));
else
delrho=sqrt(((mat(2,m)-mat(2,n))^2+(mat(3,m)-mat(3,n))^2));
z(m,n)=delL*besselh(0,2,k*delrho);
end
end
end
rmn=inv(z);
%%%%%%%Exitation vector%%%%%%%%%%%%
phiI=input('Enter incident angle in degrees=')
for m=1:1:N*6
b(m,1)=4*exp(i*k*(((mat(2,m)*cos(phiI*3.14159/180))+(mat(3,m)*sin(phiI*3.14159/180)))));
end
%%%%%current calculation%%%
current=(rmn*b);
% %%%%%%%%%%%%PLOT%%%%%%%%%%
r=1:1:N*6;
figure(3)
plot(r,abs(current'))
% figure(4)
% polar(mat(1,:),current')
```