?? mom.m
字號:
%C*****************************************************************
%C MoM
%C*****************************************************************
%C This is a MATLAB based program using
%C
%C I. POCKLINGTON'S INTEGRAL EQUATION (8-24)
%C II. HALLEN'S INTEGRAL EQUATION (8-27)
%C
%C That computes:
%C
%C A. CURRENT DISTRIBUTION
%C B. INPUT IMPEDANCE
%C C. NORMALIZED AMPLITUDE RADIATION PATTERN
%C
%C Of a linear symmetrically-excited dipole.
%C
%C THIS PROGRAM USES PULSE EXPANSION FOR THE ELECTRIC CURRENT MODE
%C AND POINT-MATCHING FOR THE ELECTRIC FIELD AT THE CENTER OF EACH
%C WIRE SEGMENT.
%C
%C DELTA-GAP FEED MODEL IS USED IN BOTH FORMULATIONS. IN ADDITION,
%C MAGNETIC-FRILL GENERATOR IS AVAILABLE IN THE POCKLINGTON'S
%C INTEGRAL EQUATION.
%C
%C OPTION I. POCKLINGTON'S INTEGRAL EQUATION
%C OPTION II. HALLEN'S INTEGRAL EQUATION
%C
%C
%C ** INPUT DATA:
%C
%C TL = TOTAL DIPOLE LENGTH (IN WAVELENGTHS)
%C RA = RADIUS OF THE WIRE (IN WAVELENGTHS)
%C NM = TOTAL NUMBER OF SUBSECTIONS (MUST BE AN ODD INTEGER)
%C IEX = OPTION TO USE EITHER MAGNETIC-FRILL GENERATOR OR DELTA GAP
%C
%C IEX = 1 : MAGNETIC-FRILL GENERATOR
%C IEX = 2 : DELTA-GAP FEED
%C
%C ** NOTE: IGNORE INPUT PARAMETER IEX WHEN CHOOSING OPTION II
%C (i.e., HALLEN'S FORMULATION)
%C *****************************************************************
clear;
close;
%*********************************************************
%*** OUTPUT DEVICE OPTION # : 1 ---> SCREEN
% 2 ---> OUTPUT FILE
%*********************************************************
disp('OUTPUT DEVICE OPTION FOR THE OUTPUT PARAMETERS');
disp(' OPTION (1):SCREEN');
disp(' OPTION (2):OUTPUT FILE');
option_a=0; option_b=0; filename=0;
while((option_a~=1)&(option_a~=2))
option_a=input('OUTPUT DEVICE =');
if(option_a==2)
filename=input('INPUT THE DESIRED OUTPUT FILENAME <in single quotes> = ','s');
end
end
%*********************************************************
%*** READING OPTION # : 1 ---> POCKLINGTON'S EQUATIONS
% 2 ---> HALLEN'S EQUATION
%*********************************************************
disp('CHOICE OF POCKLINGTON''S OR HALLEN''S EQN.');
disp(' OPTION (1):POCKLINGTON''S EQN.');
disp(' OPTION (2):HALLEN''S EQN.');
option=input('OPTION NUMBER =');
switch option
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%***********************************************************************
%*** %POCKLINGTON'S INTEGRAL EQUATIONS ***
%***********************************************************************
case 1
%******************************************************
% SOME CONSTANTS AND ARRAYS
%******************************************************
nmax=200;
nmt=2*nmax-1;
cgaw=zeros(nmax);
cga=cgaw(1,1:nmax);
zmnw=zeros(nmt);
zmn=zmnw(1,1:nmt);
waw=zeros(nmt);
wa=waw(1,1:nmt);
etm=zeros(361);
etmm=etm(1,1:361);
pi=3.14159265;
rad=pi/180;
beta=2.0*pi;
eta=120*pi;
%******************************************************
%*** INPUT DATA
%******************************************************
nm=input('NUMBER OF SUBDIVISIONS <ODD NUMBER> =');
tl=input('TOTAL DIPOLE LENGTH <WAVELENGTHS> =');
ra=input('RADIUS OF DIPOLE <WAVELENGTHS> =');
%*********************************************************
%*** EXITATION OPTION # : 1 ---> MAGNETIC-FRILL
% 2 ---> DELTA-GAP
%**********************************************************
disp('EXCITATION :');
disp(' OPTION (1):MAGNETIC-FRILL');
disp(' OPTION (2):DELTA-GAP');
iex=input('OPTION NUMBER :');
hl=tl/2;
nmh=0.5*(nm+1);
dz=2*hl/nm;
zm=hl-0.5*dz;
b=0.5*dz;
a=-0.5*dz;
n=79;
hd=(b-a)/(n+1);
%*********************************************************************
% THE IMPEDANCE MATRIX HAS A TOEPLITZ PROPERTY, THEREFORE ONLY
% NM ELEMENTS NEED TO BE COMPUTED, AND THE MATRIX IS FILLED IN A
% FORM THAT CAN BE SOLVED BY A TOEPLITZ MATRIX SOLVING
%**********************************************************************
for I = 1: nm
zn=hl-(I-0.5)*dz;
za1=zn-zm+a;
%******************************************************************
% FAST ALGORITHM FORM OF THE SIMPSON'S INTEGRAL ROUTINE
%******************************************************************
recgp=sqrt(ra*ra+za1*za1);
cgp1=exp(-j*beta*recgp)*((1.0+j*beta*recgp)*(2.0*recgp*recgp-3.0*ra*ra)+(beta*ra*recgp)^2)/(2.0*beta*recgp^5);
zb1=zn-zm+b;
roc=sqrt(ra*ra+zb1*zb1);
cgp2=exp(-j*beta*roc)*((1.0+j*beta*roc)*(2.0*roc*roc-3.0*ra*ra)+(beta*ra*roc)^2)/(2.0*beta*roc^5);
crt=cgp1+cgp2;
for k = 1: n
xk=a+k*hd;
zx1=zn-zm+xk;
r=sqrt(ra*ra+zx1*zx1);
cgp3=exp(-j*beta*r)*((1.0+j*beta*r)*(2.0*r*r-3.0*ra*ra)+(beta*ra*r)^2)/(2.0*beta*r^5);
if mod(k,2)~=0
crt=crt+4.0*cgp3;
else
crt=crt+2.0*cgp3;
end
end
crt=crt*hd*0.33333;
zmn(I)=crt;
if I~=1
zmn(nm+I-1)=crt;
end
end
rb=2.3*ra;
tlab=2.0*log(2.3);
for i = 1: nm
zi=hl-(i-0.5)*dz;
r1=beta*sqrt(zi*zi+ra*ra);
r2=beta*sqrt(zi*zi+rb*rb);
j=0.0+j*1.0;
if iex==1
cga(i)=-j*beta^2/(eta*tlab)*(expm(-j*r1)/r1-expm(-j*r2)/r2);
else
if i~=nmh
cga(i)=0;
else
cga(i)=-j*beta/(eta*dz);
end
end
end
lolo = cga;
%******************************************************************
% INPUT:
% (C)A(2*M - 1) THE FIRST ROW OF THE T-MATRIX FOLLOWED BY
% ITS FIRST COLUMN BEGINNING WITH THE SECOND
% ELEMENT. ON RETURN A IS UNALTERED.
% (C)B(M) THE RIGHT HAND SIDE VECTOR B.
% (C)WA(2*M-2) A WORK AREA VECTOR
% (I)M ORDER OF MATRIX A.
% OUTPUT:
% (C)B(M) THE SOLUTION VECTOR.
% PURPOSE:
% SOLVE A SYSTEM OF EQUATIONS DESCRIBED BY A TOEPLITZ MATRIX.
% A * X = B
% ******************************************************************
t2=length(zmn);
t1=length(wa);
m=nm;
a=zmn;
b=cga;
wa=wa;
m=nm;
a1=a;
a2=a(m+1:t2);
t=b;
c1=wa;
c2=wa(m-1:t1);
r1=a1(1);
t(1)=b(1)/r1;
if m ~=1
for n=2:m
n1=n-1;
n2=n-2;
r5=a2(n1);
r6=a1(n);
if n ~= 2
c1(n1)=r2;
for i1=1:n2
i2=n-i1;
r5=r5+a2(i1)*c1(i2);
r6=r6+a1(i1+1)*c2(i1);
end
end
r2=-r5/r1;
r3=-r6/r1;
r1=r1+r5*r3;
if n~=2
r6=c2(1);
c2(n1)=0.0+j*0.0;
for i1=2:n1
r5=c2(i1);
c2(i1)=c1(i1)*r3+r6;
c1(i1)=c1(i1)+r6*r2;
r6=r5;
end
end
c2(1)=r3;
r5=0.0+j*0.0;
for i1=1:n1
i2=n-i1;
r5=r5+a2(i1)*t(i2);
end
r6=(b(n)-r5)/r1;
for i1=1:n1
t(i1)=t(i1)+c2(i1)*r6;
end
t(n)=r6;
end
end
cga=t;
%******************************************************************
% Output files for curr-MoM.dat
%******************************************************************
fid=fopen('Curr-MoM_m.dat','wt');
fprintf(fid,'CURRENT DISTRIBUTION ALONG ONE HALF OF THE DIPOLE\n');
fprintf(fid,'POSITION Z: MAGNITUDE: REAL PART: IMAGINARY PART: SECTION #:\n\n');
%*****************************************************************************
% OUTPUT THE CURRENT DISTRIBUTION ALONG OF THE DIPOLE
%*****************************************************************************
for i=1:nmh
xi=hl-(i-0.5)*dz;
yi=abs(cga(i));
table2=[xi,yi,real(cga(i)),imag(cga(i)),i];
fprintf(fid,'%1.4f %1.6f %1.6f %1.6f %2.1f \n\n',table2);
end
%******************************************************************
% Output figure for curr-MoM.dat
%******************************************************************
i=[1:nmh];
xi=hl-(i-0.5)*dz;
yi=abs(cga(i));
stem(xi,yi)
grid on;
xlabel('Position (z)')
ylabel('Magnitude')
title('CURRENT DISTRIBUTION ALONG ONE HALF OF THE DIPOLE')
figure;
fclose(fid);
%******************************************************
% COMPUTATION OF THE INPUT IMPEDANCE
%******************************************************
zin=1.0/cga(nmh);
%******************************************************************
% COMPUTATION OF AMPLITUDE RADIATION PATTERN OF THE ANTENNA
%******************************************************************
for i=1:181
theta=(i-1.0)*rad;
cth=cos(theta);
sth=sin(theta);
if abs(cth)<0.001
ft=1;
else
ft=sin(beta*dz*cth*0.5)/(beta*dz*cth*0.5);
end
crt=0;
for m=1:nm
zm=hl-(m-0.5)*dz;
crt=crt+exp(j*beta*zm*cth)*ft*cga(m)*dz;
end
ptt=abs(crt)*sth*sth*eta*0.5;
etmm(i)=ptt;
end
amax=etmm(1);
for i=2:181
if etmm(i)>=amax
amax=etmm(i);
end
end
for i=1:181
ptt=etmm(i)/amax;
if ptt<=0.00001
ptt=0.00001;
etmm(i)=20*log10(ptt);
end
etmm(i)=20*log10(ptt);
end
%******************************************************************
% Output files for Patt-MoM_m.dat
%******************************************************************
pok=fopen('Patt-MoM_m.dat','wt');
fprintf(pok,'RADIATION PATTERN vs OBSERVATION ANGLE THETA\n');
fprintf(pok,'THETA (in degrees) MAGNITUDE(in dB)\n');
for i=1:181
xi=i-1;
table3=[xi,etmm(i)];
fprintf(pok,'%3.1f %3.3f \n',table3);
end
%******************************************************************
% Output figure for Patt-MoM_m.dat
%******************************************************************
i=[1:181];
xi=i-1;
% Polar Plot
etmm1=[etmm(1:181),fliplr(etmm(1:180))];
q=polar_dB([0:360],etmm1,-40,0,4,'-');
set(q,'linewidth',1.5);
% plot(xi,etmm(i),'linewidth',2);
% axis ([0 180 -60 0]);
% grid on;
% xlabel('Theta(degrees)');
% ylabel('Magnitude(dB)');
title('RADIATION PATTERN vs OBSERVATION ANGLE')
for i=182:361
xi=i-1;
table4=[xi,etmm(362-i)];
fprintf(pok,'%3.1f %3.3f \n',table4);
end
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -