plotsmithchart.m

使用matlab画的smith圆图

   %The following conditional tree creates the internal bounding between the two types of curves for variable resolution
   if(X<.2)
      if(mod(X,.1)==0) 
         max_bound = (-1+X^2+2^2)/(X^2+2^2+2*2+1);
      elseif(mod(X,.02)==0)
         max_bound = (-1+X^2+.5^2)/(X^2+.5^2+2*.5+1);
      else
         max_bound = (-1+X^2+.2^2)/(X^2+.2^2+2*.2+1);
      end
   elseif(X<1)
      if(mod(X,.2)==0) 
         max_bound = (-1+X^2+5^2)/(X^2+5^2+2*5+1);
      elseif(mod(X,.1)==0) 
         max_bound = (-1+X^2+2^2)/(X^2+2^2+2*2+1); 
      elseif(X<.5) 
         max_bound = (-1+X^2+.5^2)/(X^2+.5^2+2*.5+1); 
      else 
         max_bound = (-1+X^2+1^2)/(X^2+1^2+2*1+1); 
      end      
   elseif(X<5)
      if(mod(X,2)==0)
         max_bound = (-1+X^2+20^2)/(X^2+20^2+2*20+1); 
      elseif(mod(X,1)==0)
         max_bound = (-1+X^2+10^2)/(X^2+10^2+2*10+1); 
      elseif(X>2)
         max_bound = (-1+X^2+5^2)/(X^2+5^2+2*5+1);
      else
         if(mod(X,.2)==0)
   	   	max_bound = (-1+X^2+5^2)/(X^2+5^2+2*5+1);
         else
            max_bound = (-1+X^2+2^2)/(X^2+2^2+2*2+1);
         end
      end
   elseif(X<10)
      if(mod(X,2)==0)
         max_bound = (-1+X^2+20^2)/(X^2+20^2+2*20+1); 
      else
         max_bound = (-1+X^2+10^2)/(X^2+10^2+2*10+1); 
      end
   else
      if(X==10|X==20)
         max_bound = (-1+X^2+50^2)/(X^2+50^2+2*50+1); 
      elseif(X==50)
         max_bound = 1;
      elseif(X<20)
         max_bound = (-1+X^2+20^2)/(X^2+20^2+2*20+1); 
      else
         max_bound = (-1+X^2+50^2)/(X^2+50^2+2*50+1); 
      end
   end 
  
   inter_index = ceil((inter_bound+1)*(MAX-1)/2+1);
   imag_index = ceil((imag_bound+1)*(MAX-1)/2+1);   
   index4 = ceil((max_bound+1)*(MAX-1)/2+1);
   
   index1 = max(inter_index,imag_index); %maximum index for c,d
   index2 = min(imag_index,inter_index); %minimum index for c,d
   
   if(condition)
      index3=imag_index;
   else
      index3=inter_index;
   end
  
   actual_value1 = Gr(index1);
   actual_value2 = Gr(index2);
   actual_value3 = Gr(index3);
   actual_value4 = Gr(index4);
   
   if((actual_value1 > inter_bound & index1 == inter_index)|(actual_value1 > imag_bound & index1 == imag_index))
      index1 = index1 - 1;
   end
   if((actual_value2 < inter_bound & index2 == inter_index)|(actual_value2 < imag_bound & index2 == imag_index))
      index2 = index2 + 1;
   end   
   if((actual_value3 < inter_bound & index3 == inter_index)|(actual_value3 < imag_bound & index3 == imag_index))
   	index3 = index3 + 1;
   end   
   if(actual_value4 > max_bound)
      index4 = index4 - 1;
   end
  
   MIN=index2;
   MAX2=index1;
   MAX3=index4;
   MIN2 = index3;
   
 %  actual_value1 = Gr(MIN);
 %  actual_value2 = Gr(MAX2);
 %  MIN=1;
 %  MAX2=MAX;
 %  MIN2=1;
      
   x_L_a = real(1/2/X*(-2+2*(1-X^2+2*X^2.*Gr(MIN2:MAX3)-X^2.*Gr(MIN2:MAX3).^2).^(1/2)));   
   x_L_b = real(1/2/X*(2-2*(1-X^2+2*X^2.*Gr(MIN2:MAX3)-X^2.*Gr(MIN2:MAX3).^2).^(1/2)));
   x_L_c= real(1/2/X*(2+2*(1-X^2+2*X^2.*Gr(MIN:MAX2)-X^2.*Gr(MIN:MAX2).^2).^(1/2)));
   x_L_d= real(1/2/X*(-2-2*(1-X^2+2*X^2.*Gr(MIN:MAX2)-X^2.*Gr(MIN:MAX2).^2).^(1/2)));   
   
   if(MIN2<MAX3)
   	x_L_c(1)=x_L_b(1);
   	x_L_d(1)=x_L_a(1);
   end
   
   check1 = abs(round(10000*1/2/X*(-2-2*(1-X^2+2*X^2*inter_bound-X^2*inter_bound^2)^(1/2))));
   check2 = abs(round(10000*(1-inter_bound^2)^(1/2)));

   if(imag_bound > -1 & check1 == check2)
      plot(Gr(MIN:MAX2),x_L_c,'g')
      plot(Gr(MIN:MAX2),x_L_d,'g') 
   end
      
   plot(Gr(MIN2:MAX3),x_L_a,'g')
   plot(Gr(MIN2:MAX3),x_L_b,'g')
   
   condition = Gr(MIN2)^2 + x_L_d(1)^2 > .985;   
   if(X<=1)      
      if(mod(X,.1)==0) 
         if(mod(X,1)==0)
            word = [num2str(X) '.0'];
         else
            word = num2str(X);
         end
         if(X==1)
            angle = 90;
         else
            angle = -atan(angle_point)*180/pi;
         end
         set(text(Gr(MIN2),x_L_d(1),word),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
         set(text(Gr(MIN2),-x_L_d(1),word),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');
       	if(mod(X,.2)==0)
            xval=X^2/(X^2+4);
            yval = 1/2/X*(-2+2*(1-X^2+2*X^2*xval-X^2*xval^2)^(1/2));
            angle = -atan(yval/(.5-xval))*180/pi;
				set(text(xval,yval,word),'Rotation',angle,'HorizontalAlignment','left','VerticalAlignment','bottom');
            set(text(xval,-yval,word),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom')
         end   
      end
   elseif(X<=2)
      if(mod(X,.2)==0)
         if(mod(X,1)==0)
            word = [num2str(X) '.0'];
         else
            word = num2str(X);
         end
         if(condition)
            angle = -atan(angle_point)*180/pi+180;
            set(text(Gr(MIN2),x_L_a(1),word),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
            set(text(Gr(MIN2),-x_L_a(1),word),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');            
         else
            angle = -atan(angle_point)*180/pi+180;            
      		set(text(Gr(MAX2),x_L_d(length(x_L_d)),word),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
      		set(text(Gr(MAX2),-x_L_d(length(x_L_d)),word),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');            
         end
      end
   elseif(X<=5)
      if(mod(X,1)==0)
         if(condition)
            angle = -atan(angle_point)*180/pi+180;
      		set(text(Gr(MIN2),x_L_a(1),[num2str(X) '.0']),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
      		set(text(Gr(MIN2),-x_L_a(1),[num2str(X) '.0']),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');            
         else
            angle = -atan(angle_point)*180/pi+180;            
      		set(text(Gr(MAX2),x_L_d(length(x_L_d)),[num2str(X) '.0']),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
      		set(text(Gr(MAX2),-x_L_d(length(x_L_d)),[num2str(X) '.0']),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');            
         end
      end
   else
      if(mod(X,10)==0)
         if(condition)
            angle = -atan(angle_point)*180/pi+180;            
      		set(text(Gr(MIN2),x_L_a(1),num2str(X)),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
      		set(text(Gr(MIN2),-x_L_a(1),num2str(X)),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');            
         else
            angle = -atan(angle_point)*180/pi+180;            
            set(text(Gr(MAX2),x_L_d(length(x_L_d)),num2str(X)),'Rotation',angle,'VerticalAlignment','bottom','HorizontalAlignment','left');
        		set(text(Gr(MAX2),-x_L_d(length(x_L_d)),num2str(X)),'Rotation',-angle+180,'HorizontalAlignment','right','VerticalAlignment','bottom');
   		end
      end
   end
end
%plot imaginary axis
plot(zeros(1,length(Gr)),Gr,'r');


wavelengths = 0:.01:.5;
angle = linspace(pi,-pi,length(wavelengths));
wave_circle = 1.05*exp(j*phaseAngle);
plot(real(wave_circle),imag(wave_circle),'r');

for i=1:length(wavelengths)-1,
   x=real(1.025*exp(j*angle(i)));
   y=imag(1.025*exp(j*angle(i)));
   if(x>0)
      rot_angle=atan(y/x)*180/pi-90;
   else
      rot_angle=atan(y/x)*180/pi+90;
   end
   if(wavelengths(i)==0)
      word = '0.00';
   elseif(mod(wavelengths(i),.1)==0)
      word = [num2str(wavelengths(i)) '0'];
   else
      word = num2str(wavelengths(i));
   end
   set(text(x,y,word),'Rotation',rot_angle,'VerticalAlignment','middle','HorizontalAlignment','center');
end

%plot reflection coefficient and line of intersection only if arguments are present

if(nargin == 2)
	radius = abs((Zl-Zo)/(Zl+Zo));
	Load_circle=radius*exp(j*phaseAngle);
	plot(real(Load_circle),imag(Load_circle),'r');

	slope = (-(1-radius^2)^(1/2)*radius)/(radius^2);

	value=1/(1+slope^2)^(1/2);
	MAX2 = ceil((value+1)*(MAX-1)/2+1);
	actual_value = Gr(MAX2);
	if(actual_value>value)
  		MAX2 = MAX2 - 1;
	end

	%plot line of intersection
	line = slope*Gr(fix(MAX/2):MAX2);
	plot(Gr(fix(MAX/2):MAX2),line,'r');
end