function [magout,phase,w] = dnichols(a,b,c,d,Ts,iu,w)
%DNICHOLS Nichols frequency response for discrete-time linear systems.
%	DNICHOLS(A,B,C,D,Ts,IU) produces a Nichols plot from the single 
%	input IU to all the outputs of the discrete 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.  Ts is the
%	sample period.  The frequency range and number of points are 
%	chosen automatically.
%
%	DNICHOLS(NUM,DEN,Ts) produces the Nichols plot for the polynomial 
%	transfer function G(z) = NUM(z)/DEN(z) where NUM and DEN contain
%	the polynomial coefficients in descending powers of z. 
%
%	DNICHOLS(A,B,C,D,Ts,IU,W) or DNICHOLS(NUM,DEN,Ts,W) uses the user-
%	supplied freq. vector W which must contain the frequencies, in 
%	rad/sec, at which the Nichols response is to be evaluated.  
%	Aliasing will occur at frequencies greater than the Nyquist 
%	frequency (pi/Ts). With left hand arguments,
%		[MAG,PHASE,W] = DNICHOLS(A,B,C,D,Ts,...)
%		[MAG,PHASE,W] = DNICHOLS(NUM,DEN,Ts,...) 
%	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,SEMILOGX,MARGIN,DBODE, and DNYQUIST.

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

if nargin==0, eval('dexresp(''dnichols'')'), return, end

error(nargchk(3,7,nargin));

% --- Determine which syntax is being used ---
if (nargin==3),		% Transfer function without frequency vector
	num = a; den = b;
	Ts=c;
	w=dfrqint2(num,den,Ts,30);
	[ny,nn] = size(num); nu = 1;

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

elseif (nargin==5),	% State space system without iu or freq. vector
	error(abcdchk(a,b,c,d));
	w=dfrqint2(a,b,c,d,Ts,30);
	[iu,nargin,mag,phase]=dmulresp('dnichols',a,b,c,d,Ts,w,nargout,1);
	if ~iu, if nargout, magout = mag; end, return, end
	[ny,nu] = size(d);

elseif (nargin==6),	% State space system, with iu but w/o freq. vector
	error(abcdchk(a,b,c,d));
	w=dfrqint2(a,b,c,d,Ts,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

% Compute frequency response
if (nargin==3)|(nargin==4)
	g=freqresp(num,den,exp(sqrt(-1)*w*Ts));
else
	g=freqresp(a,b,c,d,iu,exp(sqrt(-1)*w*Ts));
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('Phase (deg)')
	ylabel('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;