function [magout,phase,w] = nichols(a,b,c,d,iu,w)
%NICHOLS Nichols frequency response for continuous-time linear systems.
%	NICHOLS(A,B,C,D,IU) produces a Nichols plot from the single
%	input IU to all the outputs of the continuous state-space system 
%	(A,B,C,D).  IU is an index into the inputs of the system and 
%	specifies which input to use for the Nichols response.  The 
%	frequency range and number of points are chosen automatically.
%
%	NICHOLS(NUM,DEN) produces the Nichols plot for the polynomial 
%	transfer function G(s) = NUM(s)/DEN(s) where NUM and DEN contain
%	the polynomial coefficients in descending powers of s. 
%
%	NICHOLS(A,B,C,D,IU,W) or NICHOLS(NUM,DEN,W) uses the user-supplied
%	freq. vector W which must contain the frequencies, in radians/sec,
%	at which the Nichols response is to be evaluated.  When invoked 
%	with left hand arguments,
%		[MAG,PHASE,W] = NICHOLS(A,B,C,D,...)
%		[MAG,PHASE,W] = NICHOLS(NUM,DEN,...) 
%	returns the frequency vector W and matrices MAG and PHASE (in 
%	degrees) with as many columns as outputs and length(W) rows.  No 
%	plot is drawn on the screen.  Draw the Nichols grid with NGRID.
%
%	See also: LOGSPACE,NGRID,MARGIN,BODE, and NYQUIST.

%	Clay M. Thompson  7-6-90
%	Revised ACWG 6-21-92
%	Copyright (c) 1986-93 by the MathWorks, Inc.

if nargin==0, eval('exresp(''nichols'')'), return, end

error(nargchk(2,6,nargin));

% --- Determine which syntax is being used ---
if (nargin==1),
	error('Wrong number of input arguments.');

elseif (nargin==2),	% Transfer function form without frequency vector
	num = a; den = b;
	w = freqint2(num,den,30);
	[ny,nn] = size(num); nu=1;

elseif (nargin==3),	% Transfer function form with frequency vector
	num = a; den = b;
	w = c;
	[ny,nn] = size(num); nu=1;

elseif (nargin==4),	% State space system, w/o iu or frequency vector
	error(abcdchk(a,b,c,d));
	w = freqint2(a,b,c,d,30);
	[iu,nargin,mag,phase]=mulresp('nichols',a,b,c,d,w,nargout,1);
	if ~iu, if nargout, magout = mag; end, return, end
	[ny,nu] = size(d);

elseif (nargin==5),	% State space system, with iu but w/o freq. vector
	error(abcdchk(a,b,c,d));
	w = freqint2(a,b,c,d,30);
	[ny,nu] = size(d);

else
	error(abcdchk(a,b,c,d));
	[ny,nu] = size(d);

end

if nu*ny==0, phase=[]; w=[]; if nargout~=0, magout=[]; end, return, end

% --- Evaluate the frequency response ---
if (nargin==2)|(nargin==3),
	g = freqresp(num,den,sqrt(-1)*w);
else
	g = freqresp(a,b,c,d,iu,sqrt(-1)*w);
end

mag = abs(g);
phase = -180+180./pi*atan2(-imag(g),-real(g));

% If no left hand arguments then plot graph.
if nargout==0
	status = ishold;
	magdb = 20*log10(mag);

	% Find the points where the phase wraps.  Plot each branch separately.
	% The following code is based on the algorithm in the unwrap function.

	cutoff = 200;
	[m, n] = size(phase);				% Assume column orientation.

	pmin = min(phase); pmin = pmin(ones(m, 1), :);% To force REM to behave.
	p = rem(phase - pmin, 360) + pmin;            % Phases modulo 360.

	dp = [p(1,:);diff(p)];			% Differentiate phases.
	jumps = (dp > cutoff) + (dp < -cutoff);	% Array of jumps locations
	if n==1, jvec = jumps; else jvec = (sum(jumps')~=0)'; end
	seg = cumsum(jvec);				% Used to identify segments

	jloc = find(jvec~=0);

	% Strip off each continuous piece and plot separately
	pj=[]; mj=[]; jj=[];
	for k=0:seg(length(seg)),
		mseg = magdb(seg==k,:); pseg = p(seg==k,:);

		% Add last point of previous segments so plots are continuous
		if ~isempty(pj), pj(jj)=pseg(1,jj); mj(jj)=mseg(1,jj); end
		pseg=[pj;pseg]; mseg=[mj;mseg];

		% Remember last point
		[np,mp] = size(pseg);
		if np~=0,
		    pj = pseg(np,:); mj = mseg(np,:); 
		    if k<length(jloc), jj=(jumps(jloc(k+1),:)~=0); end
		end

		if np>1,
		    plot(pseg,mseg), 		% Plot segment
		    hold on
		end
	end

	axstate = axis('state');
	if axstate(1)=='a'  % axis is set to auto range
		set(gca,'xlim',[-360 0]);
	end
	set(gca,'xtick',[-360, -270, -180, -90 0]);
	xlabel('Open-Loop Phase (deg)')
	ylabel('Open-Loop Gain (db)')

	if ~status, hold off, end	% Return hold to previous status
	return % Suppress output
end

% Uncomment the following line for decibels, but be warned that the 
% MARGIN function will not work with decibel data.

% mag = 20*log10(mag);
magout = mag;