% function vplot([axis],vmat1,vmat2,vmat3, ...)
% function vplot([axis],vmat1,'linetype1',vmat2,'linetype2',...)
%
%   Plot one or more varying matrices.  The syntax is
%   the same as the MATLAB plot command except that all data
%   is contained in VMATi, and the axes are specified by AXIS.
%
%   The (optional) axis argument must be one of:
%
%      'iv,d'       matin .vs. independent variable (default option)
%      'iv,m'       magnitude .vs. independent variable
%      'iv,lm'      log(magnitude) .vs. independent variable
%      'iv,p'       phase .vs. independent variable
%      'liv,d'      matin .vs. log(independent variable)
%      'liv,m'      magnitude .vs. log(independent variable)
%      'liv,lm'     log(magnitude) .vs. log(independent variable)
%      'liv,p'      phase .vs. log(independent variable)
%      'ri'         real .vs. imaginary  (parametrized by indep variable)
%      'nyq'        real .vs. imaginary  (parametrized by indep variable)
%      'nic'        Nichols chart
%      'bode'       Bode magnitude and phase plots
%
%
%    See also: PLOT

%====================================================================
%
%	MAGNITUDE/PHASE/FREQUENCY UNIT DEFINITION
%
%	Choose only one of these commands to set up the
%	phase units.  The others must be commented out.

function vplot(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17, ...
	a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33)

phasefunc = '180/pi*angle(';	    % degrees (-180 to 180)
%phasefunc = 'angle(';		    % radians (-pi to pi)
%phasefunc = '180/pi*negangle(';    % degrees (-360 to 0)
%phasefunc = 'negangle(';	    % radians (-2pi to 0)

%	Choose the appropriate set of commands to setup the
%	magnitude scaling for log magnitude commands. 

%magfunc = '20*log10(';		% Scale in dB
%mplt1 = 'semilogx(';		% You must select all 3
%mplt2 = 'plot(';		% commands here.

magfunc = 'abs(';		% Scale in log magnitude
mplt1 = 'loglog(';		% Again all 3 commands are
mplt2 = 'semilogy(';		% required.

%	Note that 20*log10(x) gives a complex result when x is
%	complex.  The real part is equal to 20*log10(abs(x)) and
%	the MATLAB plot function plots only the real part.
%
%	The following section customizes the 'liv,..' plots.
%	In almost all cases a logarithmic independent variable
%	refers to frequency -- this section chooses between displaying
%	in radians/second and Hertz.  If you intend to use log independent
%	variable plots for something other than frequency then you should
%	select the default scaling of 1.  Note that the output
%	from frsp is in radians/second.  Choosing this throughout is
%	very strongly encouraged.

ivscale = 1;			% frequency in radians/sec
%ivscale = 1/(2*pi);		% frequency in Hertz
%ivscale = 1/pi;		% discrete normalized frequency

%	The following section customizes the Bode plot labels.
%	Make sure that they agree with your choices for the
%	log magnitude, phase, and frequency units.

%toplabel = 'Magnitude (dB)';
toplabel = 'Log Magnitude';

bottomlabel = 'Phase (degrees)';
%bottomlabel = 'Phase (radians)';

ivlabel = 'Frequency (radians/sec)';
%ivlabel = 'Frequency (Hertz)';
%ivlabel = 'Normalized Frequency';

%===================================================================

n = nargin;
if n == 0,
	disp('usage: vplot([axis],vmat1,vmat2,...)')
	disp('       vplot([axis],vmat1,linetype1,...)')
	return
	end

if n > 33,
	disp('    Sorry, vplot is limited to 33 arguments')
	return
	end

%	Attempt to see if the first argument is the axis string.
%	If it isn't use the default axis setting.  argstart denotes
%	the start of the data arguments.

if isstr(a1),
	argstart = 2;
	axis = a1;
else
	axis = 'iv,d';
	argstart = 1;
	end

%	Determine the required plot type and what functions have
%	to be applied to the data.  Three types of plot are
%	distinguished by the variable plttype:
%
%		1:	function(data) vs. independent variable
%		2:	function(data) vs. function(data)
%		3:	customized double plot (each of type 1)
%

if strcmp(axis,'iv,d'),
	func = '(';
	ivscale = 1;
	pltfunc = 'plot(';
	plttype = 1;
elseif strcmp(axis,'iv,lm'),
	func = magfunc;
	ivscale = 1;
	pltfunc = mplt2;
	plttype = 1;
elseif strcmp(axis,'iv,m'),
	func = 'abs(';
	ivscale = 1;
	pltfunc = 'plot(';
	plttype = 1;
elseif strcmp(axis,'iv,p'),
	func = phasefunc;
	ivscale = 1;
	pltfunc = 'plot(';
	plttype = 1;
elseif strcmp(axis,'liv,d'),
	func = '(';
	pltfunc = 'semilogx(';
	plttype = 1;
elseif strcmp(axis,'liv,lm'),
	func = magfunc;
	pltfunc = mplt1;
	plttype = 1;
elseif strcmp(axis,'liv,m'),
	func = 'abs(';
	pltfunc = 'semilogx(';
	plttype = 1;
elseif strcmp(axis,'liv,p'),
	func = phasefunc;
	pltfunc = 'semilogx(';
	plttype = 1;
elseif strcmp(axis,'ri');
	func1 = 'imag(';
	func2 = 'real(';
	pltfunc = 'plot(';
	plttype = 2;
elseif strcmp(axis,'nyq');
	func1 = 'imag(';
	func2 = 'real(';
	pltfunc = 'plot(';
	plttype = 2;
elseif strcmp(axis,'nic');
	func1 = '20*log10(';
	func2 = '360/(2*pi)*negangle(';
	pltfunc = 'plot(';
	plttype = 2;
elseif strcmp(axis,'bode'),
	tfunc = magfunc;
	bfunc = phasefunc;
	tpltfunc = mplt1;
	bpltfunc = 'semilogx(';
	plttype = 3;
else
	errstr = ['    unrecognized axis specification: ',axis];
	error(errstr)
	end

%	For the particular plttype parse the argument list and
%	create the string to pass as an argument list to a standard
%	MATLAB plot routine.  The functions determined above (eg,
%	abs, angle, etc) are applied to the appropriate data.
%	If a string is found it is assumed to specify a line type
%	and is included directly in the argument list.  Constant
%	matrices are assumed to be constant over independent variables.
%	The minimum and maximum independent variables are tracked so
%	that constants can be plotted out over the full range of the
%	graph.

if plttype == 1,
    ivc = [];
    for i = argstart:n,
	eval(['strtest=isstr(a',num2str(i),');'])
	if strtest,
	    str = [str,',a',num2str(i)];
	else
	    eval(['a = a',num2str(i),';']);
	    [nt,nr,nc,npts] = minfo(a);
	    if nt == 'vary',
		[dat,ptr,iv] = vunpck(a);
		dat = dat.';
		dat = reshape(dat,nr*nc,npts).';
		if npts == 1,
			dat = [dat;dat];
			iv = [iv;iv];
			end
		eval(['dat',num2str(i),' = ',func,'dat);']);
		eval(['iv',num2str(i),' = ivscale*iv;']);
		if isempty(ivc),
		    ivc = [max(ivscale*iv); min(ivscale*iv)];
		else
		    ivc = [max(ivc(1),max(ivscale*iv));...
				min(ivc(2),min(ivscale*iv))];
		    end
		str = [str,',iv',num2str(i),',dat',num2str(i)];
	    elseif nt == 'syst',
		error('cannot plot system matrix')
		break
	    else 
		dat = a.';
		dat = reshape(dat,nr*nc,1).';
		dat = [dat; dat];
		eval(['dat',num2str(i),' = ',func,'dat);']);
		str = [str,',ivc,dat',num2str(i)];
		end
	    end
	end

    str = str(2:length(str));
    str = [pltfunc, str, ')'];
%    disp(str)
    eval(str)
elseif plttype == 2,
    for i = argstart:n,
	eval(['strtest=isstr(a',num2str(i),');'])
	if strtest,
	    str = [str,',a',num2str(i)];
	else
	    eval(['a = a',num2str(i),';']);
	    [nt,nr,nc,npts] = minfo(a);
	    if nt == 'vary' | nt == 'cons',
		[dat,ptr,iv] = vunpck(a);
		dat = dat.';
		if npts == 0,
		    npts = 1;
		    end
		dat = reshape(dat,nr*nc,npts).';
		if npts == 1,
			dat = [dat;dat];
			iv = [iv;iv];
			end
		eval(['dat',num2str(i),' = ',func1,'dat);']);
		eval(['iv',num2str(i),' = ',func2,'dat);']);
		str = [str,',iv',num2str(i),',dat',num2str(i)];
	    else nt == 'syst',
		error('cannot plot system matrix')
		break
		end
	    end
	end

    str = str(2:length(str));
    str = [pltfunc, str, ')'];
%    disp(str)
    eval(str)
else,
    ivc = [];
    for i = argstart:n,
	eval(['strtest=isstr(a',num2str(i),');'])
	if strtest,
	    str1 = [str1,',a',num2str(i)];
	    str2 = [str2,',a',num2str(i)];
	else
	    eval(['a = a',num2str(i),';']);
	    [nt,nr,nc,npts] = minfo(a);
	    if nt == 'vary',
		[dat,ptr,iv] = vunpck(a);
		dat = dat.';
		dat = reshape(dat,nr*nc,npts).';
		if npts == 1,
			dat = [dat;dat];
			iv = [iv;iv];
			end
		eval(['tdat',num2str(i),' = ',tfunc,'dat);']);
		eval(['bdat',num2str(i),' = ',bfunc,'dat);']);
		eval(['iv',num2str(i),' = ivscale*iv;']); 
		if isempty(ivc),
		    ivc = [max(ivscale*iv); min(ivscale*iv)];
		else
		    ivc = [max(ivc(1),max(ivscale*iv));...
				min(ivc(2),min(ivscale*iv))];
		    end
		str1 = [str1,',iv',num2str(i),',tdat',num2str(i)];
		str2 = [str2,',iv',num2str(i),',bdat',num2str(i)];
	    elseif nt == 'syst',
		error('cannot plot system matrix')
		break
	    else 
		dat = a.';
		[nr,nc] = size(a);
		dat = reshape(dat,nr*nc,1).';
		dat = [dat; dat];
		eval(['tdat',num2str(i),' = ',tfunc,'dat);']);
		eval(['bdat',num2str(i),' = ',bfunc,'dat);']);
		str1 = [str1,',ivc,tdat',num2str(i)];
		str2 = [str2,',ivc,bdat',num2str(i)];
		end
	    end
	end

    str1 = str1(2:length(str1));
    str1 = [tpltfunc, str1, ')'];
    str2 = str2(2:length(str2));
    str2 = [bpltfunc, str2, ')'];
%    disp(str1)
%    disp(str2)
    subplot(211)
    eval(str1)
    ylabel(toplabel)
    xlabel(ivlabel)
    subplot(212)
    eval(str2)
    ylabel(bottomlabel)
    xlabel(ivlabel)
    subplot(111)
    end

%------------------------------------------------------------------------

%
% Copyright MUSYN INC 1991,  All Rights Reserved
