史密斯圆图的matlab程序

%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

%Victor Aprea Cornell University 6/27/02
%
%Usage:		plotsmithchart（ZlZo）
%		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（2550）
%		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（ZlZo）;
constant = linspace（0105）;
phaseAngle = linspace（02*pi50）;
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（-11MAX）;
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（1length（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
         max_bound = （-1+.2^2+R^2）/（.2^2+R^2+2*R+1）;
         if（R==.05 | （R<.151 & R>.149））
            min_bound2 = （-1+.5^2+R^2）/（.5^2+R^2+2*R+1）;
            max_bound2 = （-1+1^2+R^2）/（1^2+R^2+2*R+1）;
         end
      end	
   elseif（R<1）
      if（mod（R.2）==0） 
         max_bound = （-1+5^2+R^2）/（5^2+R^2+2*R+1）;
      elseif（mod（R.1）==0） 
         max_bound = （-1+2^2+R^2）/（2^2+R^2+2*R+1）; 
      elseif（R==.25 | R==.35 | R==.45）
         temp = （-1+.5^2+R^2）/（.5^2+R^2+2*R+1）;
         min_bound2 = max（min_bound1 temp）;
         max_bound = （-1+1^2+R^2）/（1^2+R^2+2*R+1）;
      elseif（R<.5） 
         max_bound = （-1+.5^2+R^2）/（.5^2+R^2+2*R+1）; 
      else 
         max_bound = （-1+1^2+R^2）/（1^2+R^2+2*R+1）; 
      end      
   elseif（R<5）
      if（mod（R2）==0）
         max_bound = （-1+20^2+R^2）/（20^2+R^2+2*R+1）; 
      elseif（mod（R1）==0）
         max_bound = （-1+10^2+R^2）/（10^2+R^2+2*R+1）; 
      elseif（R>2）
         max_bound = （-1+5^2+R^2）/（5^2+R^2+2*R+1）;
      else
         if（mod（R.2）==0）
	         max_bound = （-1+5^2+R^2）/（5^2+R^2+2*R+1）;
         else