image thumbnail

Scattering from hexagonal PEC cylinder of infinite length

by

 

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')

Contact us