function out = MultiSys(arg1,arg2,flag,parms,funflag)
% MULTISYS  illustrates the response of a compound dynamic system
% consisting of first- and/or second-order dynamic systems coupled
% together so that the output of one is the input of the next.
% The startup response to the default input is shown, and also the
% steady response to a harmonic input.  All the parameters of the
% dynamic systems may be edited.
% MULTISYS(t,F) uses columns of time t and forcing function F(t),
% piecewise linear data to replace the default input.
% MULTISYS(tspan,'F') uses the user-defined function 'F' for the
% forcing function.  tspan should be either a two-element array of
% the form [t0 tfinal], or else an array of all the times at which
% display points are desired.  This is consistent with ODE45.

if nargin==5	% ODEFile functionality
	t=arg1;		r=arg2;
	rp=zeros(size(r));
	first=1;
	for k=1:parms.Nstage
		% interstage linkage system
		if k==1		% input to first system
			if funflag		% function input
				Fk=feval(parms.F,t);
			else			% piecewise linear input
				Fk=interp1(parms.t,parms.F,t);
			end
		else		% previous output coupled to
			Fk=r(first-parms.order(k-1));	% current stage
		end	
		% dynamic equations
		if parms.order(k)==1		% First-order system
			rp(first)=(parms.K(k)*Fk - r(first))/parms.tau(k);
		else			% Second-order system
			fp1=first+1;
			rp(first)=r(fp1);
			rp(fp1)=[parms.A(k),parms.B(k),parms.C(k)]* ...
				[r(fp1);r(first);Fk];
		end
		first=first+parms.order(k);
	end
	out=rp;

elseif (nargin==1 & isstr(arg1))	% GUI callback functionality
	action=arg1;
	q=get(gcbf,'UserData');
	switch (action)
	case 'popup'
		switch q.mode
		case 1					% Exit General mode
			legend off
			cla
			set(q.genclump,'Visible','off')
		case 2					% Exit Step mode
			legend off
			cla
			set(q.stepclump,'Visible','off')
		case 3					% Exit Harmonic mode
			legend off
			cla
			set(q.harmclump,'Visible','off')
		case 4					% Exit Design mode
			axes(q.magax),	legend off,		cla
			axes(q.phaseax),	legend off,	cla
			axes(q.poleax),		cla
			set(q.desclump,'Visible','off')
		end
		q.mode=get(gcbo,'Value');
		switch q.mode
		case 1					% Enter General mode
			set(q.genclump,'Visible','on')
			q.lhs=gendisp(q);
		case 2					% Enter Step mode
			set(q.stepclump,'Visible','on')
			stepdisp(q);
		case 3					% Enter Harmonic mode
			set(q.harmclump,'Visible','on')
			harmdisp(q);
		case 4					% Enter Design mode
			set(q.desclump,'Visible','on')
			desdisp(q);
		end
	case 'return'	% build output format query window
		choice=questdlg('Select Output Time Basis','Numeric Output', ...
			'Resample','Display','Display');
		switch (choice)
		case 'Display'	% get output already calculated for line using
			t=get(q.lhs(1),'XData')';	% ode45 internal time basis
			yc=get(q.lhs,'YData');
			y=zeros(length(yc{1}),length(yc));
			for k=1:length(yc)
				y(:,k)=yc{k}';		% extract y from cell array
			end
		case 'Resample'	% recalculate output using new sample rate
			T=q.t(end)-q.t(1);
			fs=(length(q.t)-1)/T;
			fstr=inputdlg('Enter Sampling Frequency','fs',1,{num2str(fs)});
			fs=str2num(fstr{1});
			trange=q.t(1):1/fs:q.t(end);
			[t,y]=solvesys(q,trange);
		end
		assignin('base','tout',t)
		assignin('base','yout',y)
		disp('Returned tout and yout to workspace')
	case {'f','A'}
	   eval(['q.',action,'=str2num(get(gcbo,''String''));'])
	   harmdisp(q)
	case {'Fi','Ff'}
	   eval(['q.',action,'=str2num(get(gcbo,''String''));'])
	   stepdisp(q)
	case 'popstage'
		q.stage=get(gcbo,'Value');
		destxt(q)
	case 'new'
		str=inputdlg('Create a New System','Number of Stages',1,{'3'});
		if size(str)
			q.Nstage=str2num(str{1});	N=q.Nstage;
			coefs=buttcoefs(2*N-1,1);		% see function below
			q.stage=1;	q.order=coefs.order;
			q.tau=coefs.tau;	q.K=coefs.K;	q.fn=coefs.fn;
			q.zeta=coefs.zeta;	q.y0=zeros(1,N);	q.yp0=q.y0;
			str='1';
			for k=2:q.Nstage	str=[str,'|',int2str(k)];	end
			set(q.popstage,'String',str,'Value',1)
		end
		desdisp(q)
	case 'edit'
		defs={'First','Second'};
		s=q.stage;
		titlestr=['Edit Stage ',int2str(s)];
		choice=questdlg('Select Order',titlestr, ...
			'First','Second',defs{q.order(s)});
		switch (choice)
		case 'First'
			prompt={'tau';'K';'y0'};
			defs={num2str(q.tau(s));num2str(q.K(s));num2str(q.y0(s))};
			str=inputdlg(prompt,titlestr,1,defs);
			if (~size(str)) break, end
			q.order(s)=1;
			q.tau(s)=eval(str{1});
			q.K(s)=eval(str{2});
			q.y0(s)=eval(str{3});
		case 'Second'
			prompt={'fn';'zeta';'K';'y0';'y''0'};
			defs={num2str(q.fn(s));num2str(q.zeta(s));num2str(q.K(s));
				num2str(q.y0(s));num2str(q.yp0(s))};
			str=inputdlg(prompt,titlestr,1,defs);
			if (~size(str)) break, end
			q.order(s)=2;
			q.fn(s)=eval(str{1});
			q.zeta(s)=eval(str{2});
			q.K(s)=eval(str{3});
			q.y0(s)=eval(str{4});
			q.yp0(s)=eval(str{5});
		end
		desdisp(q)
	case 'q'
		out=q;
	case 'help'
    	helpwin('aboutmultisys')
	otherwise
		disp(['callback ',action,' unknown'])
	end
	set(gcbf,'UserData',q)

%	Start-ups
elseif (nargin==0 | (nargin==2 & isnumeric(arg1) & isnumeric(arg2)) | ...
		(nargin==2 & isstr(arg2)))
	funflag=0;
	if nargin==0					% easy startup
		t=[0 0.499 0.5 0.999 1 1.499 1.5 1.999 2 2.499 2.5 3]';
      F=[1 1 -1 -1 1 1 -1 -1 1 1 -1 -1]';
      % user data startup
   elseif 	nargin==2 & isnumeric(arg1) & isnumeric(arg2)
      t=arg1(:);	F=arg2(:);	% make columns
		if any(diff(t)<=0)
			error('all t must be strictly increasing')
		end
		if (length(t)~=length(F))
			error('t and F must be of equal length')
		end
	else		% user function startup
		t=arg1;		F=arg2;		% this F is a string
		funflag=1;
	end
	coefs=buttcoefs(5,1);
	q=struct('Nstage',coefs.Nstage, ...
		'mode',1, ...
		'stage',1, ...
		'order',coefs.order, ...
		'tau',coefs.tau, ...
		'K',coefs.K, ...
		'fn',coefs.fn, ...
		'zeta',coefs.zeta, ...
		'y0',zeros(1,coefs.Nstage), ...
		'yp0',zeros(1,coefs.Nstage), ...
		'A',1, ...
		'f',3, ...
		't',t, ...
		'F',F, ...
		'funflag',funflag, ...
		'Fi',1,'Ff',2, ...
		'dispax',0,'lhs',[0 0 0], ...
		'magax',0,'phaseax',0,'poleax',0,'popstage',0, ...
		'trivia',0,'consttext',0,'inittext',0, ...
		'genclump',[0;0],'stepclump',zeros(7,1), ...
		'harmclump',zeros(7,1),'desclump',zeros(10,1));

    a=figure('Units','normalized','Position',[0.17 0.07 0.8 0.8], ...
      'Numbertitle','off','Name','MultiSys','Tag','MultiSysBR549');
	set(0,'ShowHiddenHandles','on')
	wh=findobj('Tag','MultiSysBR549');
	s=length(wh);
	if (s>1)
		set(a,'Name',['MultiSys #',num2str(s)]);
	end
	set(0,'ShowHiddenHandles','off')
	currdir=pwd;
	homedir=pathto('mainswitch.m');
	cd(homedir)
	d=dir('fontpref.mat');
	if (length(d))
		load fontpref
		set(a,'DefaultAxesFontName',axfont.FontName, ...
			'DefaultAxesFontUnits',axfont.FontUnits, ...
			'DefaultAxesFontSize',axfont.FontSize, ...
			'DefaultAxesFontWeight',axfont.FontWeight, ...
			'DefaultAxesFontAngle',axfont.FontAngle, ...
			'DefaultUicontrolFontName',uifont.FontName, ...
			'DefaultUicontrolFontUnits',uifont.FontUnits, ...
			'DefaultUicontrolFontSize',uifont.FontSize, ...
			'DefaultUicontrolFontWeight',uifont.FontWeight, ...
			'DefaultUicontrolFontAngle',uifont.FontAngle)
	end
	cd(currdir)
	color=get(a,'Color');
	% display objects for general mode
	q.dispax=axes('Parent',a,'Units','normalized', ...
		'Position',[0.33 0.12 0.6 0.80],'box','on');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.94 0.195 0.05],'BackgroundColor',color, ...
	   'String','Window Mode');
	str='Line Segments';
	if funflag str=F;	end
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'Position',[0.015 0.88 0.195 0.05],'BackgroundColor',color, ...
	   'String',[str,'|Step Response|Steady Harmonic|Design'], ...
	   'Callback','multisys popup','Value',1);
	rbutton=uicontrol('Parent',a,'Style','pushbutton', ...
	   'Units','normalized','Position',[0.015 0.015 0.15 0.06], ...
	   'String','Return','Callback','multisys return');
	q.genclump=[q.dispax;rbutton];
	% display objects for step mode
	steptitl=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.767 0.195 0.05],'BackgroundColor',color, ...
	   'String','Step Input','Visible','off');
	Fiedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.707 0.15 0.06],'String',num2str(q.Fi), ...
	   'callback','multisys Fi','Visible','off');
	Filabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.707 0.03 0.05],'BackgroundColor',color, ...
	   'String','Fi','Visible','off');
	Ffedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.626 0.15 0.06],'String',num2str(q.Ff), ...
	   'callback','multisys Ff','Visible','off');
	Fflabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.626 0.03 0.05],'BackgroundColor',color, ...
	   'String','Ff','Visible','off');
	q.trivia=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0 0.192 0.21 0.2],'BackgroundColor',color, ...
	   'String',{'B = xxx';'M = xxx';'delay = zzz'},'Visible','off');
	q.stepclump=[q.dispax;steptitl;Fiedit;Filabel;Ffedit;Fflabel;q.trivia];
	% display objects for steady harmonic mode
	harmtitl=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.767 0.195 0.05],'BackgroundColor',color, ...
	   'String','Harmonic Input','Visible','off');
	Aedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.707 0.15 0.06],'String',num2str(q.A), ...
	   'callback','multisys A','Visible','off');
	Alabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.707 0.03 0.05],'BackgroundColor',color, ...
	   'String','A','Visible','off');
	fedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.626 0.15 0.06],'String',num2str(q.f), ...
	   'callback','multisys f','Visible','off');
	flabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.626 0.03 0.05],'BackgroundColor',color, ...
	   'String','f','Visible','off');
	q.harmclump=[q.dispax;harmtitl;Aedit;Alabel;fedit;flabel;q.trivia];
	% display objects for design mode
	q.magax=axes('Parent',a,'Units','normalized', ...
		'Position',[0.55 0.58 0.42 0.38],'box','on','Visible','off');
	q.phaseax=axes('Parent',a,'Units','normalized', ...
		'Position',[0.55 0.08 0.42 0.38],'box','on','Visible','off');
	q.poleax=axes('Parent',a,'Units','normalized', ...
		'Position',[0.07 0.08 0.40 0.38],'box','on','Visible','off');
	newb=uicontrol('Parent',a,'Style','pushbutton', ...
	   'Units','normalized','Position',[0.275 0.87 0.195 0.06], ...
	   'Visible','off','String','New Design','Callback','multisys new');
	editb=uicontrol('Parent',a,'Style','pushbutton','Visible','off', ...
	   'Units','normalized','Position',[0.275 0.789 0.195 0.06], ...
	   'String','Edit Stage','Callback','multisys edit');
	stagelabel=uicontrol('Parent',a,'Style','text', ...
	   'Units','normalized','Position',[0.015 0.789 0.1 0.06], ...
	   'BackgroundColor',color,'String','Stage','Visible','off');
	q.popstage=uicontrol('Parent',a,'Style','popupmenu', ...
	   'Units','normalized','Position',[0.115 0.789 0.095 0.06], ...
	   'BackgroundColor',color,'String',['1|2|3'],'Visible','off', ...
	   'Callback','multisys popstage','Value',1);
	stageframe=uicontrol('Parent',a,'Style','frame', ...
	   'Units','normalized','Position',[0.015 0.55 0.455 0.218], ...
	   'BackgroundColor',color,'Visible','off');
	q.consttext=uicontrol('Parent',a,'Style','text','Visible','off', ...
	   'Units','normalized','Position',[0.025 0.56 0.213 0.188], ...
	   'String',{'fn = x.xxxx';'zeta = x.xxxx';'K = 1'}, ...
	   'BackgroundColor',color,'HorizontalAlignment','left');
	q.inittext=uicontrol('Parent',a,'Style','text','Visible','off', ...
	   'Units','normalized','Position',[0.248 0.56 0.212 0.188], ...
	   'String',{'ORDER 2';'y0 = x.xxxx';'yp0 = x.xxxx'}, ...
	   'BackgroundColor',color,'HorizontalAlignment','left');
	q.desclump=[q.magax;q.phaseax;q.poleax;newb;editb;stagelabel;
		q.popstage;stageframe;q.consttext;q.inittext];
	
	b = uimenu('Parent',a,'Label','MultiSys','Tag','uimenu1');
	c1 = uimenu('Parent',b,'Label','Help','Enable','on', ...
		'Tag','helpmenu','Callback','multisys help');
	c2 = uimenu('Parent',b,'Label','Visible','Enable','on', ...
		'Callback',['set(gcf,''Handlevisibility'',''on'');', ...
			'set(get(gcbo,''UserData''),''Enable'',''on'');', ...
			'set(gcbo,''Enable'',''off'')'],'Tag','onmenu');
	c3 = uimenu('Parent',b,'Label','Callback','Enable','off', ...
		'Callback',['set(gcf,''Handlevisibility'',''callback'');', ...
			'set(get(gcbo,''UserData''),''Enable'',''on'');', ...
			'set(gcbo,''Enable'',''off'')'],'UserData',c2);
	set(c2,'UserData',c3)
% 	c=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
% 	   'Position',[0.015 0.09 0.05 0.06],'String','q', ...
% 	   'Callback','q=multisys(''q'')');
	% start program
	q.lhs=gendisp(q);
	set(gcf,'HandleVisibility','callback')
	set(a,'UserData',q)

else
	error('Argument Error!')
end

function lhs=gendisp(q)
axes(q.dispax),	cla,	hold on
[t,y]=solvesys(q,[min(q.t) max(q.t)]);
if q.funflag
	F=feval(q.F,t);
	title(['System Response to Forcing Function ',q.F,'(t)'])
else
	F=interp1(q.t,q.F,t);
	title('System Response to Line Segments Input')
end
Fh=plot(t,F,'k','LineWidth',1.5);
lhs=zeros(1,q.Nstage);
lstr=cell(1,q.Nstage+1);
lstr{1}='F';
co=get(q.dispax,'colororder');
lco=size(co,1);
for k=1:q.Nstage
	lhs(k)=plot(t,y(:,k),'color',co(mod(k-1,lco)+1,:));
	lstr{k+1}=int2str(k);
end
xlabel('time'),	ylabel('y')
legend([Fh,lhs],lstr{:})

function stepdisp(q)
axes(q.dispax),	cla,	hold on
T=zeros(q.Nstage,1);	% approx time scale for display limits
for k=1:q.Nstage
	if q.order(k)==1
		T(k)=q.tau(k);
	else
		if q.zeta(k)<=1
			T(k)=1/(q.zeta(k)*2*pi*q.fn(k));
		else
			T(k)=1/((q.zeta(k)-sqrt(q.zeta(k)^2-1))*2*pi*q.fn(k));
		end
	end
end
T=sum(T);
q2=q;	% set up modified q structure for step solution
q2.y0=q.Fi*cumprod(q.K);	% steady solution at Fi
q2.yp0=zeros(1,q.Nstage);
q2.t=[0 4*T];	q2.F=[1 1]*q.Ff;	q2.funflag=0;
[t,y]=solvesys(q2,q2.t);
Fh=plot([-T 0 0 4*T],[q.Fi q.Fi q.Ff q.Ff],'k','LineWidth',1.5);
lhs=zeros(1,q.Nstage);
lstr=cell(1,q.Nstage+1);
lstr{1}='F';
co=get(q.dispax,'colororder');
lco=size(co,1);
for k=1:q.Nstage
	lhs(k)=plot([-T;t],[q2.y0(k);y(:,k)],'color',co(mod(k-1,lco)+1,:));
	lstr{k+1}=int2str(k);
end
y=y(:,q.Nstage);
start=q.Fi*prod(q.K);	targ=q.Ff*prod(q.K);
index=1:201;	crit=abs(targ-start)/10;
rise=min(index(abs(y-targ)<crit))-1;
settle=max(index(abs(y-targ)>crit));
if (settle<rise) settle=rise; end
threshr=targ+crit*sign(start-targ);
riset=interp1(y(rise:rise+1),t(rise:rise+1),threshr);
threshs=targ+crit*sign(y(settle)-targ);
settlet=interp1(y(settle:settle+1),t(settle:settle+1),threshs);
xd=[[1 1]*riset NaN [1 1]*settlet];
yd=[threshr (start+targ)/2 NaN threshs (start+targ*2)/3];
plot(xd,yd,'k')
xd=[riset/2 (settlet+4*T)/2 NaN riset/2 (settlet+4*T)/2];
yd=[[1 1]*targ+crit NaN [1 1]*targ-crit];
plot(xd,yd,'k--')
th=get(gca,'Title');
font=get(th,'FontName');	fsz=get(th,'FontSize')+3;
th=text(riset,(start+targ)/2,'Rise');
set(th,'FontSize',fsz,'FontName',font,'HorizontalAlignment','right')
th=text(settlet,(start+targ*2)/3,'Settle');
set(th,'FontSize',fsz,'FontName',font,'HorizontalAlignment','left')
title('Multiple System Response to Step Input')
xlabel('time'),	ylabel('y')
legend([Fh,lhs],lstr{:})
res=200*10^floor(log10(4*T));
riset=round(riset*res)/res;
settlet=round(settlet*res)/res;
trivia={['Rise Time = ',num2str(riset)];['Settle Time = ',num2str(settlet)]};
set(q.trivia,'String',trivia)

function harmdisp(q)
w=2*pi*q.f;
Ms=ones(1,q.Nstage);	phis=Ms;	yhs=Ms;
for k=1:q.Nstage
	if q.order(k)==1
		Ms(k)=1/sqrt(1+(w*q.tau(k))^2);
		phis(k)=-atan(w*q.tau(k));
	else
		wn=2*pi*q.fn(k);
		Ms(k)=1/sqrt((1-(w/wn)^2)^2 + (2*q.zeta(k)*w/wn)^2);
		phis(k)=atan2(-2*q.zeta(k)*w/wn, 1-(w/wn)^2);
	end
end
delay=sum(phis)/w;
xb=1/q.f-delay;
t=[0:0.005:1]*xb;
F=q.A*sin(w*t);
axes(q.dispax),	cla,	hold on
co=get(q.dispax,'colororder');
lco=size(co,1);
lstr=cell(1,q.Nstage+1);
lstr{1}='F';
Fh=plot(t,F,'k','LineWidth',1.5);
for k=1:q.Nstage
	ys=q.A*prod(Ms(1:k).*q.K(1:k))*sin(w*t + sum(phis(1:k)));
	yhs(k)=plot(t,ys,'color',co(mod(k-1,lco)+1,:));
	lstr{k+1}=int2str(k);
end
B=q.A*prod(Ms.*q.K);	M=prod(Ms);	ord=sum(q.order);
h=max(abs([q.A B]))*sign(B);
xd=[-delay+[1 1.3]*1/4/q.f NaN, ...
	 [1 1]/2/q.f NaN [1 1]/2/q.f-delay 1/2/q.f+[-delay -0.1/q.f] NaN, ...
	 1/2/q.f-[1 0 1]*delay/3.5 NaN 1/2/q.f-delay+[1 0 1]*delay/3.5];
yd=[[1 1]*B NaN, ...
	 0 -0.8*h NaN 0 -0.8*h,[h h]*-0.6 NaN, ...
	 -0.6*h+[1 0 -1]*h/20 NaN -0.6*h+[1 0 -1]*h/20];
plot(xd,yd,'k')
th=get(gca,'Title');
font=get(th,'FontName');	fsz=get(th,'FontSize')+3;
th=text(-delay+1/4/q.f*1.35,B,'B');
set(th,'FontName',font,'FontSize',fsz)
th=text(0.4/q.f,-0.6*h,['\beta_',int2str(ord)],'fontsize',12);
set(th,'FontName',font,'FontSize',fsz,'horizontalalignment','right')
trivia={['B = ',num2str(B)];['M = ',num2str(M)];
	['delay = ',num2str(delay)]};
set(q.trivia,'String',trivia)
xlabel('time'),	ylabel('y')
title('Steady Response to Sinusoidal Input')
legend([Fh,yhs],lstr{:})

function desdisp(q)
axes(q.magax), hold off,	cla
co=get(q.dispax,'colororder');
lco=size(co,1);
lstr=cell(1,q.Nstage+1);
fc=zeros(q.Nstage,1);
for k=1:q.Nstage
	if q.order(k)==1
		fc(k)=1/(2*pi*q.tau(k));
	else
		fc(k)=q.fn(k);
	end
end
f=exp(linspace(log(min(fc)/10),log(max(fc)*10),201));
w=2*pi*f;
Ms=zeros(q.Nstage+1,201);	phis=Ms;	yhs=zeros(q.Nstage+1,1);
for k=1:q.Nstage
	if q.order(k)==1
		Ms(k,:)=1./sqrt(1+(w*q.tau(k)).^2);
		phis(k,:)=-180/pi*atan(w*q.tau(k));
	else
		wn=2*pi*q.fn(k);
		Ms(k,:)=1./sqrt((1-(w/wn).^2).^2 + (2*q.zeta(k)*w/wn).^2);
		phis(k,:)=180/pi*atan2(-2*q.zeta(k)*w/wn, 1-(w/wn).^2);
	end
end
Ms(q.Nstage+1,:)=prod(Ms(1:q.Nstage,:),1);
phis(q.Nstage+1,:)=sum(phis(1:q.Nstage,:),1);
for k=1:q.Nstage+1
	yhs(k)=loglog(f,Ms(k,:),'color',co(mod(k-1,lco)+1,:));
	lstr{k}=int2str(k);	hold on
end
lstr{q.Nstage+1}='Total';
xlabel('Frequency (Hz)'),	ylabel('Magnification')
title('Magnitude Response')
legend(yhs,lstr{:}),	grid on
axes(q.phaseax),	hold off,	cla
co=get(q.phaseax,'colororder');
lco=size(co,1);
for k=1:q.Nstage+1
	yhs(k)=semilogx(f,phis(k,:),'color',co(mod(k-1,lco)+1,:));	hold on
end
xlabel('Frequency (Hz)'),	ylabel('Phase Shift (Degrees)')
title('Phase Response')
legend(yhs,lstr{:}),	grid on
axes(q.poleax),	cla
px=zeros(1,sum(q.order));	py=px;
first=1;
for k=1:q.Nstage
	if q.order(k)==1
		px(first)=-1/(2*pi*q.tau(k));
		py(first)=0;
	else
		if q.zeta(k)<1			% under-damped
			px(first:first+1)=-q.fn(k)*q.zeta(k)*[1 1];
			py(first:first+1)=q.fn(k)*sqrt(1-q.zeta(k)^2)*[1 -1];
		elseif q.zeta(k)==1		% critical damped
			px(first:first+1)=-q.fn(k)*[1 1];
			py(first:first+1)=[0 0];
		else					% over-damped
			px(first:first+1)=-q.fn(k)*(1+sqrt(q.zeta(k)^2-1)*[1 -1]);
			py(first:first+1)=[0 0];
		end
	end
	first=first+q.order(k);
end
plot(px,py,'rx','markersize',10), grid on
drawnow
xl=get(q.poleax,'XLim');
if xl(2)<0	set(q.poleax,'XLim',[xl(1) -0.2*xl(1)]),	end
xlabel('Real (Hz)'),	ylabel('Imaginary Frequency (Hz)')	
title('System Pole Plot (S plane)')
destxt(q)

function destxt(q)
if q.order(q.stage)==1
	str={['tau = ',num2str(q.tau(q.stage))];
		['K = ',num2str(q.K(q.stage))]};
	set(q.consttext,'String',str)
	str={'ORDER 1';['y0 = ',num2str(q.y0(q.stage))]};
	set(q.inittext,'String',str)
else
	str={['fn = ',num2str(q.fn(q.stage))];
		['zeta = ',num2str(q.zeta(q.stage))];
		['K = ',num2str(q.K(q.stage))]};
	set(q.consttext,'String',str)
	str={'ORDER 2';
		['y0 = ',num2str(q.y0(q.stage))];
		['y''0 = ',num2str(q.yp0(q.stage))]};
	set(q.inittext,'String',str)
end

function [t,y]=solvesys(q,trange)
parms=struct('Nstage',q.Nstage,'order',q.order, ...
	'tau',q.tau,'K',q.K, ...
	'A',-4*pi*q.zeta.*q.fn, ...
	'B',-(2*pi*q.fn).^2, ...
	'C',(2*pi*q.fn).^2.*q.K, ...
	't',q.t,'F',q.F);
r0=zeros(sum(q.order),1);
first=1;		% set up initial condition array
for k=1:q.Nstage
	r0(first)=q.y0(k);
	if q.order(k)==2	r0(first+1)=q.yp0(k);	end
	first=first+q.order(k);
end
[t,r]=ode45('multisys',trange,r0,[],parms,q.funflag);
y=zeros(length(t),q.Nstage);
first=1;
for k=1:q.Nstage
	y(:,k)=r(:,first);
	first=first+q.order(k);
end

% design coupled dynamic systems implementing
% a Butterworth LPF of order N
function coefs=buttcoefs(N,fc)
if mod(N,2)		% odd-order filter
	rkr=ones(1,(N-1)/2);
	k=1:(N-1)/2;
	rkr(k+1)=cos(pi*k/N);	% real parts of the system poles
	coefs=struct('Nstage',(N+1)/2, ...
		'order',[1 2*ones(1,(N-1)/2)], ...
		'tau',[1/(2*pi*fc) ones(1,(N-1)/2)], ...
		'K',ones(1,(N+1)/2), ...
		'fn',[1 fc*ones(1,(N-1)/2)], ...
		'zeta',rkr);
else			% even-order filter
	l=N/2;
	rkr=cos(pi/N*(-0.5 + 1:l));
	coefs=struct('Ntage',N/2, ...
		'order',2*ones(1,N/2), ...
		'tau',ones(1,N/2), ...
		'K',ones(1,N/2), ...
		'fn',fc*ones(1,N/2), ...
		'zeta',rkr);
end

function dest=pathto(fname)
r=which(fname);
[tok1,r]=strtok(r,filesep);
dest=tok1;
[tok2,r]=strtok(r,filesep);
while ~isempty(r)
	tok1=tok2;
	dest=[dest,filesep,tok1];
	[tok2,r]=strtok(r,filesep);
end

% by Tim Conover 2/2/2000
