gain2h.m

用matlab编写的电磁学方面的函数和图像

% gain2h.m - gain of 2D array of non-identical linear antennas with Hallen currents
%
% Usage: [ge,gh,th] = gain2h(L,d,I,N,ph0)   
%        [ge,gh,th] = gain2h(L,d,I,N)      (equivalent to ph0=0)
%
% L   = antenna lengths in wavelengths, L = [L1,L2,...,LK]
% d   = [x,y] or [x] locations of the K antennas, d must be Kx2 or Kx1 or 1xK
% I   = matrix of antenna currents, I is (2M+1)xK, see notes
% N   = number of azimuthal and polar angles over [0,2*pi]
% ph0 = azimuthal direction for E-plane pattern (in degrees)
%
% ge,gh = E-plane/H-plane gains at (N+1) polar or azimuthal angles over [0,2*pi]
% th    = (N+1) equally-spaced polar or azimuthal angles over [0,2*pi] in radians
%
% notes: the columns of I = [I1,I2,...,IK] are the currents on the K antennas,
%        evaluated at the equally-spaced points z(m)=m*Dz, m=-M:M, where h=L/2,  
%        and subject to the constraint that I(M,:)=I(-M,:)=0. The m-th current 
%        sample on the p-th antenna is Ip(m) = I(m,p). 
%
%        The current matrix I must be obtained from HALLEN3 or HALLEN4, as follows:
%        I = hallen3(L,a,d,V,M) or I = hallen4(L,a,d,V,M)
%
%        d is the matrix of the [x,y] locations of the antennas and is Kx2, that is, 
%        d = [x1,y1; x2,y2; ...; xK,yK]. If the antennas are along the x-axis then 
%        d is the vector of x-coordinates only and can be entered either as a column 
%        or row vector, d=[x1,x2,...,xK]
%
%        E-plane gain is evaluated at phi = ph0 for 0 <= theta <= 2*pi. The range 
%        [0,pi] corresponds to the forward ph0-direction and the range [pi,2*pi] to the
%        backward (ph0+pi)-direction. The E-plane gain must be plotted with DBP2 or ABP2.
%
%        H-plane gain is evaluated at theta = pi/2 for 0 <= phi <= 2*pi and must be
%        plotted with DBZ2 or ABZ2.
%
% see also GAIN and GAIN2 (which assumes sinusoidal currents)

% S. J. Orfanidis - 1999 - www.ece.rutgers.edu/~orfanidi/ewa

function [ge,gh,th] = gain2h(L,d,I,N,ph0)

if nargin==0, help gain2h; return; end
if nargin==4, ph0=0; end

K = length(L);

if max(size(d))~=K,
    error('d must have size Kx2 or Kx1 or 1xK');
end
                                                    
if min(size(d))==1,
    d = [d(:),zeros(K,1)];                          % make d into [x,y] pairs
end

ph0 = ph0*pi/180;

th  = 0 : 2*pi/N : 2*pi;                    

for i=1:N+1,                                
    ge(i) = U(L,d,I,th(i),ph0);
    gh(i) = U(L,d,I,pi/2,th(i));                    % here th is the azimuthal angle
end

ge = ge/max(ge);
gh = gh/max(gh);


% ---------------------------------------------------------------

function G = U(L,d,I,th,phi)                        % radiation intensity U(th,phi)

k = 2*pi;   

[N,K] = size(I);
M = (N-1)/2;
Dz = L(:)/N;                                        % K-dim column vector of z-spacings

kx = k*sin(th)*cos(phi);                            % (-th,phi) equivalent to (th,phi+pi)
ky = k*sin(th)*sin(phi);                            % if 0<th<pi, then -th varies over
kz = k*cos(th);                                     % pi<2*pi-th<2*pi

x = d(:,1);                                         % d is entered as [x,y] pairs
y = d(:,2);
z = (-M:M)'*Dz';                                    % NxK matrix of sampled z-points   

F = (exp(j*kz*z).*I) * (exp(j*kx*x) .* exp(j*ky*y) .* Dz .* sinc(Dz*cos(th))); 

Nz = sum(F);                                        % z-component of radiation vector

G = abs(Nz)^2 * sin(th)^2;                          % proportional to radiation intensity


% The factor (exp(j*kz*z).*I) forms the NxK matrix exp(j*kz*z(m,p))*I(m,p), which then
% acts on the column vector of x,y phase and sinc factors, that is, summing over p 
% exp(j*kz*z(m,p))*I(m,p) * exp(j*kx*x(p)+j*ky*y(p))*Dz(p)*sinc(Dz(p)*cos(th)) 
% Then the sum(F) performs the sum over m.

% The sinc factor is essentially flat and arises because of the assumed pulse expansion
% for the currents, that is, I(m) is constant over the interval [z(m)-Dz/2, z(m)+Dz/2]
% MATLAB's sinc is defined as sinc(x) = sin(pi*x)/(pi*x).