plotsmithchart.m

使用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
         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(R,2)==0)
         max_bound = (-1+20^2+R^2)/(20^2+R^2+2*R+1); 
      elseif(mod(R,1)==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
            max_bound = (-1+2^2+R^2)/(2^2+R^2+2*R+1);       
         end
      end
   elseif(R<10)
      if(mod(R,2)==0)
         max_bound = (-1+20^2+R^2)/(20^2+R^2+2*R+1); 
      else
         max_bound = (-1+10^2+R^2)/(10^2+R^2+2*R+1); 
      end
   else
      if(R==10|R==20)
         max_bound = (-1+50^2+R^2)/(50^2+R^2+2*R+1); 
      elseif(R==50)
         max_bound = 1;
      elseif(R<20)
         max_bound = (-1+20^2+R^2)/(20^2+R^2+2*R+1); 
      else
         max_bound = (-1+50^2+R^2)/(50^2+R^2+2*R+1); 
      end
   end      
   
   index = ceil((min_bound1+1)*(MAX-1)/2+1);
   actual_value = Gr(index);
   if(actual_value<min_bound1)
      index = index + 1;
   end
   MIN=index;  
   
   index = ceil((max_bound+1)*(MAX-1)/2+1);
   actual_value = Gr(index);
   if(actual_value>max_bound)
      index = index - 1;
   end
   
   MIN2 = ceil((min_bound2+1)*(MAX-1)/2+1);
   actual_value = Gr(MIN2);
   if(actual_value<min_bound2	)
      MIN2 = MIN2 + 1;
   end 
   
   MAX2 = ceil((max_bound2+1)*(MAX-1)/2+1);
   actual_value = Gr(MAX2);
   if(actual_value<max_bound2	)
      MAX2 = MAX2 + 1;
	end	

   r_L_a=1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN:index)+R.*Gr(MIN:index).^2-1+Gr(MIN:index).^2)).^(1/2);
   r_L_b=-1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN:index)+R.*Gr(MIN:index).^2-1+Gr(MIN:index).^2)).^(1/2);
   r_L_b(1)=0;
   r_L_a(1)=0;
      
   r_L_a2=1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN2:MAX2)+R.*Gr(MIN2:MAX2).^2-1+Gr(MIN2:MAX2).^2)).^(1/2);
   r_L_b2=-1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN2:MAX2)+R.*Gr(MIN2:MAX2).^2-1+Gr(MIN2:MAX2).^2)).^(1/2);
   
   r_L_a3=1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN2:index)+R.*Gr(MIN2:index).^2-1+Gr(MIN2:index).^2)).^(1/2);
   r_L_b3=-1/(R+1)*(-(R+1)*(R-2*R.*Gr(MIN2:index)+R.*Gr(MIN2:index).^2-1+Gr(MIN2:index).^2)).^(1/2);
   
   %fix resolution issues in .2-.5 range
   
   if(~(R>.2 & R<.5 & ~(mod(R,.02)==0)))
		if(R==1)
         color = 'r';
      else
         color ='b';
      end            
      plot(Gr(MIN:index),r_L_a(1:index-MIN+1),color,Gr(MIN:index), r_L_b(1:index-MIN+1),color);   
      if(R<=1)
      	if(mod(R,1)==0)
         	word = [num2str(R) '.0'];
         else
            word = num2str(R);
         end
         
         if(mod(R,.1)==0) 
         	set(text(Gr(MIN),0,word),'Rotation',90,'HorizontalAlignment','left','VerticalAlignment','bottom');
      	end
   	elseif(R<=2)
         if(mod(R,.2)==0)
            if(mod(R,1)==0)
               word = [num2str(R) '.0'];
            else
               word = num2str(R);
            end
         	set(text(Gr(MIN),0,word),'Rotation',90,'HorizontalAlignment','left','VerticalAlignment','bottom');
      	end
   	elseif(R<=5)
      	if(mod(R,1)==0)
         	set(text(Gr(MIN),0,[num2str(R) '.0']),'Rotation',90,'HorizontalAlignment','left','VerticalAlignment','bottom');
      	end
   	else
      	if(mod(R,10)==0)
         	set(text(Gr(MIN),0,num2str(R)),'Rotation',90,'HorizontalAlignment','left','VerticalAlignment','bottom');
      	end
      end   
      elseif(R==.25 | R==.35 | R==.45)
      plot(Gr(MIN2:index),r_L_a3,'b');
      plot(Gr(MIN2:index),r_L_b3,'b');		
   end
   
   if(R==.05 | (R>.149 & R<.151))
      plot(Gr(MIN2:MAX2),r_L_a2(length(Gr(MIN2:MAX2))-length(r_L_a2)+1:length(r_L_a2)),'b');
      plot(Gr(MIN2:MAX2),r_L_b2(length(Gr(MIN2:MAX2))-length(r_L_b2)+1:length(r_L_b2)),'b');		
   end
   
end

%equations were derived using the symbolic toolbox as follows
%solve('2*Gi/((1-Gr)^2+Gi^2)=x','Gi')
%bound was derived as follows
%solve('1-X^2+2*X^2*Gr-X^2*Gr^2=0','Gr')
%solve('1/2/X*(2+2*(1-X^2+2*X^2*Gr-X^2*Gr^2)^(1/2))=(1-Gr^2)^(1/2)','Gr')
for X = interval,
   inter_bound = (-1+X^2)/(X^2+1); %intersection with unit circle: all values must be less than this\
   imag_bound =  (-1+X)/X; %boundary of imagination: all values must be greater than this
   angle_point = 0;
   if(inter_bound ~= 0)
	   angle_point = sqrt(1-inter_bound^2)/inter_bound;
   end
   
   imag_bound_y =  1/2/X*(-2+2*(1-X^2+2*X^2.*inter_bound-X^2.*inter_bound.^2).^(1/2));
   
   imag_rad = (imag_bound^2 + imag_bound_y^2)^(1/2);
   condition = imag_rad < 1;
   if(inter_bound > 1)
      inter_bound = 1;
   elseif(inter_bound < -1)
      imag_bound=-1;
   end
   
   if(imag_bound > 1)
      imag_bound = 1;
   elseif(imag_bound < -1)
      imag_bound=-1;
   end
   
   %used solve function to find intersection of appropriate circle with corresponding hyperbolics
   %solve('-1/(R+1)*(-(R+1)*(R-2*R*Gr+R*Gr^2-1+Gr^2))^(1/2)=1/2/X*(-2+2*(1-X^2+2*X^2*Gr-X^2*Gr^2)^(1/2))','Gr')