Arquivos_Balanis_Matlab

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

Capitulo 15/LEIA-ME.doc
%**********************************************************************
%*******************************************************************
% Reflector (Refletor)
%*******************************************************************
% O programa MATLAB calcula a radiação de um refletor parabólico. 
% O programa pode ser usado para investigar a(o):
%
1. Diretividade
% 
2. Nível de lóbulo secundário 
% 
3. Polarização 
% 
4. Campos nas zonas próxima e distante
% *******************************************************************
% O programa é baseado na formulação e no correspondente programa de computador 
% em Fortran do seguinte artigo:
%
% B. Houshmand, B. Fu, and Y. Rahmat-Samii, “Reflector Antenna Analysis Software”
% Chapter 11, CAEME, Center of Excellence for Multimedia Education, 
% Software Book, Vol. II, Editor by M. F. Iskander, CAEME Center, 
% Electrical Engineering Dept., University of Utah, Salt Lake City, UT.
%
% O leitor é convidado a ler este artigo, para detalhes técnicos. 
% 
%*********************************************************************
% Dados de Entrada Típicos 
 
 Valores 
% ....................................................................................................................................
% Freqüência (GHz) 
3
% ....................................................................................................................................
% Informação sobre o refletor:
% ....................................................................................................................................
% Foco (metro) 
6
% Diâmetro (metro) 
6
% Altura do deslocamento (offset) (metro) 
0
% .....................................................................................................................................
% Informação sobre o alimentador:
% ......................................................................................................................................
% Ponto de alimentação no eixo x (metro) 
0
% Ponto de alimentação no eixo y (metro) 
0
% Ponto de alimentação no eixo z (metro) 
6
% Ângulo de inclinação do alimentador (degrees) 
0
% Número de pontos de integração por comprimento de onda (tipicamente, 4) 
4
% parâmetro q para o diagrama de alimentação [f()=cosq()] 
 
2
% ........................................................................................................................................
% Informação sobre o observador
% .....................................................................
% Distância radial ao ponto de observação (tipicamente, 1.000 m) 
1000
% Mínimo ângulo teta (graus) 
 
0
% Máximo ângulo teta (graus) 
30
% Ângulo fi (graus) 
 
0
% Número de ângulos de observação 
31
% [tipicamente =teta(máx)-teta(mín)+1 para incrementos de 1 grau] 
% .........................................................................................................................................
% Informação sobre a polarização 
% .........................................................................................................................................
% Exemplo: polarização circular
% X POL 
1
% Y POL 
1
% Diferença de fase PSI entre X POL e Y POL (graus) 
90
% ........................................................................................................................................
% ********************************************************************
Capitulo 15/parameters.txt
% TYPICAL INPUT DATA Values 
% ...................................................................
% Frequency (GHz) 3
% ...................................................................
% REFLECTOR INFORMATION:
% ...................................................................
% Focus (meters) 6
% Diameter (meters) 6
% Offset height (meters) 0
% ....................................................................
% FEED INFORMATION:
% ....................................................................
% Feed location of x axis (meters) 0
% Feed location of y axis (meters) 0
% Feed location of z axis (meters) 6
% Feed tilt angle (degrees) 0
% Number of integration points per wavelength (typically 4) 4
% q parameter for feed pattern [f(theta)=cos(theta)^q] 2.3
% .....................................................................
% OBSERVATION INFORMATION
% .....................................................................
% Radial distance to observation point (typically 1,000 m) 1000
% Minimum theta angle (degrees) 0
% Maximum theta angle (degrees) 30
% Phi angle (degrees) 0
% Number of observation angles 31
% [typically =theta(max)-theta(min)+1 for one-degree increments] 
% .....................................................................
% POLARIZATION INFORMATION
% .....................................................................
% X POL 1
% Y POL 1
% Phase difference PSI between X POL and Y POL (degrees) 90
% .....................................................................
Capitulo 15/polar_dB.m
%************************************************************************
% polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
%************************************************************************
%	POLAR_DB is a MATLAB function that plots 2-D patterns in
%	polar coordinates where:
%		 0 <= THETA (in degrees) <= 360
%		-infinity < RHO (in dB) < +infinity
%
%	Input Parameters Description
%	----------------------------
%	- theta (in degrees) must be a row vector from 0 to 360 degrees
%	- rho (in dB) must be a row vector
%	- rmin (in dB) sets the minimum limit of the plot (e.g., -60 dB)
%	- rmax (in dB) sets the maximum limit of the plot (e.g., 0 dB)
%	- rticks is the # of radial ticks (or circles) desired. (e.g., 4)
%	- linestyle is solid (e.g., '-') or dashed (e.g., '--')
%*************************************************************************
%	Credits:
%		S. Bellofiore
%		S. Georgakopoulos
%		A. C. Polycarpou
%		C. Wangsvick
%		C. Bishop
%
%	Tabulate your data accordingly, and call polar_dB to provide the
%	2-D polar plot
%
%	Note: This function is different from the polar.m (provided by
%	 MATLAB) because RHO is given in dB, and it can be negative
%-----------------------------------------------------------------------------
function hpol = polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
% Convert degrees into radians
theta = theta * pi/180;
% Font size, font style and line width parameters
font_size = 16;
font_name = 'Times';
line_width = 1.5;
if nargin < 5
	error('Requires 5 or 6 input arguments.')
elseif nargin == 5 
	if isstr(rho)
		line_style = rho;
		rho = theta;
		[mr,nr] = size(rho);
		if mr == 1
			theta = 1:nr;
		else
			th = (1:mr)';
			theta = th(:,ones(1,nr));
		end
	else
		line_style = 'auto';
	end
elseif nargin == 1
	line_style = 'auto';
	rho = theta;
	[mr,nr] = size(rho);
	if mr == 1
		theta = 1:nr;
	else
		th = (1:mr)';
		theta = th(:,ones(1,nr));
	end
end
if isstr(theta) | isstr(rho)
	error('Input arguments must be numeric.');
end
if any(size(theta) ~= size(rho))
	error('THETA and RHO must be the same size.');
end
% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle = get(cax, 'DefaultTextFontAngle');
fName = get(cax, 'DefaultTextFontName');
fSize = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
set(cax, 'DefaultTextFontAngle', get(cax, 'FontAngle'), ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', font_size, ...
	'DefaultTextFontWeight', get(cax, 'FontWeight') )
% only do grids if hold is off
if ~hold_state
% make a radial grid
	hold on;
 % v returns the axis limits
 % changed the following line to let the y limits become negative
	hhh=plot([0 max(theta(:))],[min(rho(:)) max(rho(:))]);
	v = [get(cax,'xlim') get(cax,'ylim')];
	ticks = length(get(cax,'ytick'));
	delete(hhh);
% check radial limits (rticks)
 	if rticks > 5 % see if we can reduce the number
 		if rem(rticks,2) == 0
 			rticks = rticks/2;
 		elseif rem(rticks,3) == 0
 			rticks = rticks/3;
 		end
 	end
% define a circle
	th = 0:pi/50:2*pi;
	xunit = cos(th);
	yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
 inds = [1:(length(th)-1)/4:length(th)];
 xunits(inds(2:2:4)) = zeros(2,1);
 yunits(inds(1:2:5)) = zeros(3,1);
	rinc = (rmax-rmin)/rticks;
% label r
 % change the following line so that the unit circle is not multiplied
 % by a negative number. Ditto for the text locations.
	for i=(rmin+rinc):rinc:rmax
 is = i - rmin;
		plot(xunit*is,yunit*is,'-','color',tc,'linewidth',0.5);
		text(0,is+rinc/20,[' ' num2str(i)],'verticalalignment','bottom' );
	end
% plot spokes
	th = (1:6)*2*pi/12;
	cst = cos(th); snt = sin(th);
	cs = [-cst; cst];
	sn = [-snt; snt];
	plot((rmax-rmin)*cs,(rmax-rmin)*sn,'-','color',tc,'linewidth',0.5);
% plot the ticks
	george=(rmax-rmin)/30; % Length of the ticks
 th2 = (0:36)*2*pi/72;
 cst2 = cos(th2); snt2 = sin(th2);
	cs2 = [(rmax-rmin-george)*cst2; (rmax-rmin)*cst2];
	sn2 = [(rmax-rmin-george)*snt2; (rmax-rmin)*snt2];
	plot(cs2,sn2,'-','color',tc,'linewidth',0.15); % 0.5
 plot(-cs2,-sn2,'-','color',tc,'linewidth',0.15); % 0.5
% annotate spokes in degrees
 % Changed the next line to make the spokes long enough
	rt = 1.1*(rmax-rmin);
	for i = 1:max(size(th))
		text(rt*cst(i),rt*snt(i),int2str(abs(i*30-90)),'horizontalalignment','center' );
		if i == max(size(th))
			loc = int2str(90);
		elseif i*30+90<=180
			loc = int2str(i*30+90);
 else
 loc = int2str(180-(i*30+90-180)); 
		end
		text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center' );
	end
% set viewto 2-D
	view(0,90);
% set axis limits
 % Changed the next line to scale things properly
	axis((rmax-rmin)*[-1 1 -1.1 1.1]);
end
% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', fSize, ...
	'DefaultTextFontWeight', fWeight );
% transform data to Cartesian coordinates.
 % changed the next line so negative rho are not plotted on the other side
 
 for i = 1:length(rho)
 if (rho(i) > rmin)
 if theta(i)*180/pi >=0 & theta(i)*180/pi <=90
 xx(i) = (rho(i)-rmin)*cos(pi/2-theta(i));
 yy(i) = (rho(i)-rmin)*sin(pi/2-theta(i));
 elseif theta(i)*180/pi >=90
 xx(i) = (rho(i)-rmin)*cos(-theta(i)+pi/2);
 yy(i) = (rho(i)-rmin)*sin(-theta(i)+pi/2);
 elseif theta(i)*180/pi < 0 
 xx(i) = (rho(i)-rmin)*cos(abs(theta(i))+pi/2);
 yy(i) = (rho(i)-rmin)*sin(abs(theta(i))+pi/2);
 end
 else
 xx(i) = 0;
 yy(i) = 0;
 end
 end
% plot data on top of grid
if strcmp(line_style,'auto')
	q = plot(xx,yy);
else
	q = plot(xx,yy,line_style);
end
if nargout > 0
	hpol = q;
end
if ~hold_state
	axis('equal');axis('off');
end
set(q,'linewidth',2);
% reset hold state
if ~hold_state, set(cax,'NextPlot',next); end
Capitulo 15/REFLECTOR.asv
% MATLAB Program Reflector
% The program computes the radiation of a parabolic reflector.
% It can be used to investigate the:
% * Directivity
% * Beamwidth
% * Sidelobe Leve
% * Polarization
% * Near-to-Far Zone Fields
% *****************************************************************
%
% TYPICAL INPUT DATA Values 
% ...................................................................
% Frequency (GHz) 3
% ...................................................................
% REFLECTOR INFORMATION:
% ...................................................................
% Focus (meters) 6
% Diameter (meters) 6
% Offset height (meters) 0
% ....................................................................
% FEED INFORMATION:
% ....................................................................
% Feed location of x axis (meters) 0
% Feed location of y axis (meters) 0
% Feed location of z axis (meters) 6
% Feed tilt angle (degrees) 0
% Number of integration points per wavelength (typically 4) 4
% q parameter for feed pattern [f(theta)=cos(theta)^q] 2.3
% .....................................................................
% OBSERVATION INFORMATION
% .....................................................................
% Radial distance to observation point (typically 1,000 m) 1000
% Minimum theta angle (degrees) 0
% Maximum theta angle (degrees) 30
% Phi angle (degrees) 0
% Number of observation angles 31
% [typically =theta(max)-theta(min)+1 for one-degree increments] 
% .....................................................................
% POLARIZATION INFORMATION
% .....................................................................
% X POL 1
% Y POL 1
% Phase difference ps
%
% Written by: Seunghwan Yoon, Arizona State University, 12/06/2002
%
% *****************************************************************
function[]=reflector;
close all;
clc;
reflant=imread('reflector.bmp');
image(reflant);
axis off;
disp(strvcat('********************************************************************'));
disp(strvcat('Input parameters Values'));
disp(strvcat('********************************************************************'));
freq= input(' Frequency(GHz) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Reflector information'));
disp(strvcat('--------------------------------------------------------------------'));
foc= input(' Focus(meters) =
');
dia= input(' Diameter(meters) = ');
xoset= input(' Offset height(meters) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Feeder information'));
disp(strvcat('--------------------------------------------------------------------'));
xf= input(' Feed location of x axis(meters) = ');
yf= input(' Feed location of y axis(meters) = ');
zf= input(' Feed location of z axis(meters) = ');
fdta= input(' Feed tilt angle(degrees) = ');
ndel= input(' Number of integration points per wavelength(typically 4)= ');
qq= input(' q parameter for feed pattern[f(theta)=cos(theta)^q] = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Observation information '));
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat(' Radial distance to observation'));
rob= input(' (typically 1000 meters for far field) = ');
thstart=input(' Minimum theta angle(degrees) = ');
thetaof=input(' Maximum theta angle(degrees) = ');
phiob= input(' Phi angle(degrees) = ');
disp(strvcat(' Number of observation angle'));
nob= input(' (typically=theta(max)-theta(min)+1 for 1 degree increments) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Polarization information'));
disp(strvcat('--------------------------------------------------------------------'));
apol= input(' X POL = ');
bpol= input(' Y POL = ');
psi= input(' Phase difference PSI between X POL and Y POL(degrees) = ');
disp(strvcat('********************************************************************'));
freq=freq*1e+9;
%degree to radian
rad=pi/180;
fdta=fdta*rad;
thetaof=thetaof*rad;
thstart=thstart*rad;
phiob=phiob*rad;
psi=psi*rad;
c0=3*1e+8; %the speed of electromagnetic wave in free space
rlam=c0/freq; %the wavelength in free space
del=rlam/ndel;
xfac=-1/2/foc;
mx=floor(dia/del);
my=mx;
%input power
qe=qq;
qh=qq;
xxx=sqrt(apol*apol+bpol*bpol);
apol=apol/xxx;
bpol=bpol/xxx;
hwk=exp(j*psi);
hwkn=conj(hwk);
hwk1=exp(j*(psi+pi));
hwkn1=conj(hwk1);
pwrin=(qe+qh+1)/(60*(2*qe+1)*(2*qh+1));
xfactor=sqrt(pwrin*30);
%Loop
h = waitbar(0,'Program is running. Please wait...');
for i=1:nob,
 thetaob=(i-1)*(thetaof-thstart)/(nob-1)+thstart;
 ctob=cos(thetaob);
 stob=sin(thetaob);
 cpob=cos(phiob);
 spob=sin(phiob);
 
 x3=rob*stob*cpob;
 y3=rob*stob*spob;
 z3=rob*ctob;
 
 hsumx=0;
 hsumy=0;
 hsumz=0;
 icount=0;
 
% Integration 
 for jxx=1:mx+1,
 x2=-dia/2+(jxx-1)*del;
 for jyy=1:my+1,
 y2=-dia/2+(jyy-1)*del+xoset;
 if (x2^2+(y2-xoset)^2)<(dia/2)^2
 
 
 z2=((x2^2+y2^2)/4/foc);
 
 y2m=y2+del/2;
 z2m=((x2^2+y2m^2)/4/foc);
 
 yphn=sqrt((del/2)^2+(z2m-z2)^2);
 yph(1)=0;
 yph(2)=del/2/yphn;
 yph(3)=(z2m-z2)/yphn;
 
 nnorm=sqrt(1+(x2^2+y2^2)*(xfac^2));
 norm(1)=x2*xfac/nnorm;
 norm(2)=y2*xfac/nnorm;
 norm(3)=1/nnorm;
 
 [xph]=vector1(yph,norm);
 
 icount=icount+1;
 
 [hx,hy,hz]=incfld(x2,y2,z2,rlam,foc,xf,yf,zf,fdta,qq,apol,bpol,psi);
 
 [jx,jy,jz]=vector(hx,hy,hz,norm(1),norm(2),norm(3));
 
 hj(1)=jx;
 hj(2)=jy;
 hj(3)=jz;
 
 [hxp]=hdot(xph,hj);
 [hyp]=hdot(yph,hj);
 [hzp]=hdot(norm,hj);
 
 p23(1)=x3-x2;
 p23(2)=y3-y2;
 p23(3)=z3-z2;
 
 dist1=sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2);
 
 [pp23(1)]=dot(xph,p23);
 
 [pp23(2)]=dot(yph,p23);
 
 [pp23(3)]=dot(norm,p23);
 if pp23(3)>=dist1
 dist1=pp23(3)+1;
 end
 
 theta2=acos(pp23(3)/(dist1+(1e-8)));
 phi2=atan2(pp23(2),pp23(1)+(1e-8));
 
 [hex,hey,hez]=patchie(dist1,theta2,phi2,hxp,hyp,hzp,del,del,rlam);
 
 hx=hex*xph(1)+hey*yph(1)+hez*norm(1);
 hy=hex*xph(2)+hey*yph(2)+hez*norm(2);
 hz=hex*xph(3)+hey*yph(3)+hez*norm(3);
 
 hsumx=hsumx+hx;
 hsumy=hsumy+hy;
 hsumz=hsumz+hz;
 else
 end
 
 end
 end
 
 
 hsumth=ctob*cpob*hsumx+ctob*spob*hsumy-stob*hsumz;
 
 hsumphi=-spob*hsumx+cpob*hsumy;
 
 hfact11=apol*hwk1*cpob+bpol*spob;
 hfact22=-apol*hwk1*spob+bpol*cpob;
 hfact33=apol*hwkn1*spob-bpol*cpob;
 hfact44=apol*hwkn1*cpob+bpol*spob;
 
 href=hsumth*conj(hfact11)+hsumphi*conj(hfact22);
 hcr=hsumth*conj(hfact33)+hsumphi*conj(hfact44);
 
 href=href*rob/xfactor/(16*pi*pi)*4*pi;
 hcr=hcr*rob/xfactor/(16*pi*pi)*4*pi;
 ex=ctob*cpob*hsumth-spob*hsumphi;
 thetaobi(i)=thetaob;
 hrefi(i)=20*log10(abs(href));
 hcri(i)=20*log10(abs(hcr));
 waitbar(i/nob,h)
 end
close(h)
 
patt = [thetaobi/rad; hrefi; hcri];
fid = fopen('Reflector_Pattern.dat','w');
fprintf(fid,'THETA(degrees) CO-POL(dB) CROSS-POL(dB)\n');
fprintf(fid,'%6.2f %6.3f %6.3f\n',patt);
fclose(fid); 
disp('************************** Program is done **************************');
disp('NOTE:');
disp('THE DIRECTIVITY PATTERN(in dB) IS STORED IN');
disp('AN OUTPUT FILE CALLED ....... Reflector_Pattern.dat');
disp('=====================================================================');
%Plot the graphs
figure;
%subplot(2,2,1)
plot(thetaobi/rad,hrefi,'r','linewidth',2)
title('PARABOLIC REFLECTOR (CO-POL)')
xlabel(['\theta',' (degrees)']),ylabel('DIRECTIVITY(dB)')
grid on;
axis([thstart/rad thetaof/rad 10*floor(min(hrefi)/10) 10*ceil(max(hrefi)/10)]);
figure;
%subplot(2,2,2)
plot(thetaobi/rad,hcri,'b','linewidth',2)
title('PARABOLIC REFLECTOR (CROSS-POL)')
xlabel(['\theta',' (degrees)']),ylabel('DIRECTIVITY(dB)')
grid on;
axis([thstart/rad thetaof/rad 10*floor(min(hcri)/10) 10*ceil(max(hcri)/10)]);
figure;
%subplot(2,2,3)
polar_dB(thetaobi/rad,hrefi,10*floor(min(hrefi)/10),10*ceil(max(hrefi)/10),...
 (10*ceil(max(hrefi)/10)-10*floor(min(hrefi)/10))/10,'r-')
hold on;
polar_dB(-thetaobi/rad,hrefi,10*floor(min(hrefi)/10),10*ceil(max(hrefi)/10),...
 (10*ceil(max(hrefi)/10)-10*floor(min(hrefi)/10))/10,'r-')
title('PARABOLIC REFLECTOR (CO-POL)')
figure;
%subplot(2,2,4)
polar_dB(thetaobi/rad,hcri,10*floor(min(hcri)/10),10*ceil(max(hcri)/10),...
 (10*ceil(max(hcri)/10)-10*floor(min(hcri)/10))/10,'b-')
hold on;
polar_dB(-thetaobi/rad,hcri,10*floor(min(hcri)/10),10*ceil(max(hcri)/10),...
 (10*ceil(max(hcri)/10)-10*floor(min(hcri)/10))/10,'b-')
title('PARABOLIC REFLECTOR (CROSS-POL)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Subroutines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%incfld
function[hx,hy,hz]=incfld(x2,y2,z2,rlam,foc,xf,yf,zf,fdta,qq,apol,bpol,psi)
qe=qq;
qh=qq;
ca=cos(fdta);
sa=sin(fdta);
x2=-x2;
xf=-xf;
xp=x2-xf;
yp=y2-yf;
zp=zf-z2;
xxf=xp;
yyf=ca*yp-sa*zp;
zzf=sa*yp+ca*zp;
r=sqrt(xxf*xxf+yyf*yyf+zzf*zzf);
thetap=acos(zzf/r);
phip=atan2(yyf,xxf+1e-8);
sp=sin(phip);
cp=cos(phip);
st=sin(thetap);
ct=cos(thetap);
hh=exp(-j*2*pi*r/rlam)/r;
hwk=apol*exp(j*psi);
ve=abs(cos(thetap))^qe;
vh=abs(cos(thetap))^qh;
heth=hh*(hwk*cp+bpol*sp)*ve;
heph=hh*(-hwk*sp+bpol*cp)*vh;
hrp=0;
htp=-heph/120/pi;
hpp=heth/120/pi;
hxp=st*cp*hrp+ct*cp*htp-sp*hpp;
hyp=st*sp*hrp+ct*sp*htp+cp*hpp;
hzp=ct*hrp-st*htp;
hhyp=ca*hyp+sa*hzp;
hhzp=-sa*hyp+ca*hzp;
hhxp=-hxp;
hhzp=-hhzp;
hx=hhxp;
hy=hhyp;
hz=hhzp;
%vector
function[jx,jy,jz]=vector(hx,hy,hz,nx,ny,nz)
jx=2*(ny*hz-nz*hy);
jy=2*(nz*hx-nx*hz);
jz=2*(nx*hy-ny*hx);
%vector1
function[hvec3]=vector1(hvec1,hvec2)
hvec3(1)=hvec1(2)*hvec2(3)-hvec1(3)*hvec2(2);
hvec3(2)=hvec1(3)*hvec2(1)-hvec1(1)*hvec2(3);
hvec3(3)=hvec1(1)*hvec2(2)-hvec1(2)*hvec2(1);
%vectorh
function[hvec3]=vectorh(hvec1,hvec2)
hvec3(1)=hvec1(2)*hvec2(3)-hvec1(3)*hvec2(2);
hvec3(2)=hvec1(3)*hvec2(1)-hvec1(1)*hvec2(3);
hvec3(3)=hvec1(1)*hvec2(2)-hvec1(2)*hvec2(1);
%hdot
function[hscalar]=hdot(hvec1,hvec2)
hscalar=hvec1(1)*hvec2(1)+hvec1(2)*hvec2(2)+hvec1(3)*hvec2(3);
%dot
function[hscalar]=dot(hvec1,hvec2)
hscalar=hvec1(1)*hvec2(1)+hvec1(2)*hvec2(2)+hvec1(3)*hvec2(3);
%patch
function[hex,hey,hez]=patchie(dist1,theta2,phi2,hxp,hyp,hzp,xd,yd,rlam)
knum=2*pi/rlam;
zz0=120*pi;
ct=cos(theta2);
st=sin(theta2);
sp=sin(phi2);
cp=cos(phi2);
hh=exp(-j*2*pi*dist1/rlam)/dist1;
hbpol=hyp;
hwk=hxp;
xx=(pi/rlam)*sp*cp*xd+(1e-8);
yy=(pi/rlam)*sp*sp*yd+(1e-8);
fact5=(sin(xx)/xx)*(sin(yy)/yy);
ctcp=ct*cp;
ctsp=ct*sp;
hfact1=sp*hwk-cp*hbpol;
hfact3=ctcp*hwk+ctsp*hbpol;
htp=hfact1*hh*fact5*xd*yd*knum;
hpp=hfact3*hh*fact5*xd*yd*knum;
heth=zz0*hpp;
heph=-zz0*htp;
hex=ct*cp*heth-sp*heph;
hey=ct*sp*heth+cp*heph;
hez=-st*heth;
%-----------------------------------------------------------------------------
% polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
%
%	POLAR_DB is a MATLAB function that plots 2-D patterns in
%	polar coordinates where:
%		 0 <= THETA (in degrees) <= 360
%		-infinity < RHO (in dB) < +infinity
%
%	Input Parameters Description
%	----------------------------
%	- theta (in degrees) must be a row vector from 0 to 360 degrees
%	- rho (in dB) must be a row vector
%	- rmin (in dB) sets the minimum limit of the plot (e.g., -60 dB)
%	- rmax (in dB) sets the maximum limit of the plot (e.g., 0 dB)
%	- rticks is the # of radial ticks (or circles) desired. (e.g., 4)
%	- linestyle is solid (e.g., '-') or dashed (e.g., '--')
%
%	Credits:
%		S. Bellofiore
%		S. Georgakopoulos
%		A. C. Polycarpou
%		C. Wangsvick
%		C. Bishop
%
%	Tabulate your data accordingly, and call polar_dB to provide the
%	2-D polar plot
%
%	Note: This function is different from the polar.m (provided by
%	 MATLAB) because RHO is given in dB, and it can be negative
%-----------------------------------------------------------------------------
function hpol = polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
% Convert degrees into radians
theta = theta * pi/180;
% Font size, font style and line width parameters
font_size = 16;
font_name = 'Times';
line_width = 1.5;
if nargin < 5
	error('Requires 5 or 6 input arguments.')
elseif nargin == 5 
	if isstr(rho)
		line_style = rho;
		rho = theta;
		[mr,nr] = size(rho);
		if mr == 1
			theta = 1:nr;
		else
			th = (1:mr)';
			theta = th(:,ones(1,nr));
		end
	else
		line_style = 'auto';
	end
elseif nargin == 1
	line_style = 'auto';
	rho = theta;
	[mr,nr] = size(rho);
	if mr == 1
		theta = 1:nr;
	else
		th = (1:mr)';
		theta = th(:,ones(1,nr));
	end
end
if isstr(theta) | isstr(rho)
	error('Input arguments must be numeric.');
end
if any(size(theta) ~= size(rho))
	error('THETA and RHO must be the same size.');
end
% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle = get(cax, 'DefaultTextFontAngle');
fName = get(cax, 'DefaultTextFontName');
fSize = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
set(cax, 'DefaultTextFontAngle', get(cax, 'FontAngle'), ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', font_size, ...
	'DefaultTextFontWeight', get(cax, 'FontWeight') )
% only do grids if hold is off
if ~hold_state
% make a radial grid
	hold on;
 % v returns the axis limits
 % changed the following line to let the y limits become negative
	hhh=plot([0 max(theta(:))],[min(rho(:)) max(rho(:))]);
	v = [get(cax,'xlim') get(cax,'ylim')];
	ticks = length(get(cax,'ytick'));
	delete(hhh);
% check radial limits (rticks)
 	if rticks > 5 % see if we can reduce the number
 		if rem(rticks,2) == 0
 			rticks = rticks/2;
 		elseif rem(rticks,3) == 0
 			rticks = rticks/3;
 		end
 	end
% define a circle
	th = 0:pi/50:2*pi;
	xunit = cos(th);
	yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
 inds = [1:(length(th)-1)/4:length(th)];
 xunits(inds(2:2:4)) = zeros(2,1);
 yunits(inds(1:2:5)) = zeros(3,1);
	rinc = (rmax-rmin)/rticks;
% label r
 % change the following line so that the unit circle is not multiplied
 % by a negative number. Ditto for the text locations.
	for i=(rmin+rinc):rinc:rmax
 is = i - rmin;
		plot(xunit*is,yunit*is,'-','color',tc,'linewidth',0.5);
		text(0,is+rinc/20,[' ' num2str(i)],'verticalalignment','bottom' );
	end
% plot spokes
	th = (1:6)*2*pi/12;
	cst = cos(th); snt = sin(th);
	cs = [-cst; cst];
	sn = [-snt; snt];
	plot((rmax-rmin)*cs,(rmax-rmin)*sn,'-','color',tc,'linewidth',0.5);
% plot the ticks
	george=(rmax-rmin)/30; % Length of the ticks
 th2 = (0:36)*2*pi/72;
 cst2 = cos(th2); snt2 = sin(th2);
	cs2 = [(rmax-rmin-george)*cst2; (rmax-rmin)*cst2];
	sn2 = [(rmax-rmin-george)*snt2; (rmax-rmin)*snt2];
	plot(cs2,sn2,'-','color',tc,'linewidth',0.15); % 0.5
 plot(-cs2,-sn2,'-','color',tc,'linewidth',0.15); % 0.5
% annotate spokes in degrees
 % Changed the next line to make the spokes long enough
	rt = 1.1*(rmax-rmin);
	for i = 1:max(size(th))
		text(rt*cst(i),rt*snt(i),int2str(abs(i*30-90)),'horizontalalignment','center' );
		if i == max(size(th))
			loc = int2str(90);
		elseif i*30+90<=180
			loc = int2str(i*30+90);
 else
 loc = int2str(180-(i*30+90-180)); 
		end
		text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center' );
	end
% set viewto 2-D
	view(0,90);
% set axis limits
 % Changed the next line to scale things properly
	axis((rmax-rmin)*[-1 1 -1.1 1.1]);
end
% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', fSize, ...
	'DefaultTextFontWeight', fWeight );
% transform data to Cartesian coordinates.
 % changed the next line so negative rho are not plotted on the other side
 
 for i = 1:length(rho)
 if (rho(i) > rmin)
 if theta(i)*180/pi >=0 & theta(i)*180/pi <=90
 xx(i) = (rho(i)-rmin)*cos(pi/2-theta(i));
 yy(i) = (rho(i)-rmin)*sin(pi/2-theta(i));
 elseif theta(i)*180/pi >=90
 xx(i) = (rho(i)-rmin)*cos(-theta(i)+pi/2);
 yy(i) = (rho(i)-rmin)*sin(-theta(i)+pi/2);
 elseif theta(i)*180/pi < 0 
 xx(i) = (rho(i)-rmin)*cos(abs(theta(i))+pi/2);
 yy(i) = (rho(i)-rmin)*sin(abs(theta(i))+pi/2);
 end
 else
 xx(i) = 0;
 yy(i) = 0;
 end
 end
% plot data on top of grid
if strcmp(line_style,'auto')
	q = plot(xx,yy);
else
	q = plot(xx,yy,line_style);
end
if nargout > 0
	hpol = q;
end
if ~hold_state
	axis('equal');axis('off');
end
% reset hold state
if ~hold_state, set(cax,'NextPlot',next); end
Capitulo 15/reflector.bmp
Capitulo 15/REFLECTOR.m
%**********************************************************************
% Reflector.m
%********************************************************************
% The program computes the radiation by a parabolic reflector.
% It can be used to investigate the:
% * Directivity
% * Beamwidth
% * Sidelobe Level
% * Polarization
% * Near-to-Far Zone Fields
% *****************************************************************
%
% TYPICAL INPUT DATA Values 
% ...................................................................
% Frequency (GHz) 3
% ...................................................................
% REFLECTOR INFORMATION:
% ...................................................................
% Focus (meters) 6
% Diameter (meters) 6
% Offset height (meters) 0
% ....................................................................
% FEED INFORMATION:
% ....................................................................
% Feed location of x axis (meters) 0
% Feed location of y axis (meters) 0
% Feed location of z axis (meters) 6
% Feed tilt angle (degrees) 0
% Number of integration points per wavelength (typically 4) 4
% q parameter for feed pattern [f(theta)=cos(theta)^q] 2.3
% .....................................................................
% OBSERVATION INFORMATION
% .....................................................................
% Radial distance to observation point (typically 1,000 m) 1000
% Minimum theta angle (degrees) 0
% Maximum theta angle (degrees) 30
% Phi angle (degrees) 0
% Number of observation angles 31
% [typically =theta(max)-theta(min)+1 for one-degree increments] 
% .....................................................................
% POLARIZATION INFORMATION
% .....................................................................
% X POL 1
% Y POL 1
% Phase difference PSI between X POL and Y POL (degrees) 90
% .....................................................................
%************************************************************************
% Written by: Seunghwan Yoon, Arizona State University, 12/06/2002
%
%************************************************************************
function[]=reflector;
close all;
clc;
reflant=imread('reflector.bmp');
imagesc(reflant); colormap(gray); 
axis off;
disp(strvcat('********************************************************************'));
disp(strvcat('Input parameters Values'));
disp(strvcat('********************************************************************'));
freq= input(' Frequency(GHz) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Reflector information'));
disp(strvcat('--------------------------------------------------------------------'));
foc= input(' Focus(meters) = ');
dia= input(' Diameter(meters) = ');
xoset= input(' Offset height(meters) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Feeder information'));
disp(strvcat('--------------------------------------------------------------------'));
xf= input(' Feed location of x axis(meters) = ');
yf= input(' Feed location of y axis(meters) = ');
zf= input(' Feed location of z axis(meters) = ');
fdta= input(' Feed tilt angle(degrees) = ');
ndel= input(' Number of integration points per wavelength(typically 4)= ');
qq= input(' q parameter for feed pattern[f(theta)=cos(theta)^q] = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Observation information '));
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat(' Radial distance to observation'));
rob= input(' (typically 1000 meters for far field) = ');
thstart=input(' Minimum theta angle(degrees) = ');
thetaof=input(' Maximum theta angle(degrees) = ');
phiob= input(' Phi angle(degrees) = ');
disp(strvcat(' Number of observation angle'));
nob= input(' (typically=theta(max)-theta(min)+1 for 1 degree increments) = ');
disp(strvcat('--------------------------------------------------------------------'));
disp(strvcat('Polarization information'));
disp(strvcat('--------------------------------------------------------------------'));
apol= input(' X POL = ');
bpol= input(' Y POL = ');
psi= input(' Phase difference PSI between X POL and Y POL(degrees) = ');
disp(strvcat('********************************************************************'));
freq=freq*1e+9;
%degree to radian
rad=pi/180;
fdta=fdta*rad;
thetaof=thetaof*rad;
thstart=thstart*rad;
phiob=phiob*rad;
psi=psi*rad;
c0=3*1e+8; %the speed of electromagnetic wave in free space
rlam=c0/freq; %the wavelength in free space
del=rlam/ndel;
xfac=-1/2/foc;
mx=floor(dia/del);
my=mx;
%input power
qe=qq;
qh=qq;
xxx=sqrt(apol*apol+bpol*bpol);
apol=apol/xxx;
bpol=bpol/xxx;
hwk=exp(j*psi);
hwkn=conj(hwk);
hwk1=exp(j*(psi+pi));
hwkn1=conj(hwk1);
pwrin=(qe+qh+1)/(60*(2*qe+1)*(2*qh+1));
xfactor=sqrt(pwrin*30);
%Loop
h = waitbar(0,'Program is running. Please wait...');
for i=1:nob,
 thetaob=(i-1)*(thetaof-thstart)/(nob-1)+thstart;
 ctob=cos(thetaob);
 stob=sin(thetaob);
 cpob=cos(phiob);
 spob=sin(phiob);
 
 x3=rob*stob*cpob;
 y3=rob*stob*spob;
 z3=rob*ctob;
 
 hsumx=0;
 hsumy=0;
 hsumz=0;
 icount=0;
 
% Integration 
 for jxx=1:mx+1,
 x2=-dia/2+(jxx-1)*del;
 for jyy=1:my+1,
 y2=-dia/2+(jyy-1)*del+xoset;
 if (x2^2+(y2-xoset)^2)<(dia/2)^2
 
 
 z2=((x2^2+y2^2)/4/foc);
 
 y2m=y2+del/2;
 z2m=((x2^2+y2m^2)/4/foc);
 
 yphn=sqrt((del/2)^2+(z2m-z2)^2);
 yph(1)=0;
 yph(2)=del/2/yphn;
 yph(3)=(z2m-z2)/yphn;
 
 nnorm=sqrt(1+(x2^2+y2^2)*(xfac^2));
 norm(1)=x2*xfac/nnorm;
 norm(2)=y2*xfac/nnorm;
 norm(3)=1/nnorm;
 
 [xph]=vector1(yph,norm);
 
 icount=icount+1;
 
 [hx,hy,hz]=incfld(x2,y2,z2,rlam,foc,xf,yf,zf,fdta,qq,apol,bpol,psi);
 
 [jx,jy,jz]=vector(hx,hy,hz,norm(1),norm(2),norm(3));
 
 hj(1)=jx;
 hj(2)=jy;
 hj(3)=jz;
 
 [hxp]=hdot(xph,hj);
 [hyp]=hdot(yph,hj);
 [hzp]=hdot(norm,hj);
 
 p23(1)=x3-x2;
 p23(2)=y3-y2;
 p23(3)=z3-z2;
 
 dist1=sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2);
 
 [pp23(1)]=dot(xph,p23);
 
 [pp23(2)]=dot(yph,p23);
 
 [pp23(3)]=dot(norm,p23);
 if pp23(3)>=dist1
 dist1=pp23(3)+1;
 end
theta2=acos(pp23(3)/(dist1+(1e-8)));
 phi2=atan2(pp23(2),pp23(1)+(1e-8));
 
 [hex,hey,hez]=patchie(dist1,theta2,phi2,hxp,hyp,hzp,del,del,rlam);
 
 hx=hex*xph(1)+hey*yph(1)+hez*norm(1);
 hy=hex*xph(2)+hey*yph(2)+hez*norm(2);
 hz=hex*xph(3)+hey*yph(3)+hez*norm(3);
 
 hsumx=hsumx+hx;
 hsumy=hsumy+hy;
 hsumz=hsumz+hz;
 else
 end
 
 end
 end
 
 
 hsumth=ctob*cpob*hsumx+ctob*spob*hsumy-stob*hsumz;
 
 hsumphi=-spob*hsumx+cpob*hsumy;
 
 hfact11=apol*hwk1*cpob+bpol*spob;
 hfact22=-apol*hwk1*spob+bpol*cpob;
 hfact33=apol*hwkn1*spob-bpol*cpob;
 hfact44=apol*hwkn1*cpob+bpol*spob;
 
 href=hsumth*conj(hfact11)+hsumphi*conj(hfact22);
 hcr=hsumth*conj(hfact33)+hsumphi*conj(hfact44);
 
 href=href*rob/xfactor/(16*pi*pi)*4*pi;
 hcr=hcr*rob/xfactor/(16*pi*pi)*4*pi;
 ex=ctob*cpob*hsumth-spob*hsumphi;
 thetaobi(i)=thetaob;
 hrefi(i)=20*log10(abs(href));
 hcri(i)=20*log10(abs(hcr));
 waitbar(i/nob,h)
 end
close(h)
 
patt = [thetaobi/rad; hrefi; hcri];
fid = fopen('Reflector_Pattern.dat','w');
fprintf(fid,'THETA(degrees) CO-POL(dB) CROSS-POL(dB)\n');
fprintf(fid,'%6.2f %6.3f %6.3f\n',patt);
fclose(fid); 
dmaxdb=max(hrefi);
dmaxdim=10^(dmaxdb/10);
disp('************************** Program is done **************************');
disp('The Maximum Directivity Is:');
disp([num2str(dmaxdim),' (dimensionless); ',num2str(dmaxdb),' (dB)']);
disp(' ');
disp('=====================================================================');
disp('NOTE:');
disp('The Directivity Pattern(in dB) Is Stored In');
disp('An Output File Called ....... Reflector_Pattern.dat');
disp('=====================================================================');
%Plot the graphs
figure;
%subplot(2,2,1)
plot(thetaobi/rad,hrefi,'r','linewidth',2)
title('Parabolic Reflector (CO-POL)')
xlabel(['\theta',' (degrees)']),ylabel('Directivity(dB)')
grid on;
axis([thstart/rad thetaof/rad 10*floor(min(hrefi)/10) 10*ceil(max(hrefi)/10)]);
figure;
%subplot(2,2,2)
plot(thetaobi/rad,hcri,'b','linewidth',2)
title('Parabolic Reflector (CROSS-POL)')
xlabel(['\theta',' (degrees)']),ylabel('Directivity(dB)')
grid on;
axis([thstart/rad thetaof/rad 10*floor(min(hcri)/10) 10*ceil(max(hcri)/10)]);
figure;
%subplot(2,2,3)
polar_dB(thetaobi/rad,hrefi,10*floor(min(hrefi)/10),10*ceil(max(hrefi)/10),...
 (10*ceil(max(hrefi)/10)-10*floor(min(hrefi)/10))/10,'r-')
hold on;
polar_dB(-thetaobi/rad,hrefi,10*floor(min(hrefi)/10),10*ceil(max(hrefi)/10),...
 (10*ceil(max(hrefi)/10)-10*floor(min(hrefi)/10))/10,'r-')
title('Absolute Pattern of Parabolic Reflector (CO-POL)')
% Normalized
figure;
%subplot(2,2,3)
nhrefi=hrefi-max(hrefi);
polar_dB(thetaobi/rad,nhrefi,10*floor(min(nhrefi)/10),10*ceil(max(nhrefi)/10),...
 (10*ceil(max(nhrefi)/10)-10*floor(min(nhrefi)/10))/10,'r-')
hold on;
polar_dB(-thetaobi/rad,nhrefi,10*floor(min(nhrefi)/10),10*ceil(max(nhrefi)/10),...
 (10*ceil(max(nhrefi)/10)-10*floor(min(nhrefi)/10))/10,'r-')
title('Normalized Pattern of Parabolic Reflector (CO-POL)');
figure;
%subplot(2,2,4)
polar_dB(thetaobi/rad,hcri,10*floor(min(hcri)/10),10*ceil(max(hcri)/10),...
 (10*ceil(max(hcri)/10)-10*floor(min(hcri)/10))/10,'b-');
hold on;
polar_dB(-thetaobi/rad,hcri,10*floor(min(hcri)/10),10*ceil(max(hcri)/10),...
 (10*ceil(max(hcri)/10)-10*floor(min(hcri)/10))/10,'b-');
title('Pattern of Parabolic Reflector (CROSS-POL)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Subroutines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%incfld
function[hx,hy,hz]=incfld(x2,y2,z2,rlam,foc,xf,yf,zf,fdta,qq,apol,bpol,psi)
qe=qq;
qh=qq;
ca=cos(fdta);
sa=sin(fdta);
x2=-x2;
xf=-xf;
xp=x2-xf;
yp=y2-yf;
zp=zf-z2;
xxf=xp;
yyf=ca*yp-sa*zp;
zzf=sa*yp+ca*zp;
r=sqrt(xxf*xxf+yyf*yyf+zzf*zzf);
thetap=acos(zzf/r);
phip=atan2(yyf,xxf+1e-8);
sp=sin(phip);
cp=cos(phip);
st=sin(thetap);
ct=cos(thetap);
hh=exp(-j*2*pi*r/rlam)/r;
hwk=apol*exp(j*psi);
ve=abs(cos(thetap))^qe;
vh=abs(cos(thetap))^qh;
heth=hh*(hwk*cp+bpol*sp)*ve;
heph=hh*(-hwk*sp+bpol*cp)*vh;
hrp=0;
htp=-heph/120/pi;
hpp=heth/120/pi;
hxp=st*cp*hrp+ct*cp*htp-sp*hpp;
hyp=st*sp*hrp+ct*sp*htp+cp*hpp;
hzp=ct*hrp-st*htp;
hhyp=ca*hyp+sa*hzp;
hhzp=-sa*hyp+ca*hzp;
hhxp=-hxp;
hhzp=-hhzp;
hx=hhxp;
hy=hhyp;
hz=hhzp;
%vector
function[jx,jy,jz]=vector(hx,hy,hz,nx,ny,nz)
jx=2*(ny*hz-nz*hy);
jy=2*(nz*hx-nx*hz);
jz=2*(nx*hy-ny*hx);
%vector1
function[hvec3]=vector1(hvec1,hvec2)
hvec3(1)=hvec1(2)*hvec2(3)-hvec1(3)*hvec2(2);
hvec3(2)=hvec1(3)*hvec2(1)-hvec1(1)*hvec2(3);
hvec3(3)=hvec1(1)*hvec2(2)-hvec1(2)*hvec2(1);
%vectorh
function[hvec3]=vectorh(hvec1,hvec2)
hvec3(1)=hvec1(2)*hvec2(3)-hvec1(3)*hvec2(2);
hvec3(2)=hvec1(3)*hvec2(1)-hvec1(1)*hvec2(3);
hvec3(3)=hvec1(1)*hvec2(2)-hvec1(2)*hvec2(1);
%hdot
function[hscalar]=hdot(hvec1,hvec2)
hscalar=hvec1(1)*hvec2(1)+hvec1(2)*hvec2(2)+hvec1(3)*hvec2(3);
%dot
function[hscalar]=dot(hvec1,hvec2)
hscalar=hvec1(1)*hvec2(1)+hvec1(2)*hvec2(2)+hvec1(3)*hvec2(3);
%patch
function[hex,hey,hez]=patchie(dist1,theta2,phi2,hxp,hyp,hzp,xd,yd,rlam)
knum=2*pi/rlam;
zz0=120*pi;
ct=cos(theta2);
st=sin(theta2);
sp=sin(phi2);
cp=cos(phi2);
hh=exp(-j*2*pi*dist1/rlam)/dist1;
hbpol=hyp;
hwk=hxp;
xx=(pi/rlam)*sp*cp*xd+(1e-8);
yy=(pi/rlam)*sp*sp*yd+(1e-8);
fact5=(sin(xx)/xx)*(sin(yy)/yy);
ctcp=ct*cp;
ctsp=ct*sp;
hfact1=sp*hwk-cp*hbpol;
hfact3=ctcp*hwk+ctsp*hbpol;
htp=hfact1*hh*fact5*xd*yd*knum;
hpp=hfact3*hh*fact5*xd*yd*knum;
heth=zz0*hpp;
heph=-zz0*htp;
hex=ct*cp*heth-sp*heph;
hey=ct*sp*heth+cp*heph;
hez=-st*heth;
%-----------------------------------------------------------------------------
% polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
%
%	POLAR_DB is a MATLAB function that plots 2-D patterns in
%	polar coordinates where:
%		 0 <= THETA (in degrees) <= 360
%		-infinity < RHO (in dB) < +infinity
%
%	Input Parameters Description
%	----------------------------
%	- theta (in degrees) must be a row vector from 0 to 360 degrees
%	- rho (in dB) must be a row vector
%	- rmin (in dB) sets the minimum limit of the plot (e.g., -60 dB)
%	- rmax (in dB) sets the maximum limit of the plot (e.g., 0 dB)
%	- rticks is the # of radial ticks (or circles) desired. (e.g., 4)
%	- linestyle is solid (e.g., '-') or dashed (e.g., '--')
%
%	Credits:
%		S. Bellofiore
%		S. Georgakopoulos
%		A. C. Polycarpou
%		C. Wangsvick
%		C. Bishop
%
%	Tabulate your data accordingly, and call polar_dB to provide the
%	2-D polar plot
%
%	Note: This function is different from the polar.m (provided by
%	 MATLAB) because RHO is given in dB, and it can be negative
%-----------------------------------------------------------------------------
function hpol = polar_dB(theta,rho,rmin,rmax,rticks,line_style) 
% Convert degrees into radians
theta = theta * pi/180;
% Font size, font style and line width parameters
font_size = 16;
font_name = 'Times';
line_width = 1.5;
if nargin < 5
	error('Requires 5 or 6 input arguments.')
elseif nargin == 5 
	if isstr(rho)
		line_style = rho;
		rho = theta;
		[mr,nr] = size(rho);
		if mr == 1
			theta = 1:nr;
		else
			th = (1:mr)';
			theta = th(:,ones(1,nr));
		end
	else
		line_style = 'auto';
	end
elseif nargin == 1
	line_style = 'auto';
	rho = theta;
	[mr,nr] = size(rho);
	if mr == 1
		theta = 1:nr;
	else
		th = (1:mr)';
		theta = th(:,ones(1,nr));
	end
end
if isstr(theta) | isstr(rho)
	error('Input arguments must be numeric.');
end
if any(size(theta) ~= size(rho))
	error('THETA and RHO must be the same size.');
end
% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle = get(cax, 'DefaultTextFontAngle');
fName = get(cax, 'DefaultTextFontName');
fSize = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
set(cax, 'DefaultTextFontAngle', get(cax, 'FontAngle'), ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', font_size, ...
	'DefaultTextFontWeight', get(cax, 'FontWeight') )
% only do grids if hold is off
if ~hold_state
% make a radial grid
	hold on;
 % v returns the axis limits
 % changed the following line to let the y limits become negative
	hhh=plot([0 max(theta(:))],[min(rho(:)) max(rho(:))]);
	v = [get(cax,'xlim') get(cax,'ylim')];
	ticks = length(get(cax,'ytick'));
	delete(hhh);
% check radial limits (rticks)
 	if rticks > 5 % see if we can reduce the number
 		if rem(rticks,2) == 0
 			rticks = rticks/2;
 		elseif rem(rticks,3) == 0
 			rticks = rticks/3;
 		end
 	end
% define a circle
	th = 0:pi/50:2*pi;
	xunit = cos(th);
	yunit = sin(th);
% now really force points on x/y axes to lie on them exactly
 inds = [1:(length(th)-1)/4:length(th)];
 xunits(inds(2:2:4)) = zeros(2,1);
 yunits(inds(1:2:5)) = zeros(3,1);
	rinc = (rmax-rmin)/rticks;
% label r
 % change the following line so that the unit circle is not multiplied
 % by a negative number. Ditto for the text locations.
	for i=(rmin+rinc):rinc:rmax
 is = i - rmin;
		plot(xunit*is,yunit*is,'-','color',tc,'linewidth',0.5);
		text(0,is+rinc/20,[' ' num2str(i)],'verticalalignment','bottom' );
	end
% plot spokes
	th = (1:6)*2*pi/12;
	cst = cos(th); snt = sin(th);
	cs = [-cst; cst];
	sn = [-snt; snt];
	plot((rmax-rmin)*cs,(rmax-rmin)*sn,'-','color',tc,'linewidth',0.5);
% plot the ticks
	george=(rmax-rmin)/30; % Length of the ticks
 th2 = (0:36)*2*pi/72;
 cst2 = cos(th2); snt2 = sin(th2);
	cs2 = [(rmax-rmin-george)*cst2; (rmax-rmin)*cst2];
	sn2 = [(rmax-rmin-george)*snt2; (rmax-rmin)*snt2];
	plot(cs2,sn2,'-','color',tc,'linewidth',0.15); % 0.5
 plot(-cs2,-sn2,'-','color',tc,'linewidth',0.15); % 0.5
% annotate spokes in degrees
 % Changed the next line to make the spokes long enough
	rt = 1.1*(rmax-rmin);
	for i = 1:max(size(th))
		text(rt*cst(i),rt*snt(i),int2str(abs(i*30-90)),'horizontalalignment','center' );
		if i == max(size(th))
			loc = int2str(90);
		elseif i*30+90<=180
			loc = int2str(i*30+90);
 else
 loc = int2str(180-(i*30+90-180)); 
		end
		text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center' );
	end
% set viewto 2-D
	view(0,90);
% set axis limits
 % Changed the next line to scale things properly
	axis((rmax-rmin)*[-1 1 -1.1 1.1]);
end
% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
	'DefaultTextFontName', font_name, ...
	'DefaultTextFontSize', fSize, ...
	'DefaultTextFontWeight', fWeight );
% transform data to Cartesian coordinates.
 % changed the next line so negative rho are not plotted on the other side
 
 for i = 1:length(rho)
 if (rho(i) > rmin)
 if theta(i)*180/pi >=0 & theta(i)*180/pi <=90
 xx(i) = (rho(i)-rmin)*cos(pi/2-theta(i));
 yy(i) = (rho(i)-rmin)*sin(pi/2-theta(i));
 elseif theta(i)*180/pi >=90
 xx(i) = (rho(i)-rmin)*cos(-theta(i)+pi/2);
 yy(i) = (rho(i)-rmin)*sin(-theta(i)+pi/2);
 elseif theta(i)*180/pi < 0 
 xx(i) = (rho(i)-rmin)*cos(abs(theta(i))+pi/2);
 yy(i) = (rho(i)-rmin)*sin(abs(theta(i))+pi/2);
 end
 else
 xx(i) = 0;
 yy(i) = 0;
 end
 end
% plot data on top of grid
if strcmp(line_style,'auto')
	q = plot(xx,yy);
else
	q = plot(xx,yy,line_style);
end
if nargout > 0
	hpol = q;
end
if ~hold_state
	axis('equal');axis('off');
end
% reset hold state
if ~hold_state, set(cax,'NextPlot',next); end
Capitulo 01/LEIA-ME.doc
Programa de Computador
Animation-Visualization Of Radiation Problems
(Animação- Visualização de Problemas de Radiação)
%***********************************************************
% Este programa MATLAB contém três problemas separados de 
% animação-visualizcao:
%
I.
Fonte-Filamentar; Pulso gaussiano: Meio ilimitado (tm_open)
% 
II.
Fonte-Filamentar; Pulso gaussiano: Cilindro quadrado de 
%
condutor elétrico perfeito (CEP) (tm_box)
% 
III.
Corneta setorial plano E: Meio ilimitado (te_horn)
% O objetivo é permitir que o usuário possa animar e visualizar a radiação 
% em função do tempo para esses três problemas de radiação.
% 
% 
I. 
Fonte-Filamentar; Pulso gaussiano: Meio ilimitado 
%
(tm_open) 
% O primeiro programa de animação-visualização trata de uma fonte 
% filamentar excitada por um único pulso gaussiano que radia em um meio 
% ilimitado; o programa é baseado no método de Diferenças Finitas – 
% Domínio do Tempo (FDTD). O meio ilimitado é simulado usando Condições de Contorno Absorventes (ABC – Absorbing Boundary Conditions), na forma de Camada Perfeitamente Casada (PML
% – Perfectly Matched Layers) de Berenger, para truncar o domínio
% computacional.
% O arquivo-m MATLAB produz a solução FDTD de uma fonte filamentar 
% infinita excitada por um único pulso gaussiano em um domínio 
% computacional bidimensional TMz. O arquivo-m produz um filme com 37 
% quadros de duração; para isso, é tirada uma fotografia do domínio 
% computacional a cada 3 passos temporais. 
%
% 
% II.
Fonte-Filamentar; Pulso gaussiano: Cilindro quadrado 
%
CEP (tm_box)
% O segundo programa de animação-visualização trata de uma fonte 
% filamentar % excitada por um único pulso gaussiano que radia no interior de um cilindro quadrado de condutor elétrico perfeito (CEP). O programa é
% baseado no método de Diferenças Finitas – Domínio do Tempo(FDTD). O % arquivo-m MATLAB produz a solução FDTD de uma fonte filamentar %infinita excitada por um pulso gaussiano em um domínio computacional bidimensional.
%
%
% 
% TMz. 
% O arquivo-m produz um filme com 70 quadros de duração; para isso, é 
% tirada uma fotografia do domínio computacional a cada 3 passos temporais. 
%
% III.
Corneta setorial plano E: meio ilimitado (te_horn)
% O terceiro programa de animação-visualização trata de uma corneta 
% setorial plano E (2-D) que radia em um meio ilimitado; o programa é 
% baseado no método de Diferenças Finitas-Domínio do Tempo. 
% O meio ilimitado é simulado usando Condições de 
% Contorno Absorventes (ABC – Absorbing Boundary Conditions), na forma 
% de Camada Perfeitamente Casada (PML – Perfectly Matched Layers) 
% de Berenger, para truncar o domínio computacional. 
% O arquivo-m MATLAB produz a solução FDTD de uma corneta setorial 
% plano E (2-D) excitada por uma tensão senoidal em um domínio 
% computacional bidimensional TMz. O arquivo-m produz um filme com 37 
% quadros de duração; para isso, é tirada uma fotografia do domínio 
% computacional a cada 3 passos temporais. 
%
% **Nota:
% Para animar e, então, visualizar esses três problemas de radiação, o 
% usuário precisa de uma versão de MATLAB e do arquivo-m MATLAB 
% encontrado no CD anexo para produzir a correspondente solução FDTD 
% de cada problema de radiação. Detalhes adicionais sobre cada um dos 
% problemas de visualização podem ser encontrados no CD que acompanha 
% o livro. 
%************************************************************
Capitulo 01/TE_HORN.m
%***********************************************************************
%
% SECTORAL HORN: UNBOUNDED MEDIUM
%
% PROGRAM AUTHOR-- WILLIAM V. ANDREW
% DEPARTMENT OF ELECTRICAL ENGINEERING
% ARIZONA
STATE UNIVERSITY
% TEMPE, ARIZONA 85287-7206
% (602) 965-5311
%
% DATE OF THIS VERSION-- Nov. 3, 1995
%
% This MATLAB M-file will produce the FDTD solution
% of a sectoral (2-D) Perfectly Electric Conducting
% (PEC) horn antenna excited by a sinusoidal voltage
% in a TEz computational domain. The computational
% domain is truncated with a Berenger Perfectly Matched
% Layer (PML) absorbing boundary condition whose depth
% in layers is set by the variable NPMLS. The PML is
% introduced to eliminate reflections from the grid
% truncation and to simulate an outgoing traveling
% wave propagating in an unbounded medium. The M-file
% can also create a movie. For example, you can create
% a movie which is 70 frames long by taking a picture of 
% the computational domain every 3rd time step.
% 
% To execute this M-file, type ``te_horn'' at the
% MATLAB prompt. The file will save each frame of
% the movie and then write the entire movie to a file
% named ``te_horn.mat''.
%
% To play the movie at any time after it has been 
% created and saved to file ``te_horn.mat'' just
% execute the following MATLAB commands:
%
% load te_horn.mat
% movie(M,n,fps) 
%
% where n is the number of times to play the movie
% and fps is the number of frames per second. 
%
% This M-file will not work with the Student edition
% of MATLAB due to the restrictions on array size.
% Therefore, this M-file will work only with the
% Professional edition of MATLAB. The movie which
% this file creates is approximately 10.6 Mbytes in
% size. Therefore, the available RAM memory or the
% swap space on whatever operating system must be
% large enough to accommodate a file this large.
%
% The horn is modeled by setting the necessary
% FDTD update equation coefficients to represent
% PEC material (sigma=infinity). 
% The cell size of the space is:
% dx = 0.0025 meters
% The time step is:
% dt = 4.23e-12 seconds
% The frequency of excitation is:
% freq = 9.84252 GHz
% The wavelength is:
% lambda = 12*dx = 0.0305 meters
%
% The flare section of the horn is staircased. As
% modeled, the horn looks like:
%
% jc-7 jc+7 
% | |
% \ | / \ | / 
% \|/ \|/
% ` ` 
%
% ic+11 |_ _| 
% ic+10 | |
% ic+9 |_ _|
% ic+8 | |
% ic+7 |_ _|
% ic+6 | |
% ic+5 |_ _|
% ic+4 | |
% ic+3 ------------> |_ _|
% ic+2 | |
% ic+1 | |
% ic | |
% ic-1 | |
% ic-2 | |
% ic-3 | |
% ic-4 | | `|' <--- Ex Field Component
% ic-5 | |
% ic-6 | | `-' <--- Ey Field Component
% ic-7 | |
% ic-8 | |
% ic-9 | | 
% ic-10--------------> |_ _ _ _| <----- Excitation Plane
% ic-11 | | 
% ic-12 | |
% ic-13---------------> |_ _ _ _|
%
% . .
% /|\ /|\
% / | \ / | \
% | |
% jc-2 jc+2
%
%***********************************************************************
clear 
REPEAT = input('How many times to repeat? ');
FPS = input('Enter frames per second... ');
%***********************************************************************
% Initialize some constants
%***********************************************************************
 npmls=8; % Depth of PML region in # of cells
 
 nmax=210; % Number of time steps
 ie=100; 
 ib=ie+1;
 ic=ie/2-20;
 ip=ie-npmls;
 
 je=100;
 jb=je+1;
 jc=je/2;
 jp=je-npmls;
 pi=4.0*atan(1.0);
 muo=4.0*pi*1.0e-7; % Permeability of free space
 epso=8.854e-12; % Permittivity of free space
 co=1.0/sqrt(muo*epso); % Speed of light in free space
 aimp=sqrt(muo/epso); % Wave impedance in free space
 freq=9.84252e+09; % Frequency of excitation
 lambda=co/freq; % Wavelength of excitation
 dx=lambda/12.0; % FDTD cell size
 dt=dx/co/2.0; % Time step size
%***********************************************************************
% .... Set up the Berenger PML ABC material constants ....
%***********************************************************************
 sigmax=-3.0*epso*co*log(1e-5)/(2.0*dx*npmls);
 rhomax=sigmax*(aimp^2);
 for m=1:npmls;
 sig(m)=sigmax*((m-0.5)/(npmls+0.5))^2;
 rho(m)=rhomax*(m/(npmls+0.5))^2;
 end;
%***********************************************************************
% .... Set up constants needed in the FDTD equations for the ....
% .... Berenger PML ABCs (exponential difference expressions)....
%***********************************************************************
 for m=1:npmls;
 re=sig(m)*dt/epso;
 rm=rho(m)*dt/muo;
 ca(m)=exp(-re);
 cb(m)=-(exp(-re)-1.0)/sig(m)/dx;
 da(m)=exp(-rm);
 db(m)=-(exp(-rm)-1.0)/rho(m)/dx;
 end;
%***********************************************************************
% Initialize all of the matrices for the field components HZ, HZX,
% HZY, EX, EY, CAEX, CAEY, DAHZX, DAHZY, CBEX, CBEY, DBHZX, and DBHZY.
%***********************************************************************
 for i=1:ie;
 for j=1:jb;
	ex(i,j)=0.0;
	caex(i,j)=1.0; % Free space
	cbex(i,j)=dt/epso/dx; % Free space
 end;
 end;
 for i=1:ib;
 for j=1:je;
	ey(i,j)=0.0;
	caey(i,j)=1.0; % Free space
	cbey(i,j)=dt/epso/dx; % Free space
 end;
 end;
 for i=1:ie;
 for j=1:je;
	hz(i,j)=0.0;
	hzx(i,j)=0.0;
	dahzx(i,j)=1.0; % Free space
	dbhzx(i,j)=dt/muo/dx; % Free space
	hzy(i,j)=0.0;
	dahzy(i,j)=1.0; % Free space
	dbhzy(i,j)=dt/muo/dx; % Free space
 end;
 end;
%*******************************************************************
% Initialize all of the matrices for the Berenger PML absorbing
% boundaries.
%*******************************************************************
%<<<<<<<<<<<<<<<<<<<<<<<<< Ex Fields >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% ..... Left and Right PML Regions .....
 for i=2:ie;
 for j=2:npmls+1;
	m=npmls+2-j;
	caex(i,j)=ca(m);
	cbex(i,j)=cb(m);
 end;
 for j=jp+1:je;
	m=j-jp;
	caex(i,j)=ca(m);
	cbex(i,j)=cb(m);
 end;
 end;
%<<<<<<<<<<<<<<<<<<<<<<<<< Ey Fields >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% ..... Back and Front PML Regions .....
 for j=2:je;
 for i=2:npmls+1;
	m=npmls+2-i;
	caey(i,j)=ca(m);
	cbey(i,j)=cb(m);
 end;
 for i=ip+1:ie;
	m=i-ip;
	caey(i,j)=ca(m);
	cbey(i,j)=cb(m);
 end;
 end;
%<<<<<<<<<<<<<<<<<<<<<<<<< Hz Fields >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% ..... Left and Right PML Regions .....
 for i=2:ie;
 for j=1:npmls;
	m=npmls+1-j;
	dahzy(i,j)=da(m);
	dbhzy(i,j)=db(m);
 end;
for j=jp+1:je;
	m=j-jp;
	dahzy(i,j)=da(m);
	dbhzy(i,j)=db(m);
 end;
 end;
% ..... Front and Back PML Regions .....
 for j=2:je;
 for i=1:npmls;
	m=npmls+1-i;
	dahzx(i,j)=da(m);
	dbhzx(i,j)=db(m);
 end;
 for i=ip+1:ie;
	m=i-ip;
	dahzx(i,j)=da(m);
	dbhzx(i,j)=db(m);
 end;
 end;
%***********************************************************************
% ..... Material coefficients PEC of horn antenna .....
%***********************************************************************
 for i=ic-13:ic+2;
 caex(i,jc-2)=-1.0;
 cbex(i,jc-2)= 0.0;
 caex(i,jc+2)=-1.0;
 cbex(i,jc+2)= 0.0;
 end;
 caex(ic+3,jc-3)=-1.0;
 cbex(ic+3,jc-3)= 0.0;
 caex(ic+3,jc+3)=-1.0;
 cbex(ic+3,jc+3)= 0.0;
 caex(ic+4,jc-3)=-1.0;
 cbex(ic+4,jc-3)= 0.0;
 caex(ic+4,jc+3)=-1.0;
 cbex(ic+4,jc+3)= 0.0;
 caex(ic+5,jc-4)=-1.0;
 cbex(ic+5,jc-4)= 0.0;
 caex(ic+5,jc+4)=-1.0;
 cbex(ic+5,jc+4)= 0.0;
 caex(ic+6,jc-4)=-1.0;
 cbex(ic+6,jc-4)= 0.0;
 caex(ic+6,jc+4)=-1.0;
 cbex(ic+6,jc+4)= 0.0;
 caex(ic+7,jc-5)=-1.0;
 cbex(ic+7,jc-5)= 0.0;
 caex(ic+7,jc+5)=-1.0;
 cbex(ic+7,jc+5)= 0.0;
 caex(ic+8,jc-5)=-1.0;
 cbex(ic+8,jc-5)= 0.0;
 caex(ic+8,jc+5)=-1.0;
 cbex(ic+8,jc+5)= 0.0;
 caex(ic+9,jc-6)=-1.0;
 cbex(ic+9,jc-6)= 0.0;
 caex(ic+9,jc+6)=-1.0;
 cbex(ic+9,jc+6)= 0.0;
 caex(ic+10,jc-6)=-1.0;
 cbex(ic+10,jc-6)= 0.0;
 caex(ic+10,jc+6)=-1.0;
 cbex(ic+10,jc+6)= 0.0;
 caex(ic+11,jc-7)=-1.0;
 cbex(ic+11,jc-7)= 0.0;
 caex(ic+11,jc+7)=-1.0;
 cbex(ic+11,jc+7)= 0.0;
 for j=jc-2:jc+1;
 caey(ic-13,j)=-1.0;
 cbey(ic-13,j)= 0.0;
 end;
 caey(ic+3,jc-3)=-1.0;
 cbey(ic+3,jc-3)= 0.0;
 caey(ic+3,jc+2)=-1.0;
 cbey(ic+3,jc+2)= 0.0;
 caey(ic+5,jc-4)=-1.0;
 cbey(ic+5,jc-4)= 0.0;
 caey(ic+5,jc+3)=-1.0;
 cbey(ic+5,jc+3)= 0.0;
 caey(ic+7,jc-5)=-1.0;
 cbey(ic+7,jc-5)= 0.0;
 caey(ic+7,jc+4)=-1.0;
 cbey(ic+7,jc+4)= 0.0;
 caey(ic+9,jc-6)=-1.0;
 cbey(ic+9,jc-6)= 0.0;
 caey(ic+9,jc+5)=-1.0;
 cbey(ic+9,jc+5)= 0.0;
 caey(ic+11,jc-7)=-1.0;
 cbey(ic+11,jc-7)= 0.0;
 caey(ic+11,jc+6)=-1.0;
 cbey(ic+11,jc+6)= 0.0;
%***********************************************************************
% .....TIME-STEPPING LOOP.....
%***********************************************************************
 for n=1:nmax;
%***********************************************************************
% .....EX FIELD UPDATE.....
%***********************************************************************
 ex(1:ie,2:je)=caex(1:ie,2:je).*ex(1:ie,2:je)+...
	cbex(1:ie,2:je).*(hz(1:ie,2:je)-hz(1:ie,1:je-1));
%***********************************************************************
% .....EY FIELD UPDATE.....
%***********************************************************************
 ey(2:ie,1:je)=caey(2:ie,1:je).*ey(2:ie,1:je)+...
 cbey(2:ie,1:je).*(hz(1:ie-1,1:je)-hz(2:ie,1:je));
 
%***********************************************************************
% ..... Hard Source ramped sinusoidal excitation .....
%***********************************************************************
 for j=jc-2:jc+1;
 ey(ic-10,j)=(1.0-exp(-((n/20.0)^2)))*...
 aimp*sin(2.0*pi*freq*n*dt);
 end;
%***********************************************************************
% .....HZ FIELD UPDATE..... 
%***********************************************************************
 hzx(1:ie,1:je)=dahzx(1:ie,1:je).*hzx(1:ie,1:je)+...
 dbhzx(1:ie,1:je).*(ey(1:ie,1:je)-ey(2:ib,1:je));
 hzy(1:ie,1:je)=dahzy(1:ie,1:je).*hzy(1:ie,1:je)+...
 dbhzy(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je));
 hz(1:ie,1:je)=hzx(1:ie,1:je)+hzy(1:ie,1:je);
%***********************************************************************
% .....Create the movie frame by frame.....
% .....Take a frame every 3rd time step.....
%***********************************************************************
 if rem(n,3)==0;
 s=int2str(n);
 n2=n/3;
 clf;
 pcolor(log10(abs(ey+0.000001)));
 axis([0 100 0 100]);
 caxis([-6 3]);
 shading interp;
 if n==3;
	M=moviein(70);
 end;
 t2=['TEz 2D Horn Antenna. Time step #',s];
 title(t2);
 hold;
 M(:,n2)=getframe;
 end;
%***********************************************************************
% End time step loop
%***********************************************************************
 end;
%***********************************************************************
% Save the movie to file ``te_horn.mat''
%***********************************************************************
 save te_horn M;
%***********************************************************************
% Replay the movie 5 times at 7 frames per second
%***********************************************************************
 movie(M,REPEAT,FPS)
�
Capitulo 01/TM_BOX.m
%***********************************************************************
%
% LINE SOURCE EXCITED WITH A TIME DERIVATIVE GAUSSIAN PULSE:
% PERFECTLY ELECTRIC CONDUCTING (PEC) CYLINDER
%
% PROGRAM AUTHOR-- WILLIAM V. ANDREW
% DEPARTMENT OF ELECTRICAL ENGINEERING
% ARIZONA STATE UNIVERSITY
% TEMPE, ARIZONA 85287-7206
% (602) 965-5311
%
% DATE OF THIS VERSION-- Nov. 3, 1995
%
% This MATLAB M-file will produce the FDTD solution of a z-directed
% line source in a two-dimensional TMz computational domain excited
% by a time derivative Gaussian pulse. The computational domain is
% truncated with a PEC material. The walls of the PEC cylinder
% introduce reflections which interfere with the outgoing waves to
% create standing waves. The M-file
% can also create a movie. For example, you can create
% a movie which is 70 frames long by taking a picture of 
% the computational domain every 3rd time step.
% 
% To execute the M-file, type ``tm_box'' at the MATLAB prompt.
% The file will save each frame of the movie, replay the movie 5
% times at 7 frames per second and then write the entire movie
% to a file named ``tm_box.mat''.
%
% To play the movie at any time after it has been created and
% saved to the file ``tm_box.mat'' just execute the following
% MATLAB commands:
%
% load tm_box.mat
% movie(M,n,fps) 
%
% where n is the number of times to play the movie and 
% fps is the number of frames per second. 
%
% This M-file will not work with the Student edition of MATLAB
% due to the restrictions on array size. Therefore, this M-file
% will work only with the Professional edition of MATLAB. The
% movie which this file creates is approximately 10.6 Mbytes in size.
% Therefore, the available RAM memory or the swap space on whatever
% operating system must be large enough to accommodate a file 
% this large.
%
%***********************************************************************
clear 
REPEAT = input('How many times to repeat? ');
FPS = input('Enter frames per second... ');
%***********************************************************************
% Initialize some constants
%***********************************************************************
 nmax=210; % Number of time steps
 ie=50; 
 ib=ie+1;
 ic=ie/2+1;
je=50;
 jb=je+1;
 jc=je/2+1;
 pi=4.0*atan(1.0);
 muo=4.0*pi*1.0e-7; % Permeability of free space
 epso=8.854e-12; % Permittivity of free space
 co=1.0/sqrt(muo*epso); % Speed of light in free space
 aimp=sqrt(muo/epso); % Wave impedance in free space
 dx=0.003; % FDTD cell size
 dt=dx/co/2.0; % Time step size
%***********************************************************************
% Initialize all of the matrices for the field components HX, HY,
% EZX, EZY, DAHX, DAHY, CAEZX, CAEZY, DBHX, DBHY, CBEZX, and CBEZY.
%***********************************************************************
 for i=1:ib;
 for j=1:jb;
	ez(i,j)=0.0;
	caez(i,j)=1.0; % Free space
	cbez(i,j)=dt/epso/dx; % Free space
 end;
 end;
 for i=1:ib;
 for j=1:je;
	hx(i,j)=0.0;
	dahx(i,j)=1.0; % Free space
	dbhx(i,j)=dt/muo/dx; % Free space
 end;
 end;
 for i=1:ie;
 for j=1:jb;
	hy(i,j)=0.0;
	dahy(i,j)=1.0; % Free space
	dbhy(i,j)=dt/muo/dx; % Free space
 end;
 end;
%***********************************************************************
% .....TIME-STEPPING LOOP.....
%***********************************************************************
 for n=1:nmax;
%***********************************************************************
% .....Set up time derivative Gaussian pulse excitation voltage.....
%***********************************************************************
 b = 25.0;
 dum = 4.0/b/dt*(n*dt-b*dt);
 voltage = 2.0*dum*exp(-(dum^2));
%***********************************************************************
% .....EZ FIELD UPDATE.....
%***********************************************************************
 ez(2:ie,2:je)=caez(2:ie,2:je).*ez(2:ie,2:je)+cbez(2:ie,2:je).*...
 (hy(2:ie,2:je)-hy(1:ie-1,2:je)+hx(2:ie,1:je-1)-hx(2:ie,2:je));
%***********************************************************************
% ..... Hard Source excitation .....
%***********************************************************************
 ez(ic,jc)=voltage/dx;
%***********************************************************************
% .....HX FIELD UPDATE..... 
%***********************************************************************
 hx(1:ib,1:je)=dahx(1:ib,1:je).*hx(1:ib,1:je)+...
 dbhx(1:ib,1:je).*(ez(1:ib,1:je)-ez(1:ib,2:jb));
%***********************************************************************
% .....HY FIELD UPDATE.....
%***********************************************************************
 hy(1:ie,1:jb)=dahy(1:ie,1:jb).*hy(1:ie,1:jb)+...
 dbhy(1:ie,1:jb).*(ez(2:ib,1:jb)-ez(1:ie,1:jb));
%***********************************************************************
% .....Create the movie frame by frame.....
% .....Take a frame every 3rd time step.....
%***********************************************************************
 if rem(n,3)==0;
 s=int2str(n);
 n2=n/3;
 clf;
 pcolor(ez);
 axis([0 50 0 50]);
 caxis([-50 50]);
 shading interp;
 if n==3;
	M=moviein(37);
 end;
 t2=['TMz Gaussian Derivative Pulse. Time step #',s];
 title(t2);
 hold;
 M(:,n2)=getframe;
 end;
%***********************************************************************
% End time step loop
%***********************************************************************
 end;
%***********************************************************************
% Save the movie to file ``tm_box.mat''
%***********************************************************************
 save tm_box M;
%***********************************************************************
% Replay the movie 5 times at 7 frames per second
%***********************************************************************
 movie(M,REPEAT,FPS)
�
Capitulo 01/TM_OPEN.m
%***********************************************************************
%
% LINE SOURCE EXCITED WITH A TIME DERIVATIVE GAUSSIAN PULSE:
% UNBOUNDED MEDIUM
%
% PROGRAM AUTHOR-- WILLIAM V. ANDREW
% DEPARTMENT OF ELECTRICAL ENGINEERING
% ARIZONA STATE UNIVERSITY
% TEMPE, ARIZONA 85287-7206
% (602) 965-5311
%
% DATE OF THIS VERSION-- Nov. 3, 1995
%
% This MATLAB M-file will produce the FDTD solution of a z-directed
% line source in a two-dimensional TMz computational domain excited
% by a time derivative Gaussian pulse. The computational domain
% is truncated with a Berenger Perfectly Matched Layer (PML)
% absorbing boundary condition whose depth in layers is set by the
% variable NPMLS. The PML is introduced to eliminate reflections
% from the grid truncation and to simulate an outgoing traveling
% wave propagating in an unbounded medium. The M-file
% can also create a movie. For example, you can create
% a movie which is 70 frames long by taking a picture of 
% the computational domain every 3rd time step.
% 
% To execute the M-file, type ``tm_open'' at the MATLAB prompt.
% The file will save each frame of the movie, replay the movie 5
% times at 7 frames per second and then write the entire movie
% to a file named ``tm_open.mat''.
%
% To play the movie at any time after it has been created and
% saved to the file ``tm_open.mat'' just execute the following
% MATLAB commands:
%
% load tm_open.mat
% movie(M,n,fps) 
%
% where n is the number of times to play the movie and 
% fps is the number of frames per second. 
%
% This M-file will not work with the Student edition of MATLAB
% due to the restrictions on array size. Therefore, this M-file
% will work only with the Professional edition of MATLAB. The
% movie which this file creates is approximately 5.6 Mbytes in size.
% Therefore, the available RAM memory or the swap space on whatever
% operating system must be large enough to accommodate a file 
% this large.
%
%***********************************************************************
clear 
REPEAT = input('How many times to repeat? ');
FPS = input('Enter frames per second... ');
%***********************************************************************
% Initialize some constants
%***********************************************************************
 npmls=6; % Depth of PML region in # of cells
 
 nmax=111; % Number of time steps
 ie=50; 
 ib=ie+1;
 ic=ie/2+1;
 ip=ie-npmls;
 
 je=50;
 jb=je+1;
 jc=je/2+1;
 jp=je-npmls;
 pi=4.0*atan(1.0);
 muo=4.0*pi*1.0e-7; % Permeability of free space
 epso=8.854e-12; % Permittivity of free space
 co=1.0/sqrt(muo*epso); % Speed of light in free space
 aimp=sqrt(muo/epso); % Wave impedance in free space
% freq=9.84252e+09;
% lambda=co/freq;
 dx=0.003; % FDTD cell size
 dt=dx/co/2.0; % Time step size
%***********************************************************************
% .... Set up the Berenger PML ABC material constants ....
%***********************************************************************
 sigmax=-3.0*epso*co*log(1e-5)/(2.0*dx*npmls);
 rhomax=sigmax*(aimp^2);
 for m=1:npmls;
 sig(m)=sigmax*((m-0.5)/(npmls+0.5))^2;
 rho(m)=rhomax*(m/(npmls+0.5))^2;
 end;
%***********************************************************************
% .... Set up constants needed in the FDTD equations