Antenna and EM Modeling with MATLAB

Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/appendix.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_1.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_10.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_2.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_3.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_4.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_5.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_6.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_7.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_8.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/chapter_9.pdf Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/impmet.mfunction [Z]= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %IMPMET Standard impedance matrix (metal surface) % % Returns the complex impedance matrix [EdgesTotal x EdgesTotal] % Uses 9 integration points for every triangle % (barycentric subdivision) % % The impedance matrix is calculated as a sum of the contributions % due to separate triangles (similar to the "face-pair" method). % See Appendix B for a detailed algorithm. % % A 9-point quadrature is used for all integrals, including % the self-coupling terms. The alternative source code with % the analytical approximation of the self-coupling terms % is given in Appendix B. The difference between two methods % is not significant. % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 2 %Memory allocation Z =zeros (EdgesTotal,EdgesTotal)+j*zeros(EdgesTotal,EdgesTotal); %Loop over integration triangles for p=1:TrianglesTotal Plus =find(TrianglePlus-p==0); Minus =find(TriangleMinus-p==0); D=Center_-repmat(Center(:,p),[1 9 TrianglesTotal]); %[3 9 TrianglesTotal] R=sqrt(sum(D.*D)); %[1 9 TrianglesTotal] g=exp(-K*R)./R; %[1 9 TrianglesTotal] gP=g(:,:,TrianglePlus); %[1 9 EdgesTotal] gM=g(:,:,TriangleMinus); %[1 9 EdgesTotal] Fi=sum(gP)-sum(gM); %[1 1 EdgesTotal] ZF= FactorFi.*reshape(Fi,EdgesTotal,1); %[EdgesTotal 1] for k=1:length(Plus) n=Plus(k); RP=repmat(RHO__Plus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1+ZF); end for k=1:length(Minus) n=Minus(k); RP=repmat(RHO__Minus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1-ZF); end end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/matlabcompiler/impmet.mfunction [Z]= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %IMPMET Standard impedance matrix (metal surface) % % Returns the complex impedance matrix [EdgesTotal x EdgesTotal] % Uses 9 integration points for every triangle % (barycentric subdivision) % % The impedance matrix is calculated as a sum of the contributions % due to separate triangles (similar to the "face-pair" method). % See Appendix B for a detailed algorithm. % % A 9-point quadrature is used for all integrals, including % the self-coupling terms. The alternative source code with % the analytical approximation of the self-coupling terms % is given in Appendix B. The difference between two methods % is not significant. % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 2 %Memory allocation Z =zeros (EdgesTotal,EdgesTotal)+j*zeros(EdgesTotal,EdgesTotal); %Loop over integration triangles for p=1:TrianglesTotal Plus =find(TrianglePlus-p==0); Minus =find(TriangleMinus-p==0); D=Center_-repmat(Center(:,p),[1 9 TrianglesTotal]); %[3 9 TrianglesTotal] R=sqrt(sum(D.*D)); %[1 9 TrianglesTotal] g=exp(-K*R)./R; %[1 9 TrianglesTotal] gP=g(:,:,TrianglePlus); %[1 9 EdgesTotal] gM=g(:,:,TriangleMinus); %[1 9 EdgesTotal] Fi=sum(gP)-sum(gM); %[1 1 EdgesTotal] ZF= FactorFi.*reshape(Fi,EdgesTotal,1); %[EdgesTotal 1] for k=1:length(Plus) n=Plus(k); RP=repmat(RHO__Plus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1+ZF); end for k=1:length(Minus) n=Minus(k); RP=repmat(RHO__Minus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1-ZF); end end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/matlabcompiler/rwg3.mfunction [] = rwg3 %RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =3e8; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/bowtie.mfunction pdemodel [pde_fig,ax]=pdeinit; pdetool('appl_cb',1); set(ax,'DataAspectRatio',[1 1.0000000000000002 1]); set(ax,'PlotBoxAspectRatio',[1.4999999999999998 0.99999999999999978 10]); set(ax,'XLim',[-0.14999999999999999 0.14999999999999999]); set(ax,'YLim',[-0.10000000000000001 0.10000000000000001]); set(ax,'XTickMode','auto'); set(ax,'YTickMode','auto'); pdetool('gridon','on'); % Geometry description: pdepoly([ -0.099000000000000005,... -0.006000000000000001,... -0.006000000000000001,... -0.099000000000000005,... 0.099000000000000005,... 0.006000000000000001,... 0.006000000000000001,... 0.099000000000000005,... ],... [ 0.10000000000000001,... 0.007000000000000001,... -0.007000000000000001,... -0.10000000000000001,... -0.10000000000000001,... -0.007000000000000001,... 0.007000000000000001,... 0.10000000000000001,... ],... 'P1'); set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','P1') % Boundary conditions: pdetool('changemode',0) pdesetbd(8,... 'dir',... 1,... '1',... '0') pdesetbd(7,... 'dir',... 1,... '1',... '0') pdesetbd(6,... 'dir',... 1,... '1',... '0') pdesetbd(5,... 'dir',... 1,... '1',... '0') pdesetbd(4,... 'dir',... 1,... '1',... '0') pdesetbd(3,... 'dir',... 1,... '1',... '0') pdesetbd(2,... 'dir',... 1,... '1',... '0') pdesetbd(1,... 'dir',... 1,... '1',... '0') % Mesh generation: setuprop(pde_fig,'trisize',0.014999999999999999); setuprop(pde_fig,'Hgrad',1.3); setuprop(pde_fig,'refinemethod','regular'); pdetool('initmesh') % PDE coefficients: pdeseteq(1,... '1.0',... '0.0',... '10.0',... '1.0',... '0:10',... '0.0',... '0.0',... '[0 100]') setuprop(pde_fig,'currparam',... ['1.0 ';... '0.0 ';... '10.0';... '1.0 ']) % Solve parameters: setuprop(pde_fig,'solveparam',... str2mat('0','1000','10','pdeadworst',... '0.5','longest','0','1E-4','','fixed','Inf')) % Plotflags and user data strings: setuprop(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]); setuprop(pde_fig,'colstring',''); setuprop(pde_fig,'arrowstring',''); setuprop(pde_fig,'deformstring',''); setuprop(pde_fig,'heightstring',''); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/bowtie.matp:[3x174 double array] t:[4x274 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/dipole.mfunction pdemodel [pde_fig,ax]=pdeinit; pdetool('appl_cb',1); set(ax,'DataAspectRatio',[1.5 1 1]); set(ax,'PlotBoxAspectRatio',[1 1 1]); set(ax,'XLim',[-1.5 1.5]); set(ax,'YLim',[-1 1]); set(ax,'XTickMode','auto'); set(ax,'YTickMode','auto'); pdetool('gridon','on'); % Geometry description: pderect([-0.0250000000000000003 0.0250000000000000003 1 -1],'R1'); set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1') % Boundary conditions: pdetool('changemode',0) pdesetbd(4,... 'dir',... 1,... '1',... '0') pdesetbd(3,... 'dir',... 1,... '1',... '0') pdesetbd(2,... 'dir',... 1,... '1',... '0') pdesetbd(1,... 'dir',... 1,... '1',... '0') % Mesh generation: setuprop(pde_fig,'trisize',0.125); setuprop(pde_fig,'Hgrad',1.05); setuprop(pde_fig,'refinemethod','regular'); pdetool('initmesh') % PDE coefficients: pdeseteq(1,... '1.0',... '0.0',... '10.0',... '1.0',... '0:10',... '0.0',... '0.0',... '[0 100]') setuprop(pde_fig,'currparam',... ['1.0 ';... '0.0 ';... '10.0';... '1.0 ']) % Solve parameters: setuprop(pde_fig,'solveparam',... str2mat('0','1000','10','pdeadworst',... '0.5','longest','0','1E-4','','fixed','Inf')) % Plotflags and user data strings: setuprop(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]); setuprop(pde_fig,'colstring',''); setuprop(pde_fig,'arrowstring',''); setuprop(pde_fig,'deformstring',''); setuprop(pde_fig,'heightstring',''); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/dipole.matp:[3x74 double array] t:[4x84 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/plate.m%PLATE Plate mesh in the xy plane - all Chapters % % Copyright 2002 AEMM. Revision 2002/04/10 clear all warning off W=2.00; %Plate width (along the x-axis) L=2.00; %Plate length (along the y-axis) Nx=20; %Discretization parameter (width) Ny=20; %Discretization parameter (length) %Set vertexes epsilon=1e-6; M=1; for i=1:Nx+1 for j=1:Ny+1 X(M)=-W/2+(i-1)/Nx*W; Y(M)=-L/2+(j-1)/Ny*L-epsilon*X(M); M=M+1; end end %Use Delaunay triangulation TRI = delaunay(X,Y,0); t=TRI'; t(4,:)=1; p=[X; Y; zeros(1,length(X))]; %Save result save plate p t; viewer plate Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/plate.matp:[3x81 double array] t:[4x128 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/platecoarse.matp:[3x81 double array] t:[4x128 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/platefine.matp:[3x289 double array] t:[4x512 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/slot.mfunction pdemodel [pde_fig,ax]=pdeinit; pdetool('appl_cb',1); set(ax,'DataAspectRatio',[1 1.125 1]); set(ax,'PlotBoxAspectRatio',[2 1.3333333333333333 1]); set(ax,'XLim',[-2 2]); set(ax,'YLim',[-1.5 1.5]); set(ax,'XTickMode','auto'); set(ax,'YTickMode','auto'); pdetool('gridon','on'); % Geometry description: pderect([-1.5 1.5 1.5 -1.5],'R1'); pderect([-0.029999999999999999 0.029999999999999999 1 0.029999999999999999],'R2'); pderect([-0.029999999999999999 0.029999999999999999 -0.029999999999999999 -1],'R3'); set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1+R2+R3') % Boundary conditions: pdetool('changemode',0) pdesetbd(6,... 'dir',... 1,... '1',... '0') pdesetbd(5,... 'dir',... 1,... '1',... '0') pdesetbd(2,... 'dir',... 1,... '1',... '0') pdesetbd(1,... 'dir',... 1,... '1',... '0') % Mesh generation: setuprop(pde_fig,'Hgrad',1.3); setuprop(pde_fig,'refinemethod','regular'); pdetool('initmesh') % PDE coefficients: pdeseteq(1,... '1.0',... '0.0',... '10.0',... '1.0',... '0:10',... '0.0',... '0.0',... '[0 100]') setuprop(pde_fig,'currparam',... ['1.0 ';... '0.0 ';... '10.0';... '1.0 ']) % Solve parameters: setuprop(pde_fig,'solveparam',... str2mat('0','1158','10','pdeadworst',... '0.5','longest','0','1E-4','','fixed','Inf')) % Plotflags and user data strings: setuprop(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]); setuprop(pde_fig,'colstring',''); setuprop(pde_fig,'arrowstring',''); setuprop(pde_fig,'deformstring',''); setuprop(pde_fig,'heightstring',''); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/slot.matp:[3x455 double array] t:[4x864 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/mesh/viewer.mfunction []=viewer(filename) %VIEWER Visualizes the structure - all Chapters % % Usage: viewer('plate.mat') or % viewer('plate') or % viewer plate % % Copyright 2002 AEMM. Revision 2002/03/05 Chapter 2 %Load the data load(filename); %Color=[0.75 0.75 0.75]; Color='g'; Cut=find (t(4,:)>1); t(:,Cut)=[]; for m=1:length(t) n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end g=fill3(X,Y,Z,Color); axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/rwg1.m%RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % Copyright 2002 AEMM. Revision 2002/03/09 Chapter 2 clear all tic; load('mesh/platefine'); [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) %All structures of this chapter have EdgeIndicator=2 EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =3e8; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/rwg4.m%RWG4 Solves MoM equations for the scattering problem % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % Direction of the incident signal in Cartesian coordinates % d(1:3); % Direction of the E-field in the incident plane wave % in Cartesian coordinates Pol(1:3); % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('impedance'); %Incident field %Example: d=[0 0 -1] means that the incident signal % is in the -z direction. %Plate - normal incidence d =[0 0 -1]; Pol =[1 0 0]; %Dipole - normal incidence %d =[0 0 -1]; %Pol =[0 1 0]; %Custom incidence (example) %d =[1 0 0] %Pol =[0 -0.0037-0.0055*j 0] k=omega/c_; kv=k*d; for m=1:EdgesTotal ScalarProduct=sum(kv.*Center(:,TrianglePlus(m))'); EmPlus =Pol.'*exp(-j*ScalarProduct); ScalarProduct=sum(kv.*Center(:,TriangleMinus(m))'); EmMinus=Pol.'*exp(-j*ScalarProduct); ScalarPlus =sum(EmPlus.* RHO_Plus(:,m)); ScalarMinus=sum(EmMinus.*RHO_Minus(:,m)); V(m)=EdgeLength(m)*(ScalarPlus/2+ScalarMinus/2); end tic; %System solution I=Z\V.'; toc %elapsed time FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','I','V','d','Pol'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter02/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index=find(t(4,:) 0 hpol = q; end if ~HoldFlag axis('equal');axis('off'); end % reset hold state if ~HoldFlag, set(CurrentAxes,'NextPlot',NextPlot); end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter03/sphere.matp:[3x252 double array] t:[3x500 double array] n:[3x500 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/efield1.m%EFIELD1 Radiated/scattered field at a point % The point is outside the metal surface % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % The following parameters need to be specified: % % Observation point ObservationPoint[X; Y; Z] (m) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 clear all %Load the data load('mesh2'); load('current'); k=omega/c_; K=j*k; for m=1:EdgesTotal Point1=Center(:,TrianglePlus(m)); Point2=Center(:,TriangleMinus(m)); DipoleCenter(:,m)=0.5*(Point1+Point2); DipoleMoment(:,m)=EdgeLength(m)*I(m)*(-Point1+Point2); end ObservationPoint=[0; 0; 100]; [E,H]=point(ObservationPoint,eta_,K,DipoleMoment,DipoleCenter); %find the sum of all dipole contributions EField=sum(E,2); HField=sum(H,2); %Common EField %Radiated/scattered electric field %(complex vector at a point, V/m) HField %Radiated/scattered magnetic field %(complex vector at a point, A/m) Poynting=0.5*real(cross(EField,conj(HField))) %Poynting vector (W/m^2) for radiated/scattered field W=norm(Poynting) %Radiation density (W/m^2) for radiated/scattered field U=norm(ObservationPoint)^2*W %Radiation intensity (W/unit solid angle) %Only scattering RCS=4*pi*(norm(ObservationPoint))^2*sum(EField.*conj(EField)); %Backscattering radar cross-section (scattering) Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/efield2.m%EFIELD2 Radiated/scattered field over a large sphere % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Uses the structure sphere.mat/sphere1.mat to display % radiation intensity distribution over the sphere surface. % The sphere doesn't intersect the radiating object. % % The following parameters need to be specified: % % Sphere radius (m) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 clear all %Load the data load('mesh2'); load('current'); load('sphere'); p=100*p; %sphere radius is 100 m k=omega/c_; K=j*k; for m=1:EdgesTotal Point1=Center(:,TrianglePlus(m)); Point2=Center(:,TriangleMinus(m)); DipoleCenter(:,m)=0.5*(Point1+Point2); DipoleMoment(:,m)=EdgeLength(m)*I(m)*(-Point1+Point2); end TotalPower=0; %Sphere series M=length(t); for m=1:M N=t(1:3,m); ObservationPoint=1/3*sum(p(:,N),2); [E,H]=point(ObservationPoint,eta_,K,DipoleMoment,DipoleCenter); ET=sum(E,2); HT=sum(H,2); Poynting(:,m)=0.5*real(cross(ET,conj(HT))); U(m)=(norm(ObservationPoint))^2*norm(Poynting(:,m)); Vector1=p(:,N(1))-p(:,N(2)); Vector2=p(:,N(3))-p(:,N(2)); Area =0.5*norm(cross(Vector1,Vector2)); TotalPower=TotalPower+norm(Poynting(:,m))*Area; %------------------------------ X(1:3,m)=[p(1,N)]'; Y(1:3,m)=[p(2,N)]'; Z(1:3,m)=[p(3,N)]'; end TotalPower GainLogarithmic =10*log10(4*pi*max(U)/TotalPower) GainLinear =4*pi*max(U)/TotalPower RadiationResistance =2*TotalPower/abs(GapCurrent)^2 FileName='gainpower.mat'; save(FileName, 'TotalPower','GainLogarithmic','GainLinear'); U=U/norm(U); C=repmat(U,3,1); h=fill3(X,Y,Z,C); colormap gray; axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/efield3.m%EFIELD3 2D Radiation patterns % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Additionally uses the value of TotalPower saved % in file gainpower.mat (script efield2.m) % % The following parameters need to be specified: % % Radius of the circle (m) R % Plane of the circle: [x y 0] or % [x 0 z] or % [0 y z] % Number of discretization points per % pattern NumPoints % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 clear all %Load the data load('mesh2'); load('current'); load('gainpower'); k=omega/c_; K=j*k; for m=1:EdgesTotal Point1=Center(:,TrianglePlus(m)); Point2=Center(:,TriangleMinus(m)); DipoleCenter(:,m)=0.5*(Point1+Point2); DipoleMoment(:,m)=EdgeLength(m)*I(m)*(-Point1+Point2); end NumPoints=100; R=100; %pattern in m for ii=1:NumPoints+1 phi(ii)=(ii-1)*pi/(NumPoints/2); y=R*cos(phi(ii)); z=R*sin(phi(ii)); ObservationPoint=[0 y z]'; [E,H]=point(ObservationPoint,eta_,K,DipoleMoment,DipoleCenter); ET=sum(E,2); HT=sum(H,2); Poynting=0.5*real(cross(ET,conj(HT))); W(ii)=norm(Poynting); U(ii)=(norm(ObservationPoint))^2*W(ii); end Polar_=10*log10(4*pi*U/TotalPower); GainLogarithmic=max(Polar_) %gain for the particular pattern! %This is the standard Matlab polar plot OFFSET=40; polar(phi,max(Polar_+OFFSET,0)); grid on; Title=strcat('Offset= ', num2str(OFFSET), ' dB'); title(Title); %This is Balanis' relative power: %Polar=10*log10(U/max(U)); OFFSET=40; polar(phi,Polar+OFFSET); grid on; %This is Galenski's polar plot: %Use with care: outputs an error if Polar is a constant function %polarhg(phi',Polar','rlim',[min(Polar) 10],'rtick',[-30 -20 -10 0 10],'tstep',90,'color','b'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/impmet.mfunction [Z]= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %IMPMET Standard impedance matrix (metal surface) % % Returns the complex impedance matrix [EdgesTotal x EdgesTotal] % Uses 9 integration points for every triangle % (barycentric subdivision) % % The impedance matrix is calculated as a sum of the contributions % due to separate triangles (similar to the "face-pair" method). % See Appendix B for a detailed algorithm. % % A 9-point quadrature is used for all integrals, including % the self-coupling terms. The alternative source code with % the analytical approximation of the self-coupling terms % is given in Appendix B. The difference between two methods % is not significant. % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 2 %Memory allocation Z =zeros (EdgesTotal,EdgesTotal)+j*zeros(EdgesTotal,EdgesTotal); %Loop over integration triangles for p=1:TrianglesTotal Plus =find(TrianglePlus-p==0); Minus =find(TriangleMinus-p==0); D=Center_-repmat(Center(:,p),[1 9 TrianglesTotal]); %[3 9 TrianglesTotal] R=sqrt(sum(D.*D)); %[1 9 TrianglesTotal] g=exp(-K*R)./R; %[1 9 TrianglesTotal] gP=g(:,:,TrianglePlus); %[1 9 EdgesTotal] gM=g(:,:,TriangleMinus); %[1 9 EdgesTotal] Fi=sum(gP)-sum(gM); %[1 1 EdgesTotal] ZF= FactorFi.*reshape(Fi,EdgesTotal,1); %[EdgesTotal 1] for k=1:length(Plus) n=Plus(k); RP=repmat(RHO__Plus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1+ZF); end for k=1:length(Minus) n=Minus(k); RP=repmat(RHO__Minus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1-ZF); end end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/mesh.matp:[3x678 double array] t:[4x788 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/monopole.m%MONOPOLE % Creates a single monopole or a monopole array % on a finite ground plane using mouse input % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 4 clear all warning off % The following parameters need to be specified: L=2.0; %Plate length (along the x-axis) W=2.0; %Plate width (along the y-axis) Nx=11; %Discretization parameter (length) Ny=11; %Discretization parameter (width) h=1.0; %Monopole height Number=7; %Number of monopole rectangles %Set vertexes of the plate epsilon=1e-6; M=1; for i=1:Nx+1 for j=1:Ny+1 X(M)=-L/2+(i-1)/Nx*L; Y(M)=-W/2+(j-1)/Ny*W-epsilon*X(M); M=M+1; end end %Identify probe feed edge(s) x=[-0.02 0.02]; y=[0 0]; X=[X x]; Y=[Y y]; C=mean(x); x1=[C C]; y1=mean(y)+2*[max(x)-C min(x)-C]; X=[X x1]; Y=[Y y1]; %Use Delaunay triangulation TRI = delaunay(X,Y,0); t=TRI'; t(4,:)=0; %plate has the index zero!!! p=[X; Y; zeros(1,length(X))]; save plate p t %Identify junction edge(s)using mouse input %Use RETURN key to stop the process clear figure hold off viewer plate; view(0,90); hold on FeedingTriangle=[]; TRI=t(1:3,:)'; while ~isempty(t) [xi,yi]=ginput(1); TriangleNumber = tsearch(p(1,:),p(2,:),TRI,xi,yi); n=t(1:3,TriangleNumber); FeedingTriangle= [FeedingTriangle TriangleNumber]; x= p(1,n); y= p(2,n); if isempty(xi|yi) break; end fill(x,y,'w') clear xi yi end clear figure %Create the monopole(s) for n=1:length(FeedingTriangle)/2 %find bottom (feeding) edge FT=[FeedingTriangle(2*n-1) FeedingTriangle(2*n)]; N=t(1:3,FT(1)); M=t(1:3,FT(2)); a=1-all([N-M(1) N-M(2) N-M(3)]); Edge_B=M(find(a)); %find top edge p=[p p(:,Edge_B(1))+[0;0;h] p(:,Edge_B(2))+[0;0;h]]; Edge_T =[length(p)-1; length(p)]; %fill the strip Edge_MM=Edge_B; for k=1:Number-1 p(:,length(p)+1)=k/Number*(p(:,Edge_T(1))-p(:,Edge_B(1)))+p(:,Edge_B(1)); p(:,length(p)+1)=k/Number*(p(:,Edge_T(2))-p(:,Edge_B(2)))+p(:,Edge_B(2)); Edge_M=[length(p)-1,length(p)]; tFeed1(:,k) =[Edge_MM(1);Edge_MM(2);Edge_M(2);1]; tFeed2(:,k) =[Edge_MM(1);Edge_M(1);Edge_M(2);1]; Edge_MM=Edge_M; end %update array t tFeed3 =[Edge_M(1);Edge_M(2);Edge_T(2);1]; tFeed4 =[Edge_M(1);Edge_T(1);Edge_T(2);1]; t=[t tFeed1 tFeed2 tFeed3 tFeed4]; end %Save result save monopole p t hold off clear figure viewer monopole Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/multi1.m%Problem 4.5 to Chapter 4 % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 4 clear all load platefine p=p*2; d=0.5; p(3,:)=-d; pbase=p; tbase=t; clear p, t; load strip2 x=p(1,:); y=p(2,:); z=p(3,:); p(1,:)=p(1,:); p(2,:)=p(2,:); p(3,:)=0; t=[tbase t+length(pbase)]; t(4,:)=1; p=[pbase p]; save dip p t viewer('dip') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/multi2.m%Problem 4.7 to Chapter 4 % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 4 clear all load platefine p=p*2; R=1.25; d=1; p(3,:)=R-d-sqrt(R^2-p(1,:).^2); pbase=p; tbase=t; clear p, t; load strip2 x=p(1,:); y=p(2,:); z=p(3,:); p(1,:)=p(1,:); p(2,:)=p(2,:); p(3,:)=0; t=[tbase t+length(pbase)]; t(4,:)=1; p=[pbase p]; save dip p t viewer('dip') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/multi3.m%Problem 4.9 to Chapter 4 % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 4 clear all load reflector d=2; p(3,:)=p(3,:)-d; pbase=p; tbase=t; clear p, t; load strip2 x=p(1,:); y=p(2,:); z=p(3,:); p(1,:)=p(1,:); p(2,:)=p(2,:); p(3,:)=0; t=[tbase t+length(pbase)]; t(4,:)=1; p=[pbase p]; save dip p t viewer('dip') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/multi4.m%Problem 4.11 to Chapter 4 % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 4 clear all load mesh %p=p*2; d=0.5; p(3,:)=-d; pbase=p; tbase=t; clear p, t; load strip2 x=p(1,:); y=p(2,:); z=p(3,:); p(1,:)=p(1,:); p(2,:)=p(2,:); p(3,:)=0; t=[tbase t+length(pbase)]; t(4,:)=1; p=[pbase p]; save dip p t viewer('dip') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/platefine.matp:[3x289 double array] t:[4x512 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/reflector.matp:[3x154 double array] t:[4x274 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/strip.m%STRIP % Creates a dipole antenna (thin strip) in the xy-plane % % The following parameters need to be specified: % Strip width (along the x-axis) W % Strip length (along the y-axis) L % Discretization parameter (width) Nx % Discretization parameter (length) Ny % % Note: the equivalent wire radius is 0.25*W % % For more general version of this script please % see strip_.m of Appendix A. % % Copyright 2002 AEMM. Revision 2002/04/12 % Chapter 4 clear all warning off W=0.02; L=2.00; Nx=1; Ny=10; %Set vertexes of the strip epsilon=1e-6; M=1; for i=1:Nx+1 for j=1:Ny+1 X(M)=-W/2+(i-1)/Nx*W; Y(M)=-L/2+(j-1)/Ny*L-epsilon*X(M); M=M+1; end end %Use Delaunay triangulation TRI = delaunay(X,Y,0); t=TRI'; t(4,:)=1; p=[X; Y; zeros(1,length(X))]; %Save mesh file save strip p t; viewer strip Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/strip2.matp:[3x234 double array] t:[4x244 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/mesh/viewer.mfunction []=viewer(filename) %VIEWER Visualizes the structure - all Chapters % % Usage: viewer('plate.mat') or % viewer('plate') or % viewer plate % % Copyright 2002 AEMM. Revision 2002/03/05 Chapter 2 %Load the data load(filename); %Color=[0.75 0.75 0.75]; Color='g'; Cut=find (t(4,:)>1); t(:,Cut)=[]; for m=1:length(t) n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end g=fill3(X,Y,Z,Color); axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/point.mfunction[EField, HField]=... point1(Point,eta_,K,DipoleMoment,DipoleCenter) %POINT Radiated/scattered field at a point of a dipole array % or a single dipole. Gives exact near- and far-fields. Outputs % individual contribution of each dipole. % % Observation point Point(1:3) % Array of dipole moments DipoleMoment(1:3,1:EdgesTotal) % Array of dipole centers DipoleCenter(1:3,1:EdgesTotal) % E-field at the observation point E(1;3,1:EdgesTotal) % H-field at the observation point H(1;3,1:EdgesTotal) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 C=4*pi; ConstantH=K/C; ConstantE=eta_/C; m=DipoleMoment; c=DipoleCenter; r =repmat(Point,[1 length(c)])-c(1:3,:); PointRM =repmat(sqrt(sum(r.*r)),[3 1]); EXP =exp(-K*PointRM); PointRM2=PointRM.^2; C=1./PointRM2.*(1+1./(K*PointRM)); D=repmat(sum(r.*m),[3 1])./PointRM2; M=D.*r; HField=ConstantH*cross(m,r).*C.*EXP; EField=ConstantE*((M-m).*(K./PointRM+C)+2*M.*C).*EXP; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/impmet.mfunction [Z]= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %IMPMET Standard impedance matrix (metal surface) % % Returns the complex impedance matrix [EdgesTotal x EdgesTotal] % Uses 9 integration points for every triangle % (barycentric subdivision) % % The impedance matrix is calculated as a sum of the contributions % due to separate triangles (similar to the "face-pair" method). % See Appendix B for a detailed algorithm. % % A 9-point quadrature is used for all integrals, including % the self-coupling terms. The alternative source code with % the analytical approximation of the self-coupling terms % is given in Appendix B. The difference between two methods % is not significant. % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 2 %Memory allocation Z =zeros (EdgesTotal,EdgesTotal)+j*zeros(EdgesTotal,EdgesTotal); %Loop over integration triangles for p=1:TrianglesTotal Plus =find(TrianglePlus-p==0); Minus =find(TriangleMinus-p==0); D=Center_-repmat(Center(:,p),[1 9 TrianglesTotal]); %[3 9 TrianglesTotal] R=sqrt(sum(D.*D)); %[1 9 TrianglesTotal] g=exp(-K*R)./R; %[1 9 TrianglesTotal] gP=g(:,:,TrianglePlus); %[1 9 EdgesTotal] gM=g(:,:,TriangleMinus); %[1 9 EdgesTotal] Fi=sum(gP)-sum(gM); %[1 1 EdgesTotal] ZF= FactorFi.*reshape(Fi,EdgesTotal,1); %[EdgesTotal 1] for k=1:length(Plus) n=Plus(k); RP=repmat(RHO__Plus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1+ZF); end for k=1:length(Minus) n=Minus(k); RP=repmat(RHO__Minus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1-ZF); end end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/mesh/strip2.matp:[3x234 double array] t:[4x244 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/rwg1.m%RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % Copyright 2002 AEMM. Revision 2002/03/09 Chapter 2 clear all tic; load('mesh/strip2'); [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) %All structures of this chapter have EdgeIndicator=2 EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =75e6; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/rwg4.m%RWG4 Solves MoM equations for the scattering problem % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % Direction of the incident signal in Cartesian coordinates % d(1:3); % Direction of the E-field in the incident plane wave % in Cartesian coordinates Pol(1:3); % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('impedance'); %Incident field %Example: d=[0 0 -1] means that the incident signal % is in the -z direction. %%Custom incidence d =[0 0 1]; Pol =[0 -0.0044-j*0.0049 0]; k=omega/c_; kv=k*d; for m=1:EdgesTotal ScalarProduct=sum(kv.*Center(:,TrianglePlus(m))'); EmPlus =Pol.'*exp(-j*ScalarProduct); ScalarProduct=sum(kv.*Center(:,TriangleMinus(m))'); EmMinus=Pol.'*exp(-j*ScalarProduct); ScalarPlus =sum(EmPlus.* RHO_Plus(:,m)); ScalarMinus=sum(EmMinus.*RHO_Minus(:,m)); V(m)=EdgeLength(m)*(ScalarPlus/2+ScalarMinus/2); end tic; %System solution I=Z\V.'; toc %elapsed time FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','I','V','d','Pol'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/receivingantenna/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index=find(t(4,:)1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =75e6; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/rwg4.m%RWG4 Solves MoM equations for the antenna radiation problem % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % The feed point position FeedPoint(1:3); % Number of feeding edges (one for the dipole; % two for the monopole) INDEX(1:2); % % Copyright 2002 AEMM. Revision 2002/03/14 % Chapter 4 %load the data load('mesh2'); load('impedance'); %Find the feeding edge(s)(closest to the origin) FeedPoint=[0; 0; 0]; for m=1:EdgesTotal V(m)=0; Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-FeedPoint; end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=INDEX(1); %Center feed - dipole %Index=INDEX(1:2); %Probe feed - monopole %Define the voltage vector V(Index)=1*EdgeLength(Index); %Solve system of MoM equations tic; I=Z\V.'; toc %elapsed time %Find the antenna input impedance GapCurrent =sum(I(Index).*EdgeLength(Index)'); GapVoltage =mean(V(Index)./EdgeLength(Index)); Impedance =GapVoltage/GapCurrent FeedPower =1/2*real(GapCurrent*conj(GapVoltage)) FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_',... 'I','V','GapCurrent','GapVoltage','Impedance','FeedPower'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter04/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index=find(t(4,:)1) t=[t [Count-2; Count-1; Count+1] [Count-2; Count; Count+1]]; end Count=Count+2; end %Close the chain t=[t [Count-2; Count-1; 2] [Count-2; 1; 2]]; %Nodes PointNumber=Count-1; for L=1:PointNumber p(1:3,L) = [X(L); Y(L); Z(L)]; end t(4,:)=1; %Save result save loop1 p t; viewer('loop1') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/loop1.matt:[4x180 double array] p:[3x180 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/loop2.m%LOOP2 % Creates triangular mesh for the helical antenna % of given radius, number of turns, spacing, and % wire thickness. % % The following parameters need to be specified: % % Turn radius in m a % Number of loop rectangles M % Width of the strip h % Number of turns N % Spacing between turns S % % Note: the equivalent wire radius is 0.25*h % % Copyright 2002 AEMM. Revision 2002/03/13 % Chapter 5 clear all %Normal mode parameters a=0.1; M=40; h=0.005; N=9; S=0.04; %Pitch angle pitch=atan(S/(2*pi*a)) factor=sin(pitch); factor1=h*cos(pitch); L=N*S %Total length of the antenna Count=1; %Point number %Create rectangles t=[]; for n=1:M*N angle=2*pi*(n-1)/M; x= a*cos(angle); y= a*sin(angle); zM= n*L/(M*N)-L/2+h/2; X(Count:Count+1)=[x x]'; Y(Count:Count+1)=[y y]'; Z(Count:Count+1)=[zM zM-h]'; if(n>1) t=[t [Count-2; Count-1; Count+1] [Count-2; Count; Count+1]]; end Count=Count+2; end %Nodes PointNumber=Count-1; for L=1:PointNumber p(1:3,L) = [X(L); Y(L); Z(L)]; end t(4,:)=1; %Save result save loop2 p t h; viewer('loop2') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/loop3.m%LOOP3 % Creates triangular mesh for the helical antenna % of given radius, number of turns, spacing, and % wire thickness. % % LOOP3 has a somewhat better mesh quality than LOOP2 % (see the text of Chapter 5). Otherwise, both the scripts % are equivalent. % % The following parameters need to be specified: % % Turn radius in m a % Number of loop rectangles M % Width of the strip h % Number of turns N % Spacing between turns S % % Note: the equivalent wire radius is 0.25*h % % Copyright 2002 AEMM. Revision 2002/03/13 % Chapter 5 clear all %Axial mode parameters a=0.0545; M=12; h=0.005; N=15; S=0.076; %pitch angle pitch=atan(S/(2*pi*a)) factor=sin(pitch); factor1=h*cos(pitch); L=N*S %Total length of the antenna Count=1; %Point number EN=1; %Element number %Create rectangles t=[]; for n=1:M*N angle =2*pi*(n-1)/M; delta =h*factor/a; %correction angle in radians Point1=[a*cos(angle) a*sin(angle) n*L/(M*N)-L/2]; Point2=[a*cos(angle+delta) a*sin(angle+delta) n*L/(M*N)-L/2-factor1]; X(Count:Count+1)=[Point1(1) Point2(1)]'; Y(Count:Count+1)=[Point1(2) Point2(2)]'; Z(Count:Count+1)=[Point1(3) Point2(3)]'; if(n>1) t=[t [Count-2; Count-1; Count+1] [Count-2; Count; Count+1]]; end Count=Count+2; end %Nodes PointNumber=Count-1; for L=1:PointNumber p(1:3,L) = [X(L); Y(L); Z(L)]; end t(4,:)=1; %Save result save loop3 p t h; viewer('loop3') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/loop4.m%LOOP4 % Creates triangular mesh for the helical tapered % antenna of given radius, number of turns, spacing, % and wire thickness. % % The following parameters need to be specified: % % Turn radius in m -Center amin % Turn radius in m -Top/Bottom amax % Number of loop rectangles M % Width of the strip h % Number of turns N % Spacing between turns S % % Note: the equivalent wire radius is 0.25*h % % Copyright 2002 AEMM. Revision 2002/03/13 % Chapter 5 clear all amin=0.1; amax=1; M=30; h=0.05; N=10; S=0.2; L=N*S; %Total length of the antenna Count=1; %Point number %Create rectangles t=[]; for n=1:M*N angle=2*pi*(n-1)/M; R=amin-abs((n-1)-M*N/2)*(amin-amax)/(M*N/2); x= R*cos(angle); y= R*sin(angle); zM= n*L/(M*N)-L/2+h/2; X1Array(Count:Count+1)=[x x]'; X2Array(Count:Count+1)=[y y]'; X3Array(Count:Count+1)=[zM zM-h]'; if(n>1) t=[t [Count-2; Count-1; Count+1] [Count-2; Count; Count+1]]; end Count=Count+2; end %Nodes PointNumber=Count-1; for L=1:PointNumber p(1:3,L) = [X1Array(L); X2Array(L); X3Array(L)]; end t(4,:)=1; %Save result save loop4 p t h; viewer('loop4') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/strip.matt:[4x180 double array] p:[3x182 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/mesh/viewer.mfunction []=viewer(filename) %VIEWER Visualizes the structure - all Chapters % % Usage: viewer('plate.mat') or % viewer('plate') or % viewer plate % % Copyright 2002 AEMM. Revision 2002/03/05 Chapter 2 %Load the data load(filename); %Color=[0.75 0.75 0.75]; Color='g'; Cut=find (t(4,:)>1); t(:,Cut)=[]; for m=1:length(t) n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end g=fill3(X,Y,Z,Color); axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/point.mfunction[EField, HField]=... point1(Point,eta_,K,DipoleMoment,DipoleCenter) %POINT Radiated/scattered field at a point of a dipole array % or a single dipole. Gives exact near- and far-fields. Outputs % individual contribution of each dipole. % % Observation point Point(1:3) % Array of dipole moments DipoleMoment(1:3,1:EdgesTotal) % Array of dipole centers DipoleCenter(1:3,1:EdgesTotal) % E-field at the observation point E(1;3,1:EdgesTotal) % H-field at the observation point H(1;3,1:EdgesTotal) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 C=4*pi; ConstantH=K/C; ConstantE=eta_/C; m=DipoleMoment; c=DipoleCenter; r =repmat(Point,[1 length(c)])-c(1:3,:); PointRM =repmat(sqrt(sum(r.*r)),[3 1]); EXP =exp(-K*PointRM); PointRM2=PointRM.^2; C=1./PointRM2.*(1+1./(K*PointRM)); D=repmat(sum(r.*m),[3 1])./PointRM2; M=D.*r; HField=ConstantH*cross(m,r).*C.*EXP; EField=ConstantE*((M-m).*(K./PointRM+C)+2*M.*C).*EXP; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/rwg1.m%RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % Copyright 2002 AEMM. Revision 2002/03/09 Chapter 2 clear all tic; load('mesh/loop1'); [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) %All structures of this chapter have EdgeIndicator=2 EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) %f=4.8e6 means that C=lambda/10 (loop with a =1 m) f =4.8e6; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/rwg4.m%RWG4 Solves MoM equations for the antenna radiation problem % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % The feed point position FeedPoint(1:3); % Number of feeding edges (one for the dipole; % two for the monopole) INDEX(1:2); % % Copyright 2002 AEMM. Revision 2002/03/14 % Chapter 4 %load the data load('mesh2'); load('impedance'); %Find the feeding edge(s)(closest to the origin) FeedPoint=[-1; 0; 0]; for m=1:EdgesTotal V(m)=0; Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-FeedPoint; end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=INDEX(1); %Center feed - dipole and loop %Define the voltage vector V(Index)=1*EdgeLength(Index); %Solve system of MoM equations tic; I=Z\V.'; toc %elapsed time %Find the antenna input impedance GapCurrent =sum(I(Index).*EdgeLength(Index)'); GapVoltage =mean(V(Index)./EdgeLength(Index)); Impedance =GapVoltage/GapCurrent FeedPower =1/2*real(GapCurrent*conj(GapVoltage)) FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_',... 'I','V','GapCurrent','GapVoltage','Impedance','FeedPower','Index'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter05/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index_=find(t(4,:)1); t(:,Cut)=[]; for m=1:length(t) n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end g=fill3(X,Y,Z,Color); axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/efield3.mfunction [PolarXY, GAIN, phi]=efield3; %EFIELD3 2D Radiation patterns -function % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Additionally uses the value of TotalPower saved % in file gainpower.mat (script efield2.m) % % The following parameters need to be specified: % % Radius of the circle (m) R % Plane of the circle: [x y 0] or % [x 0 z] or % [0 y z] % Number of discretization points per % pattern NumPoints % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 6/Loop 2 clear all %Load the data load('mesh2'); load('current'); k=omega/c_; K=j*k; for m=1:EdgesTotal Point1=Center(:,TrianglePlus(m)); Point2=Center(:,TriangleMinus(m)); DipoleCenter(:,m)=0.5*(Point1+Point2); DipoleMoment(:,m)=EdgeLength(m)*I(m)*(-Point1+Point2); end NumPoints=100; R=1000; %pattern in m for ii=1:NumPoints+1 phi(ii)=(ii-1)*pi/(NumPoints/2); x=R*cos(phi(ii)); y=R*sin(phi(ii)); ObservationPoint=[x y 0]'; [E,H]=point(ObservationPoint,eta_,K,DipoleMoment,DipoleCenter); ET=sum(E,2); HT=sum(H,2); Poynting=0.5*real(cross(ET,conj(HT))); W(ii)=norm(Poynting); U(ii)=(norm(ObservationPoint))^2*W(ii); end PolarXY=10*log10(4*pi*U/sum(FeedPower)); GAIN=max(PolarXY); %gain for the particular pattern! Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/impmet.mfunction [Z]= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %IMPMET Standard impedance matrix (metal surface) % % Returns the complex impedance matrix [EdgesTotal x EdgesTotal] % Uses 9 integration points for every triangle % (barycentric subdivision) % % The impedance matrix is calculated as a sum of the contributions % due to separate triangles (similar to the "face-pair" method). % See Appendix B for a detailed algorithm. % % A 9-point quadrature is used for all integrals, including % the self-coupling terms. The alternative source code with % the analytical approximation of the self-coupling terms % is given in Appendix B. The difference between two methods % is not significant. % % Copyright 2002 AEMM. Revision 2002/03/12 % Chapter 2 %Memory allocation Z =zeros (EdgesTotal,EdgesTotal)+j*zeros(EdgesTotal,EdgesTotal); %Loop over integration triangles for p=1:TrianglesTotal Plus =find(TrianglePlus-p==0); Minus =find(TriangleMinus-p==0); D=Center_-repmat(Center(:,p),[1 9 TrianglesTotal]); %[3 9 EdgesTotal] R=sqrt(sum(D.*D)); %[1 9 TrianglesTotal] g=exp(-K*R)./R; %[1 9 TrianglesTotal] gP=g(:,:,TrianglePlus); %[1 9 EdgesTotal] gM=g(:,:,TriangleMinus); %[1 9 EdgesTotal] Fi=sum(gP)-sum(gM); %[1 1 EdgesTotal] ZF= FactorFi.*reshape(Fi,EdgesTotal,1); for k=1:length(Plus) n=Plus(k); RP=repmat(RHO__Plus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1+ZF); end for k=1:length(Minus) n=Minus(k); RP=repmat(RHO__Minus(:,:,n),[1 1 EdgesTotal]); %[3 9 EdgesTotal] A=sum(gP.*sum(RP.*RHO_P))+sum(gM.*sum(RP.*RHO_M)); Z1= FactorA.*reshape(A,EdgesTotal,1); Z(:,n)=Z(:,n)+EdgeLength(n)*(Z1-ZF); end end Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/loop.m%LOOP2 Loop versus a phase delay for the linear array % of dipoles. % Outputs the radiation patterns in the azimuthal plane. % The array axis is the x-axis (cf. Fig. 6.10 of Chapter 6) % The dipole axis is the z-axis (cf. Fig. 6.10 of Chapter 6) % % The following parameters need to be specified: % % Number of sampling points for phase sweep K % Spacing between array elements (m) d % Number of dipoles in the array N % % Copyright 2002 AEMM. Revision 2002/03/16 % Chapter 6 clear all K=20; Step=pi/K; %rad d=1.0; N=2; multilinear(d,N); rwg1; rwg2; rwg3; for k=1:K+1 k phase(k) =-(k-1)*Step; Power1 =rwg4(phase(k)); [PolarXY, GAIN, phi] =efield3; PolarXYLoop (k,:) =PolarXY; GainLoop (k) =GAIN; PowerLoop (k) =Power1; Radius =max(PolarXY,0); x1(k,:)=Radius.*cos(phi); y1(k,:)=Radius.*sin(phi); z1(k,:)=-phase(k)*ones(length(x1),1)'; end surf(x1,y1,z1); colormap gray; rotate3d; xlabel('gain, dB'); ylabel('gain, dB'); zlabel('-phase, rad'); view(-12,20) save loop x1 y1 z1 K phi PolarXYLoop GainLoop PowerLoop phase Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/multilinear.mfunction []=multilinear(d,N); %MULTILINEAR - function (Chapter 6/Loop 2) % Creates a linear array of N dipoles separated by d % % Copyright 2002 AEMM. Revision 2002/03/06 % Chapter 6 %Load the dipole structure load strip1 p(3,:)=p(2,:); p(2,:)=0; pbase=p; tbase=t; Length=(N-1)*d; if(N==1) save array p t; break; end Step=d; p(1,:)=p(1,:)-Length/2; %to position the array Feed(1:3,1)=[mean(p(1,:)),0,0]'; %Cloning algorithm for k=1:N-1 pbase1(1,:)=pbase(1,:)-Length/2+k*Step; pbase1(2,:)=pbase(2,:); pbase1(3,:)=pbase(3,:); Feed(1:3,k+1)=[mean(pbase1(1,:)),0,0]'; p=[p pbase1]; t=[t tbase+k*length(pbase1)]; t(4,:)=1; end save array p t Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/point.mfunction[EField, HField]=... point1(Point,eta_,K,DipoleMoment,DipoleCenter) %POINT Radiated/scattered field at a point of a dipole array % or a single dipole. Gives exact near- and far-fields % % Observation point Point(1:3) % Array of dipole moments DipoleMoment(1:3,1:EdgesTotal) % Array of dipole centroids DipoleCenter(1:3,1:EdgesTotal) % E-field at the observation point E(1;3,1:EdgesTotal) % H-field at the observation point H(1;3,1:EdgesTotal) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 C=4*pi; ConstantH=K/C; ConstantE=eta_/C; m=DipoleMoment; c=DipoleCenter; r =repmat(Point,[1 length(c)])-c(1:3,:); PointRM =repmat(sqrt(sum(r.*r)),[3 1]); EXP =exp(-K*PointRM); PointRM2=PointRM.^2; C=1./PointRM2.*(1+1./(K*PointRM)); D=repmat(sum(r.*m),[3 1])./PointRM2; M=D.*r; HField=ConstantH*cross(m,r).*C.*EXP; EField=ConstantE*((M-m).*(K./PointRM+C)+2*M.*C).*EXP; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/rwg1.mfunction []=rwg1 %RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 tic; load('array'); [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/rwg2.mfunction []=rwg2 %RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/rwg3.mfunction []=rwg3 %RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =75e6; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/rwg4.mfunction [PowerSum]=rwg4(phase); %RWG4 Solves MoM equations for the antenna radiation problem - % antenna arrays % % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % The feed point positions Feed(1:3,1:N); % Number of feeding edges (one for the dipole; % two for the monopole) INDEX(1:2); % % Copyright 2002 AEMM. Revision 2002/03/16 % Chapter 6 %Load the data load('mesh2'); load('impedance'); V(1:EdgesTotal)=0; %Find feeding edges closest to the array Feed N=length(Feed(1,:)); Index=[]; for k=1:N for m=1:EdgesTotal Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-Feed(:,k); end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=[Index INDEX(1)]; %Center feed - dipole %Index=[Index INDEX(1:2)]; %Probe feed - monopole end V(Index)=1.0*EdgeLength(Index); N=length(Index); %Identify feeding voltages-linear array (dipole array only!) for n=1:N nn=Index(n); V(nn)=V(nn)*exp(j*phase*(n-1)); end %Solve system of MoM equations tic; I=Z\V.'; toc %elapsed time %Terminal impedance (dipole array only!) for n=1:N nn=Index(n); GapCurrent(n)=I(nn)*EdgeLength(nn); GapVoltage(n)=V(nn)/EdgeLength(nn); Impedance (n)=GapVoltage(n)/GapCurrent(n); %this is the terminal impedance FeedPower (n)=1/2*real(GapCurrent(n)*conj(GapVoltage(n))); end PowerSum=sum(FeedPower); FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_',... 'I','V','GapCurrent','GapVoltage','Impedance','FeedPower'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/loop2/strip1.matt:[4x80 double array] p:[3x82 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/array.matt:[4x160 double array] p:[3x164 double array] Feed:[3x2 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/bowtie.matt:[4x30 double array] p:[3x26 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/multibowtie.m%MULTIBOWTIE % Creates a planar array of NxN bowties (or other antenna elements) % over a finite ground plane % % The following parameters need to be specified: % % Antenna element filename filename % Number of elements, N, in the NxN array N % Maximum size of the element (m) SizeE % Size of the square ground plane (m) SizeP % Height of elements above ground plane (m) Height % % Note that the following domain numbers are assigned: % Triangles of the metal ground plane t(4,:)=0 % Triangles of the antenna elements t(4,:)=1 % % Copyright 2002 AEMM. Revision 2002/03/06 % Chapter 6 clear all filename='bowtie'; N=8; SizeE=0.5; SizeP=6; Height=0.5; load(filename); t(4,:)=1; tbase=t; p=p/max(max(p))*SizeE/2; %Set centers of array elements h=SizeP/N; Count=1; for i=1:N x(i)=h/2+(i-1)*h-SizeP/2; for j=1:N y(j)=h/2+(j-1)*h-SizeP/2; Feed(1:2,Count)=[x(i),y(j)]'; Feed(3:3,Count)=Height; Indicator(Count)=(-1)^i; Count=Count+1; end end %Clone array elements P=[]; T=[]; X=p(1,:); Y=p(2,:); for k=1:Count-1 if(Indicator(k)==1) pbase(1,:)=X+Feed(1,k); pbase(2,:)=Y+Feed(2,k); else pbase(1,:)=Y+Feed(1,k); pbase(2,:)=X+Feed(2,k); end pbase(3,:)=p(3,:)+Feed(3,k); P=[P pbase]; T=[T tbase+(k-1)*length(pbase)]; T(4,:)=1; end %Attach the finite ground plane clear p t; L=SizeP; %Plate length (along the x-axis) W=SizeP; %Plate width (along the y-axis) Nx=16; %Discretization parameter (length) Ny=16; %Discretization parameter (width) plate(L,W,Nx,Ny); load plate; t_=t+length(P); t_(4,:)=0; p=[P p]; t=[T t_]; TrianglesTotal=length(t) save array p t Feed viewer('array') grid on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/multicircular.m%MULTICIRCLULAR % Creates a circular array of N dipoles % equally spaced over the circle of radius R % % The following parameters need to be specified: % % Number of elements, N, in the array N % Radius of the circular array (m) R % % Copyright 2002 AEMM. Revision 2002/03/06 % Chapter 6 clear all %Load the dipole structure load strip1 N=16; %Number of array elements R=3; %Radius of circular array p(3,:)=p(2,:); p(2,:)=0; pbase=p; tbase=t; Theta =2*pi/N; p(1,:)=p(1,:)-R; Feed(1:3,1)=[mean(p(1,:)),mean(p(2,:)),0]'; %Cloning algorithm for k=1:N-1 pbase1(1,:)=pbase(1,:)+R*cos(-pi+k*Theta); pbase1(2,:)=pbase(2,:)+R*sin(-pi+k*Theta); pbase1(3,:)=pbase(3,:); Feed(1:3,k+1)=[mean(pbase1(1,:)),mean(pbase1(2,:)),0]'; p=[p pbase1]; t=[t tbase+k*length(pbase1)]; t(4,:)=1; end save array p t Feed viewer('array') grid on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/multilinear.m%MULTILINEAR % Creates a linear array of N dipoles separated by % distance d % % The following parameters need to be specified: % % Number of elements, N, in the array N % Separation distance (m) d % % Copyright 2002 AEMM. Revision 2002/03/06 % Chapter 6 clear all N=2; %Number of array elements d=2.000; %Separation distance between elements %Load the dipole structure load strip1 p(3,:)=p(2,:); p(2,:)=0; pbase=p; tbase=t; Length=(N-1)*d; if(N==1) save array p t; break; end Step=d; p(1,:)=p(1,:)-Length/2; %to position the array Feed(1:3,1)=[mean(p(1,:)),0,0]'; %Cloning algorithm for k=1:N-1 pbase1(1,:)=pbase(1,:)-Length/2+k*Step; pbase1(2,:)=pbase(2,:); pbase1(3,:)=pbase(3,:); Feed(1:3,k+1)=[mean(pbase1(1,:)),0,0]'; p=[p pbase1]; t=[t tbase+k*length(pbase1)]; t(4,:)=1; end save array p t Feed viewer('array') grid on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/multimonopole.m%MULTIMONOPOLE % Creates an array of N base-fed % monopoles on a finite ground plane % % Uses mouse input to identify the junction edges % % Note that the following domain numbers are assigned: % Triangles of the metal ground plane t(4,:)=0 % Triangles of the monopoles t(4,:)=1 % % Copyright 2002 AEMM. Revision 2002/03/06 % Chapter 4 clear all warning off % The following parameters need to be specified: L=2.0; %Plate length (along the x-axis) W=2.0; %Plate width (along the y-axis) Nx=16; %Discretization parameter (length) Ny=16; %Discretization parameter (width) h=1.0; %Monopole height Number=7; %Number of monopole rectangles plate(L,W,Nx,Ny); viewer plate; view(0,90); load plate; hold on FeedingTriangle=[]; TRI=t(1:3,:)'; %Mouse input while ~isempty(t) [xi,yi]=ginput(1); TriangleNumber = tsearch(p(1,:),p(2,:),TRI,xi,yi); n=t(1:3,TriangleNumber); FeedingTriangle= [FeedingTriangle TriangleNumber]; x= p(1,n); y= p(2,n); if isempty(xi|yi) break; end fill(x,y,'w') clear xi yi end clear figure %Create the monopole(s) for n=1:length(FeedingTriangle)/2 %find feeding edge FT=[FeedingTriangle(2*n-1) FeedingTriangle(2*n)]; N=t(1:3,FT(1)); M=t(1:3,FT(2)); a=1-all([N-M(1) N-M(2) N-M(3)]); Edge_B=M(find(a)); Feed(:,n)=0.5*(p(:,Edge_B(1))+p(:,Edge_B(2))); %set top edge p=[p p(:,Edge_B(1))+[0;0;h] p(:,Edge_B(2))+[0;0;h]]; Edge_T =[length(p)-1; length(p)]; %fill the strip Edge_MM=Edge_B; for k=1:Number-1 p(:,length(p)+1)=k/Number*(p(:,Edge_T(1))-p(:,Edge_B(1)))+p(:,Edge_B(1)); p(:,length(p)+1)=k/Number*(p(:,Edge_T(2))-p(:,Edge_B(2)))+p(:,Edge_B(2)); Edge_M=[length(p)-1,length(p)]; tFeed1(:,k) =[Edge_MM(1);Edge_MM(2);Edge_M(2);1]; tFeed2(:,k) =[Edge_MM(1);Edge_M(1);Edge_M(2);1]; Edge_MM=Edge_M; end %update array t tFeed3 =[Edge_M(1);Edge_M(2);Edge_T(2);1]; tFeed4 =[Edge_M(1);Edge_T(1);Edge_T(2);1]; t=[t tFeed1 tFeed2 tFeed3 tFeed4]; end save array p t Feed hold off clear figure viewer array grid on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/plate.mfunction []=plate(W,L,Nx,Ny); %PLATE Plate mesh in the xy-plane -function % % W Plate width (along the x-axis) % L Plate length (along the y-axis) % Nx Discretization parameter (width) % Ny Discretization parameter (length) % % Copyright 2002 AEMM. Revision 2002/03/16 % Chapter 6; see also the equivalent script plate.m % of Chapter 2 and Appendix A. if nargin < 4 error('function plate requires four input arguments.') end warning off %Set vertexes epsilon=1e-6; M=1; for i=1:Nx+1 for j=1:Ny+1 X(M)=-W/2+(i-1)/Nx*W; Y(M)=-L/2+(j-1)/Ny*L-epsilon*X(M); M=M+1; end end %Use Delaunay triangulation TRI = delaunay(X,Y,0); t=TRI'; t(4,:)=0; p=[X; Y; zeros(1,length(X))]; %Save result save plate p t; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/strip1.matt:[4x80 double array] p:[3x82 double array] Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/mesh/viewer.mfunction []=viewer(filename) %VIEWER Visualizes the structure % Shows the ground plane and array elements % using different colors % % Usage: viewer('plate.mat') or % viewer('plate') or % viewer plate % % Copyright 2002 AEMM. Revision 2002/03/07 % Chapter 6 load(filename); Color0='g'; %[0.45 0.45 0.45]; Color1='y'; %[0.75 0.75 0.75]; TrianglesTotal=length(t) GroundPlane =find(t(4,:)==0); AntennaElements =find(t(4,:)==1); for k=1:length(GroundPlane) m=GroundPlane(k); n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end for k=1:length(AntennaElements) m=AntennaElements(k); n=t(1:3,m); X_(1:3,m)=p(1,n)'; Y_(1:3,m)=p(2,n)'; Z_(1:3,m)=p(3,n)'; end if(~isempty(GroundPlane)) g=fill3(X,Y,Z,Color0,X_,Y_,Z_,Color1); else g=fill3(X_,Y_,Z_,Color1); end axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/point.mfunction[EField, HField]=... point1(Point,eta_,K,DipoleMoment,DipoleCenter) %POINT Radiated/scattered field at a point of a dipole array % or a single dipole. Gives exact near- and far-fields. Outputs % individual contribution of each dipole. % % Observation point Point(1:3) % Array of dipole moments DipoleMoment(1:3,1:EdgesTotal) % Array of dipole centers DipoleCenter(1:3,1:EdgesTotal) % E-field at the observation point E(1;3,1:EdgesTotal) % H-field at the observation point H(1;3,1:EdgesTotal) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 C=4*pi; ConstantH=K/C; ConstantE=eta_/C; m=DipoleMoment; c=DipoleCenter; r =repmat(Point,[1 length(c)])-c(1:3,:); PointRM =repmat(sqrt(sum(r.*r)),[3 1]); EXP =exp(-K*PointRM); PointRM2=PointRM.^2; C=1./PointRM2.*(1+1./(K*PointRM)); D=repmat(sum(r.*m),[3 1])./PointRM2; M=D.*r; HField=ConstantH*cross(m,r).*C.*EXP; EField=ConstantE*((M-m).*(K./PointRM+C)+2*M.*C).*EXP; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/rwg1.m%RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all tic; load('mesh/array'); [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 2 clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =75e6; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/rwg4.m%RWG4 Solves MoM equations for the antenna radiation problem - % antenna arrays % % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % The feed point positions for N array elements % Feed(1:3,1:N); % Number of feeding edges (one for the dipole; % two for the monopole) INDEX(1:2); % % Copyright 2002 AEMM. Revision 2002/03/16 % Chapter 6 clear all %Load the data load('mesh2'); load('impedance'); V(1:EdgesTotal)=0; %Find feeding edges closest to the array Feed N=length(Feed(1,:)); Index=[]; for k=1:N for m=1:EdgesTotal Distance(m)=norm(0.5*sum(p(:,Edge_(:,m)),2)-Feed(:,k)); end [Y,INDEX]=sort(Distance); Index=[Index INDEX(1)]; %Center feed - dipole %Index=[Index INDEX(1:2)]; %Probe feed - monopole end V(Index)=1.0*EdgeLength(Index); N=length(Index); %Identify phase shift %The progressive phase shift is zero for the broadside array: phase=0; %The progressive phase shift can be any number for %the end-fire array: %phase=-2*pi/3; % or -pi/2 or etc. %Identify feeding voltages-linear array (dipole array only!) for n=1:N nn=Index(n); V(nn)=V(nn)*exp(j*phase*(n-1)); end %Solve system of MoM equations tic; I=Z\V.'; toc %elapsed time %Terminal impedance (dipole array only!) for n=1:N nn=Index(n); GapCurrent(n)=I(nn)*EdgeLength(nn); GapVoltage(n)=V(nn)/EdgeLength(nn); Impedance (n)=GapVoltage(n)/GapCurrent(n); %this is the terminal impedance FeedPower (n)=1/2*real(GapCurrent(n)*conj(GapVoltage(n))); end Impedance FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_',... 'I','V','GapCurrent','GapVoltage','Impedance','FeedPower'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter06/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index=find(t(4,:)1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter07/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter07/rwg3.m%RWG3-FREQUENCY LOOP % Calculates the impedance matrix using function IMPMET % and solves MoM equations % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Number of frequency steps NumberOfSteps % Lower frequency FreqStart % Upper frequency FreqStop % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/25 % Chapter 7 clear all %Load the data load('mesh2'); %Frequency series parameters NumberOfSteps=200; FreqStart =25e6; %in Hz FreqStop =500e6; %in Hz step=(FreqStop-FreqStart)/(NumberOfSteps-1); %EM parameters epsilon_ =8.854e-012; mu_ =1.257e-006; %speed of light c_=1/sqrt(epsilon_*mu_); %free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - metal impedance matrix for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end %Frequency series T0=cputime; for FF=1:NumberOfSteps FF f(FF) =FreqStart+step*(FF-1); omega =2*pi*f(FF); k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; FactorA =FactorA.'; FactorFi =FactorFi.'; Z = impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %Find the feeding edge(s)(closest to the origin) FeedPoint=[0; 0; 0]; for m=1:EdgesTotal V(m)=0; Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-FeedPoint; end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=INDEX(1); %Center feed - dipole %Index=INDEX(1:2); %Probe feed - monopole %Solution of MoM equations V(Index)=1*EdgeLength(Index); I=Z\V.'; CURRENT(:,FF)=I(:); %Impedance GapCurrent(FF) =sum(I(Index).*EdgeLength(Index)'); GapVoltage(FF) =mean(V(Index)./EdgeLength(Index)); Impedance(FF) =GapVoltage(FF)/GapCurrent(FF); FeedPower(FF) =1/2*real(GapCurrent(FF)*conj(GapVoltage(FF))); Imp =Impedance(FF) T=cputime-T0 end %Save result FileName='current.mat'; save(FileName, 'f','NumberOfSteps','FreqStart','FreqStop','step',... 'omega','mu_','epsilon_','c_', 'eta_',... 'CURRENT','GapCurrent','GapVoltage','Impedance','FeedPower','Index'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter07/rwg5single.m%RWG5SINGLE Visualizes the surface current magnitude % at a particular frequency % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, current.mat, % from RWG3 as inputs. % % The following parameters need to be specified: % % Frequency value [Hz] FreqToPlot % % Copyright 2002 AEMM. Revision 2002/03/23 % Chapter 7 clear all load('mesh2'); load('current'); FreqToPlot=75e6 [dummy,FF]=min(abs(FreqToPlot-f)); f(FF) I(:,1)=CURRENT(:,FF); Index=find(t(4,:)1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter08/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter08/rwg31.m%RWG31 FREQUENCY LOOP % Calculates the impedance matrix using function IMPMET % and solves MoM equations % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Number of frequency steps NumberOfSteps % Lower frequency FreqStart % Upper frequency FreqStop % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % This script is equivalent to the script RWG3 from Chapter 7 % % Copyright 2002 AEMM. Revision 2002/03/25 % Chapter 8 clear all %Load the data load('mesh2'); %Frequency series parameters NumberOfSteps=500; FreqStart =12.5e6; %in Hz FreqStop =6250e6; %in Hz step=(FreqStop-FreqStart)/(NumberOfSteps-1); %EM parameters epsilon_ =8.854e-012; mu_ =1.257e-006; %speed of light c_=1/sqrt(epsilon_*mu_); %free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - metal impedance matrix for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end %Frequency series T0=cputime; for FF=1:NumberOfSteps FF f(FF) =FreqStart+step*(FF-1); omega =2*pi*f(FF); k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; FactorA =FactorA.'; FactorFi =FactorFi.'; Z = impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %Find the feeding edge(s)(closest to the origin) FeedPoint=[0; 0; 0]; for m=1:EdgesTotal V(m)=0; Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-FeedPoint; end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=INDEX(1); %Center feed - dipole %Index=INDEX(1:2); %Probe feed - monopole %Solution of MoM equations V(Index)=1*EdgeLength(Index); I=Z\V.'; CURRENT(:,FF)=I(:); %Impedance GapCurrent(FF) =sum(I(Index).*EdgeLength(Index)'); GapVoltage(FF) =mean(V(Index)./EdgeLength(Index)); Impedance(FF) =GapVoltage(FF)/GapCurrent(FF); FeedPower(FF) =1/2*real(GapCurrent(FF)*conj(GapVoltage(FF))); Imp =Impedance(FF) T=cputime-T0 end %Save result FileName='current.mat'; save(FileName, 'f','NumberOfSteps','FreqStart','FreqStop','step',... 'omega','mu_','epsilon_','c_', 'eta_',... 'CURRENT','GapCurrent','GapVoltage','Impedance','FeedPower','Index'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter08/rwg32.m%RWG32 FREQUENCY LOOP % Calculates radiated field at a point as a function % of frequency % The point is outside the metal surface % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG31 as inputs. % % The following parameters need to be specified: % % Observation point ObservationPoint[X; Y; Z] (m) % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 8 clear all %Load the data load('mesh2'); load('current'); ObservationPoint=[0 0 1]'; %Observation point %FRequency series for FF=1:NumberOfSteps FF omega=2*pi*f(FF); k=omega/c_; K=j*k; for m=1:EdgesTotal Point1=Center(:,TrianglePlus(m)); Point2=Center(:,TriangleMinus(m)); DipoleCenter(:,m)=0.5*(Point1+Point2); DipoleMoment(:,m)=EdgeLength(m)*CURRENT(m,FF)*(-Point1+Point2); end [EP,HP]=point(ObservationPoint,eta_,K,DipoleMoment,DipoleCenter); EField=sum(EP,2); HField=sum(HP,2); E(1:3,FF)=EField; end FileName=strcat('radiatedfield.mat'); save(FileName,'f','E','ObservationPoint'); plot(f,abs(E(2,:)),f,abs(E(1,:)),f,abs(E(3,:))) grid on xlabel('f,Hz') ylabel('Electric field magnitude (three components, V/m)') b=figure; plot(f,unwrap(angle(E(2,:)))) grid on xlabel('f,Hz') ylabel('Phase of Ey, rad') Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter08/rwg33.m%RWG33 FREQUENCY LOOP FOR THE SCATTERING PROBLEM % Calculates the voltage in the feed of the receiving % antenna if the incident signal (E-field) is given. % Performs calculations at all required frequencies % % Uses mesh2.mat (created by RWG2), current.mat (created by RWG31), % and radiatedfield.mat (created by RWG32) as inputs % % The following parameters need to be specified prior to % calculations: % % Number of frequency steps NumberOfSteps % Lower frequency FreqStart % Upper frequency FreqStop % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % Wave vector kv % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 8 %load the data load('mesh2'); load('current') load('radiatedfield'); %Frequency series parameters NumberOfSteps=500; FreqStart =12.5e6; %in Hz FreqStop =6250e6; %in Hz step=(FreqStop-FreqStart)/(NumberOfSteps-1); %EM parameters epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - metal impedance matrix for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end %Frequency series T0=cputime; for FF=1:NumberOfSteps FF f(FF) =FreqStart+step*(FF-1); omega =2*pi*f(FF); k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; FactorA =FactorA.'; FactorFi =FactorFi.'; Z = impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); d=[0 0 -1]; kv=k*d; Pol=E(1:3,FF).'; for m=1:EdgesTotal ScalarProduct =sum(kv.*Center(:,TrianglePlus(m))'); EmPlus =Pol.'*exp(-j*ScalarProduct); ScalarProduct =sum(kv.*Center(:,TriangleMinus(m))'); EmMinus =Pol.'*exp(-j*ScalarProduct); ScalarPlus =sum(EmPlus.* RHO_Plus(:,m)); ScalarMinus =sum(EmMinus.*RHO_Minus(:,m)); V(m) =EdgeLength(m)*(ScalarPlus/2+ScalarMinus/2); end %Solution of MoM equations I=Z\V.'; %Find the antenna feed position (receiving antenna) FeedPoint=[0; 0; 0]; for m=1:EdgesTotal V(m)=0; Distance(:,m)=0.5*sum(p(:,Edge_(:,m)),2)-FeedPoint; end [Y,INDEX]=sort(sum(Distance.*Distance)); Index=INDEX(1); %Center feed - dipole %Index=INDEX(1:2); %Probe feed - monopole %Find the voltage across the antenna feed of the receiving %antenna Imp=Impedance(FF); FeedCurReceived =I(Index)*EdgeLength(Index); FeedVolReceived =FeedCurReceived*Imp; OutputVoltage(FF) =FeedVolReceived; PowerConjMatch(FF) =1/8*(FeedVolReceived*conj(FeedVolReceived))/real(Imp); T=cputime-T0 end %Save result FileName='receivedfield.mat'; save(FileName, 'f','NumberOfSteps','FreqStart','FreqStop','OutputVoltage','PowerConjMatch',... 'Impedance','Index'); plot(f,PowerConjMatch); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter08/rwg34.m%RWG34 Calculation of the received pulse using % antenna-to-antenna transfer function % % Uses receivedfield.mat (created by RWG33)as an input % % The following parameters need to be specified prior to % calculations: % % Duration of the primary Gaussian voltage pulse (s) Duration % Time offset to plot the received and the (delayed) % transmitted pulse on the same figure including the % "pulse tail" Offset % % Note: This script is not optimized and is rather slow % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 8 clear all load('receivedfield') %Identify the pulse parameters: Duration=1.00e-9; %in seconds sigma=Duration/7; %Gaussian pulse Offset=500*sigma; %Offset %Sampling time Ts=1/(2*max(f)) %Length of DFT N=2*NumberOfSteps; %Obtain the pulse form i=1; Time=-Duration; t(1)=Time; p(1)=t(1)*exp(-t(1)^2/(2*sigma^2))/sigma^2; while (Time1); t(:,Cut)=[]; for m=1:length(t) n=t(1:3,m); X(1:3,m)=p(1,n)'; Y(1:3,m)=p(2,n)'; Z(1:3,m)=p(3,n)'; end g=fill3(X,Y,Z,Color); axis('equal') rotate3d on Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter09/point.mfunction[EField, HField]=... point1(Point,eta_,K,DipoleMoment,DipoleCenter) %POINT Radiated/scattered field at a point of a dipole array % or a single dipole. Gives exact near- and far-fields. Outputs % individual contribution of each dipole. % % Observation point Point(1:3) % Array of dipole moments DipoleMoment(1:3,1:EdgesTotal) % Array of dipole centers DipoleCenter(1:3,1:EdgesTotal) % E-field at the observation point E(1;3,1:EdgesTotal) % H-field at the observation point H(1;3,1:EdgesTotal) % % Copyright 2002 AEMM. Revision 2002/03/11 % Chapter 3 C=4*pi; ConstantH=K/C; ConstantE=eta_/C; m=DipoleMoment; c=DipoleCenter; r =repmat(Point,[1 length(c)])-c(1:3,:); PointRM =repmat(sqrt(sum(r.*r)),[3 1]); EXP =exp(-K*PointRM); PointRM2=PointRM.^2; C=1./PointRM2.*(1+1./(K*PointRM)); D=repmat(sum(r.*m),[3 1])./PointRM2; M=D.*r; HField=ConstantH*cross(m,r).*C.*EXP; EField=ConstantE*((M-m).*(K./PointRM+C)+2*M.*C).*EXP; Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter09/rwg1.m%RWG1 Geometry calculations - all Chapters % Uses the structure mesh file, e.g. platefine.mat, % as an input. % % Creates the RWG edge element for every inner edge of % the structure. The total number of elements is EdgesTotal. % Outputs the following arrays: % % Edge first node number Edge_(1,1:EdgesTotal) % Edge second node number Edge_(2,1:EdgesTotal) % Plus triangle number TrianglePlus(1:EdgesTotal) % Minus triangle number TriangleMinus(1:EdgesTotal) % Edge length EdgeLength(1:EdgesTotal) % Edge element indicator EdgeIndicator(1:EdgesTotal) % % Also outputs areas and midpoints of separate triangles: % Triangle area Area(1:TrianglesTotal) % Triangle center Center(1:TrianglesTotal) % % This script may handle surfaces with T-junctions % including monopoles over various metal surfaces and % certain metal meshes % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all tic; load('mesh/strip'); Feed=[]; [s1 s2]=size(p); if(s1==2) p(3,:)=0; %to convert 2D to 3D end %Eliminate unnecessary triangles Remove=find(t(4,:)>1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter09/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % The modification for Chapters 6 and 7 passes the array % Feed (see Chapter 6) % % Copyright 2002 AEMM. Revision 2002/03/14 Chapter 6 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Feed Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter09/rwg3.m%RWG3-FREQUENCY LOOP % Calculates the impedance matrix using function IMPMET % and solves MoM equations. % Takes into account lumped elements (L, C, R) - antenna load % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Number of frequency steps NumberOfSteps % Lower frequency FreqStart % Upper frequency FreqStop % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Number of lumped elements LNumber % Lumped element locations [x y z] LoadPoint % Vector of L, C, and R [L C R] LoadValue % "Direction" of lumped element [x y z] LoadDir % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 9 clear all %Load the data load('mesh2'); %Frequency series parameters NumberOfSteps=200; FreqStart =25e6; %in Hz FreqStop =500e6; %in Hz step=(FreqStop-FreqStart)/(NumberOfSteps-1); %EM parameters epsilon_ =8.854e-012; mu_ =1.257e-006; %speed of light c_=1/sqrt(epsilon_*mu_); %free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - metal impedance matrix for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end %Frequency series T0=cputime; for FF=1:NumberOfSteps FF f(FF) =FreqStart+step*(FF-1); omega =2*pi*f(FF); k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; FactorA =FactorA.'; FactorFi =FactorFi.'; Z = impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); %Lumped impedance format: % LoadPoint Lumped element locations % LoadValue Vector of L, C, and R % LoadDir "Direction" of lumped element %This block introduces a few lumped elements %When LNumber is zero then skip loading LNumber=2; LoadPoint(1:3,1)=[0 0.50 0]'; LoadValue(1:3,1)=[0 1e+64 100]'; %LCR LoadDir (1:3,1)=[0 1 0]'; LoadPoint(1:3,2)=[0 -0.50 0]'; LoadValue(1:3,2)=[0 1e+64 100]'; %LCR LoadDir (1:3,2)=[0 1 0]'; for k=1:LNumber DeltaZ(k)=j*omega*LoadValue(1,k)+1/(j*omega*LoadValue(2,k))+... LoadValue(3,k); end %Impedance matrix modification L=0; for k=1:LNumber for m=1:EdgesTotal n1=Edge_(1,m); n2=Edge_(2,m); Dist(m) =norm(0.5*sum(p(:,Edge_(:,m)),2)-LoadPoint(:,k)); EdgeVector =(p(:,n1)-p(:,n2))/EdgeLength(m); Orien (m) =abs(dot(EdgeVector,LoadDir(:,k))); end [Y,INDEX1] = sort(Dist); for n=1:EdgesTotal Index1 =INDEX1(n); q =Orien(Index1); if(q2); %Here's the difference of Chapter 10 t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have two common points n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); % 0 - metal ground plane % 1 - feeding bottom edges % 2 - connecting strip (probe feed) % 3 - feeding top edges % 4 - patch area %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center ... h %This is new (patch height) %View result (make sure everything is ok) cd mesh viewer('patch') cd .. Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter10/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/17 % Chapter 2/Chapter 10 clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... h %This is new (patch height) Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/chapter10/rwg3.m%RWG3 FREQUENCY LOOP % Calculates the impedance matrix using function IMPMET % and solves MoM equations % Uses the mesh file from RWG2, mesh2.mat, as an input. % Includes three additional impedance matrixes: % Dielectric-to-dielectric impedance matrix ZDD % Dielectric-to-metal impedance matrix ZDS % Metal-to-dielectric impedance matrix ZSD % % The following parameters need to be specified prior to % calculations: % % Number of frequency steps NumberOfSteps % Lower frequency FreqStart % Upper frequency FreqStop % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % Relative dielectric constant of % the dielectric material epsilon_R % % The algorithm of Chapter 10 uses the transpose impedance matrix % Z.' instead of Z for the metal structure. To better see what is % the physical difference between Z and Z.' one should again return % to the basic paper Rao, Wilton, Glisson, IEEE Trans. Antennas and % Propagation, vol. AP-30, No 3, pp. 409-418, 1982. It can be shown % that Z utilizes 1 integration point for the outer integral (m) and % 9 integration points for the inner integral (n) whereas its transpose % implies 9 integration points for the outer integral and 1 point for % the inner integral. The difference between Z and Z.' is really % insignificant for the pure metallic antennas. Either Z or Z.' can % be used. % However, for the metal-dielectric structures, Z.' may be preferred % since it better matches other mixed integrals (metal-dielectric). % % Copyright 2002 AEMM. Revision 2002/03/23 % Chapter 10 clear all %Load the data load('mesh2'); %Frequency series parameters NumberOfSteps=21; FreqStart =1.0e9; %in Hz FreqStop =6.0e9; %in Hz step=(FreqStop-FreqStart)/(NumberOfSteps-1); %EM parameters epsilon_ =8.854e-012; epsilon_R =1.0; mu_ =1.257e-006; %speed of light c_=1/sqrt(epsilon_*mu_); %free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - metal impedance matrix for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end %Contemporary variables - dielectric/metal impedance matrix DP =find(t(4,:)==0); M =length(DP); %first M triangles form the dielectric plate delta =[0;0;1e-12]; %introduce to avoid singularity 0/0 Middle=[0; 0; h/2]; for m=1:M N=t(1:3,m); Point(:,m)=Center(:,m) +Middle; IMT(:,:,m)=Center_(:,:,m) +repmat(Middle,[1,9]); end %Frequency series for FF=1:NumberOfSteps tic; FF f(FF) =FreqStart+step*(FF-1); omega =2*pi*f(FF); k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; FactorA =FactorA.'; FactorFi =FactorFi.'; %Metal-to-metal impedance matrix (SS) ZSS= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi); ZSS=ZSS.'; toc; tic; %Dielectric-to-dielectric impedance matrix (DD) ZDD = zeros(M,M)+j*zeros(M,M); for m=1:M OP =Point(:,m)+delta; IP =Point; E =point_(OP,K,k,Constant2,Area(1:M),h,IP); ZDD(m,:)=E(3,:); %Non-diagonal terms IP =IMT(:,:,m); E =point_(OP,K,k,Constant2,ones(1,9)*Area(m)/9,h,IP); ZDD(m,m)=sum(E(3,:)); %Diagonal terms end ZDD=j*omega*epsilon_*(epsilon_R-1)*ZDD; %Non-symmetric matrix since areas are different toc; tic; %Dielectric-to-metal impedance matrix (DS) ZDS = zeros(EdgesTotal,M)+j*zeros(EdgesTotal,M); for m=1:EdgesTotal OPPlus =Center(:,TrianglePlus(m))+delta; OPMinus =Center(:,TriangleMinus(m))+delta; EP =point_(OPPlus,K,k,Constant2,Area(1:M),h,Point); EM =point_(OPMinus,K,k,Constant2,Area(1:M),h,Point); ScalarPlus =sum(EP.* repmat(RHO_Plus(:,m), [1 M])); ScalarMinus =sum(EM.* repmat(RHO_Minus(:,m),[1 M])); ZDS(m,:) =EdgeLength(m)*(ScalarPlus/2+ScalarMinus/2); end toc; tic; %Metal-to-dielectric impedance matrix (SD) ZSD = zeros(M,EdgesTotal)+j*zeros(M,EdgesTotal); C1=1/(2*epsilon_)/(-j*omega); C2=j*omega*epsilon_*(epsilon_R-1); C=C1*C2; for m=1:M Q=sum(abs([Center(1,m)-Center(1,:); Center(2,m)-Center(2,:)])); T=find(Q1); t(:,Remove)=[]; TrianglesTotal=length(t); %Find areas of separate triangles for m=1:TrianglesTotal N=t(1:3,m); Vec1=p(:,N(1))-p(:,N(2)); Vec2=p(:,N(3))-p(:,N(2)); Area(m) =norm(cross(Vec1,Vec2))/2; Center(:,m)=1/3*sum(p(:,N),2); end %Find all edge elements "Edge_" with at least two %adjacent triangles Edge_=[]; n=0; for m=1:TrianglesTotal N=t(1:3,m); for k=m+1:TrianglesTotal M=t(1:3,k); a=1-all([N-M(1) N-M(2) N-M(3)]); if(sum(a)==2) %triangles m and k have common edge n=n+1; Edge_=[Edge_ M(find(a))]; TrianglePlus(n)=m; TriangleMinus(n)=k; end; end end EdgesTotal=length(Edge_); %This block is only meaningful for T junctions %It leaves only two edge elements at a junction Edge__=[Edge_(2,:); Edge_(1,:)]; Remove=[]; for m=1:EdgesTotal Edge_m=repmat(Edge_(:,m),[1 EdgesTotal]); Ind1=any(Edge_ -Edge_m); Ind2=any(Edge__ -Edge_m); A=find(Ind1.*Ind2==0); if(length(A)==3) %three elements formally exist at a junction Out=find(t(4,TrianglePlus(A))==t(4,TriangleMinus(A))); Remove=[Remove A(Out)]; end end Edge_(:,Remove) =[]; TrianglePlus(Remove) =[]; TriangleMinus(Remove) =[]; EdgesTotal=length(Edge_) %All structures of this chapter have EdgeIndicator=2 EdgeIndicator=t(4,TrianglePlus)+t(4,TriangleMinus); %Find edge length for m=1:EdgesTotal EdgeLength(m)=norm(p(:,Edge_(1,m))-p(:,Edge_(2,m))); end toc %Save result save mesh1 p ... t ... Edge_ ... TrianglesTotal ... EdgesTotal ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... Center Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/m_appendix_b/rwg2.m%RWG2 Geometry calculations - all Chapters % Uses the mesh file from RWG1, mesh1.mat, % as an input. % % Creates the following parameters of the RWG edge elements: % % Position vector rho_c_plus from the free vertex of % the "plus" triangle to its center % RHO_Plus(1:3,1:EdgesTotal) % Position vector rho_c_minus from the center of the "minus" % triangle to its free vertex % RHO_Minus(1:3,1:EdgesTotal) % % In addition to these parameters creates the following % arrays for nine subtriangles (barycentric subdivision): % % Midpoints of nine subtriangles % Center_(1:3,1:9,1:TrianglesTotal) % Position vectors rho_c_plus from the free vertex of % the "plus" triangle to nine subtriangle midpoints % RHO__Plus(1:3,1:9,1:EdgesTotal) % Position vectors rho_c_minus from nine subtriangle midpoints % to the free vertex of the "minus" triangle % RHO__Minus(1:3,1:9,1:EdgesTotal) % % See Rao, Wilton, Glisson, IEEE Trans. Antennas and Propagation, % vol. AP-30, No 3, pp. 409-418, 1982. % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 2/Appendix B clear all %load the data load('mesh1') %Find nine sub-triangle midpoints IMT=[]; for m=1:TrianglesTotal n1=t(1,m); n2=t(2,m); n3=t(3,m); M=Center(:,m); r1= p(:,n1); r2= p(:,n2); r3= p(:,n3); r12=r2-r1; r23=r3-r2; r13=r3-r1; C1=r1+(1/3)*r12; C2=r1+(2/3)*r12; C3=r2+(1/3)*r23; C4=r2+(2/3)*r23; C5=r1+(1/3)*r13; C6=r1+(2/3)*r13; a1=1/3*(C1+C5+r1); a2=1/3*(C1+C2+M); a3=1/3*(C2+C3+r2); a4=1/3*(C2+C3+M); a5=1/3*(C3+C4+M); a6=1/3*(C1+C5+M); a7=1/3*(C5+C6+M); a8=1/3*(C4+C6+M); a9=1/3*(C4+C6+r3); Center_(:,:,m)=... [a1 a2 a3 a4 a5 a6 a7 a8 a9]; %Analytical expression for the self-coupling terms %See Eibert and Hansen, IEEE Trans. Antennas and Propagation %vol. AP-43, No 12, pp. 1499-1502, 1995. a=sum(r13.*r13); b=sum(r13.*r23); c=sum(r23.*r23); d=a-2*b+c; N1=(a-b+sqrt(a)*sqrt(d))*(b+sqrt(a)*sqrt(c)); D1=(-a+b+sqrt(a)*sqrt(d))*(-b+sqrt(a)*sqrt(c)); N2=(-b+c+sqrt(c)*sqrt(d))*(b+sqrt(a)*sqrt(c)); D2=(b-c+sqrt(c)*sqrt(d))*(-b+sqrt(a)*sqrt(c)); N3=(a-b+sqrt(a)*sqrt(d))*(-b+c+sqrt(c)*sqrt(d)); D3=(b-c+sqrt(c)*sqrt(d))*(-a+b+sqrt(a)*sqrt(d)); Integral(m)=1/6*(1/sqrt(a)*log(N1/D1) +1/sqrt(c)*log(N2/D2) +1/sqrt(d)*log(N3/D3)); Integral(m)=4*Integral(m); end %PLUS for m=1:EdgesTotal NoPlus=TrianglePlus(m); n1=t(1,NoPlus); n2=t(2,NoPlus); n3=t(3,NoPlus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Plus(:,m) =+Center(:,NoPlus)-FreeVertex; %Nine rho's of the "plus" triangle RHO__Plus(:,:,m) =... +Center_(:,:,NoPlus)-repmat(FreeVertex,[1 9]); end %MINUS for m=1:EdgesTotal NoMinus=TriangleMinus(m); n1=t(1,NoMinus); n2=t(2,NoMinus); n3=t(3,NoMinus); if((n1~=Edge_(1,m))&(n1~=Edge_(2,m))) NODE=n1; end; if((n2~=Edge_(1,m))&(n2~=Edge_(2,m))) NODE=n2; end; if((n3~=Edge_(1,m))&(n3~=Edge_(2,m))) NODE=n3; end; FreeVertex=p(:,NODE); RHO_Minus(:,m) =-Center(:,NoMinus) +FreeVertex; %Nine rho's of the "minus" triangle RHO__Minus(:,:,m)=... -Center_(:,:,NoMinus)+repmat(FreeVertex,[1 9]); end %Save result save mesh2 p ... t ... TrianglesTotal ... EdgesTotal ... Edge_ ... TrianglePlus ... TriangleMinus ... EdgeLength ... EdgeIndicator ... Area ... RHO_Plus ... RHO_Minus ... RHO__Plus ... RHO__Minus ... Center ... Center_ ... Integral Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/m_appendix_b/rwg3.m%RWG3 Calculates the impedance matrix using function IMPMET % Uses the mesh file from RWG2, mesh2.mat, as an input. % % The following parameters need to be specified prior to % calculations: % % Frequency (Hz) f % Dielectric constant (SI) epsilon_ % Magnetic permeability (SI) mu_ % % Copyright 2002 AEMM. Revision 2002/03/26 % Chapter 2/Appendix B clear all %Load the data load('mesh2'); %EM parameters (f=3e8 means that f=300 MHz) f =3e8; epsilon_ =8.854e-012; mu_ =1.257e-006; %Speed of light c_=1/sqrt(epsilon_*mu_); %Free-space impedance eta_=sqrt(mu_/epsilon_); %Contemporary variables - introduced to speed up %the impedance matrix calculation omega =2*pi*f; k =omega/c_; K =j*k; Constant1 =mu_/(4*pi); Constant2 =1/(j*4*pi*omega*epsilon_); Factor =1/9; FactorA =Factor*(j*omega*EdgeLength/4)*Constant1; FactorFi =Factor*EdgeLength*Constant2; for m=1:EdgesTotal RHO_P(:,:,m)=repmat(RHO_Plus(:,m),[1 9]); %[3 9 EdgesTotal] RHO_M(:,:,m)=repmat(RHO_Minus(:,m),[1 9]); %[3 9 EdgesTotal] end FactorA=FactorA.'; FactorFi=FactorFi.'; %Impedance matrix Z tic; %start timer Z= impmet( EdgesTotal,TrianglesTotal,... EdgeLength,K,... Center,Center_,... TrianglePlus,TriangleMinus,... RHO_P,RHO_M,... RHO__Plus,RHO__Minus,... FactorA,FactorFi,Integral); toc %elapsed time %Save result FileName='impedance.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','Z'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/m_appendix_b/rwg4.m%RWG4 Solves MoM equations for the scattering problem % Uses the mesh file from RWG2, mesh2.mat, and % the impedance file from RWG3, impedance.mat, % as inputs. % % Also calculates the "voltage" vector V (the right- % hand side of moment equations) % V(1:EdgesTotal) % % The following parameters need to be specified: % % Direction of the incident signal in Cartesian coordinates % d(1:3); % Direction of the E-field in the incident plane wave % in Cartesian coordinates Pol(1:3); % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('impedance'); %Incident field %Example: d=[0 0 -1] means that the incident signal % is in the -z direction. %Plate - normal incidence d =[0 0 -1]; Pol =[1 0 0]; %Dipole - normal incidence %d =[0 0 -1]; %Pol =[0 1 0]; %Custom incidence (example) %d =[1 0 0] %Pol =[0 -0.0037-0.0055*j 0] k=omega/c_; kv=k*d; for m=1:EdgesTotal ScalarProduct=sum(kv.*Center(:,TrianglePlus(m))'); EmPlus =Pol.'*exp(-j*ScalarProduct); ScalarProduct=sum(kv.*Center(:,TriangleMinus(m))'); EmMinus=Pol.'*exp(-j*ScalarProduct); ScalarPlus =sum(EmPlus.* RHO_Plus(:,m)); ScalarMinus=sum(EmMinus.*RHO_Minus(:,m)); V(m)=EdgeLength(m)*(ScalarPlus/2+ScalarMinus/2); end tic; %System solution I=Z\V.'; toc %elapsed time FileName='current.mat'; save(FileName, 'f','omega','mu_','epsilon_','c_', 'eta_','I','V','d','Pol'); Wiley_-_Antenna_and_EM_Modeling_with_MATLAB/Wiley - Antenna and EM Modeling with MATLAB/m_files/m_appendix_b/rwg5.m%RWG5 Visualizes the surface current magnitude % Uses the mesh file from RWG2, mesh2.mat, and % the file containing surface current coefficients, % current.mat, from RWG4 as inputs. % % Copyright 2002 AEMM. Revision 2002/03/05 % Chapter 2 clear all %load the data load('mesh2'); load('current'); Index=find(t(4,:)