mcirc.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function c=mcirc(m)
%MCIRC  Calculates M-circles
%       MCIRC(M) returns a complex matrix with
%       one circle per column.
%       M is the vector of M-circle values in dB.
%       See also NCIRC,CIRC,ARC

%       Values of M within 2*eps of 0dB are treated as 0dB

%       Dr M.P. Ford 4th August 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
% MRN026

if min(size(m))~=1
   error('M must be a vector')
end
m=10 .^(m/20);
angone=[(1:1:16)./8 2.5 (7:5:22)./2 20 30 60 100];
ly=length(angone);
ir=ly:-1:1;
minusang=-(angone(ir));
%
% Note all these three angle vectors must have the same number of elements
% 0dB and +,- dB's treated separately to ensure correct plotting on
% Nichols chart.
%
angone=[minusang,0,angone].';

angplus=[1599 (+1500:-100:+500) (+400:-40:+80) (+40:-10:-30)...
	  (-40:-40:-360) (-400:-100:-1500) -1599].';
angplus=angplus.*(pi/1600);
angminus=angplus(51:-1:1);

lm=length(m);
lang=length(angone);
c=zeros(lang,lm);
for i=1:lm
 if (m(i)<1+2*eps)&(m(i)>1-2*eps)
    %   calculate M=1 lines
    c(:,i)=-0.5+sqrt(-1)*angone;
 elseif m(i)>1
   c(:,i)=(-m(i)^2/(m(i)^2-1))+(m(i)/(m(i)^2-1)).*exp(sqrt(-1).*angminus);
 else  % m(i)<1
   c(:,i)=(-m(i)^2/(m(i)^2-1))+(m(i)/(m(i)^2-1)).*exp(sqrt(-1).*angplus);
 end;
end;