MATLAB绘制Smith圆图及反射系数

%Victor Aprea Cornell University 6/27/02
%
%Usage:  plotsmithchart(Zl,Zo)
%  where Zl is the Load Impedence (possibly complex)
%  and Zo is the Characteristic Line Impedence 
%  Plots a smith chart, along with the reflection coefficient circle
%  and the line of intersection with resistive component equal to 1.

%  plotsmithchart
%  Without any parameters draws a blank smith chart.
%  Wavelengths toward the generator are labeled around the perimeter
%
%For example: plotsmithchart(25,50)
%  Draws a smithchart, calculates and plots the reflection coefficient
%  for a characteristic impedence of 50 ohms and a load impedence of 25 ohms,
%  and draws the line of intersection with the R=1 circle.

function answer = plotSmithChart(Zl,Zo);
constant = linspace(0,10,5);
phaseAngle = linspace(0,2*pi,50);
unitGamma = exp(j*phaseAngle);
%plot the unit circle in the complex plane
hold on;
plot(real(unitGamma),imag(unitGamma),'r');
%set(gcf,'Position',[0 0 1280 990]);
axis square
zoom on
axis([-1.1 1.1 -1.1 1.1]);

MAX=2001;

bound2=0;
bound3=0;
min_bound1=0;
min_bound2=0;
max_bound2=0;
H=0;
word=0;

Gr = linspace(-1,1,MAX);
hold on;

interval = [[.01:.01:.2],[.22:.02:.5],[.55:.05:1],[1.1:.1:2],[2.2:.2:5],[6:10],[12:2:20],[30:10:50]];
interval2= [[.01:.01:.5],[.55:.05:1],[1.1:.1:2],[2.2:.2:5],[6:10],[12:2:20],[30:10:50]];

%plot real axis
plot(Gr, zeros(1,length(Gr)),'r');

%equations were derived using the symbolic toolbox as follows
%solve('R=(1-Gr^2-Gi^2)/((1-Gr)^2+Gi^2)','Gi')
%bound was derived as follows
%solve('1/(R+1)*(-(R+1)*(R-2*R*Gr+R.*Gr^2-1+Gr^2))^(1/2)=0','Gr')
for R = interval2,
   min_bound1 = (R-1)/(R+1);
   
   if(R<.2)        
      if(mod(R,.1)==0) 
         max_bound = (-1+2^2+R^2)/(2^2+R^2+2*R+1);
      elseif(mod(R,.02)==0)
         max_bound = (-1+.5^2+R^2)/(.5^2+R^2+2*R+1);
      else