function out = FirstOrd(arg1,arg2,flag,parms)
% FIRSTORD  illustrates first order system responses, including
% step response of a lumped thermal parameter model, Example 3.3,
% and sinusoidal forcing function response, Example 3.6
% FIRSTORD(t,F) uses columns of time data t and force data F(t), 
% to define a forcing function made up of line segments.
% In this case the response is evaluated numerically using the 
% Matlab function 0DE45() instead of using an analytic solution.

if nargin==4	% ODEFile functionality
	t=arg1;	y=arg2;		% (flag is not used)
	f=interp1(parms.t,parms.F,t);
	out=(parms.K*f-y)/parms.tau;
	
elseif (nargin==1 & isstr(arg1))	% GUI callback functionality
	action=arg1;
	p=get(gcbf,'UserData');
	switch (action)
	case 'popup'
		switch p.mode
		case 'step'
			set(p.stepclump,'Visible','off')
			set(p.tauedit,'String',num2str(p.tau))
		case 'ramp'
			set(p.rampclump,'Visible','off')
		case 'general'
			set(p.genclump,'Visible','off')
		case 'start'
			set(p.startclump,'Visible','off')
			delete(p.dimline),	delete(p.text1)
		case 'steady'
			set(p.steadyclump,'Visible','off')
			set(p.trivia,'Position',[0.015 0.03 0.15 0.16])
			delete(p.dimline),	delete(p.text1),	delete(p.text2)
		end
		val=get(gcbo,'Value');
		modes={'step';'ramp';'general';'start';'steady'};
		p.mode=modes{val};
		switch p.mode
		case 'step'
			set(p.stepclump,'Visible','on')
			p.tau=p.m*p.c/p.h/p.As;
			xlabel('time (sec)'),	ylabel('Temperature (\circC)')
			title('Response of Thermal System to Step Input')
			doit(p)
		case 'ramp'
			set(p.rampclump,'Visible','on')
			xlabel('t'),	ylabel('y')
			title(['Response of First-Order System to F=at+b'])
			doit(p)
		case 'general'
			set(p.genclump,'Visible','on')
			xlabel('t'),	ylabel('y')
			title(['Response of First-Order System to Line Segments F(t)'])
			doit(p)
		case 'start'
			set(p.startclump,'Visible','on')
			xlabel('t'),	ylabel('y')
			title(['Startup of First-Order System with ', ...
				'F=A sin(2\pift)'])
			p.dimline=line('Parent',gca,'XData',0,'YData',p.y0, ...
				'Color','k','Tag','dimline');
			th=get(gca,'Title');
			font=get(th,'FontName');	fsz=get(th,'FontSize')+3;
			p.text1=text('Parent',gca,'Position',[0 p.y0],'String','', ...
				'FontName',font,'FontSize',fsz);
			doit(p)
		case 'steady'
			set(p.steadyclump,'Visible','on')
			set(p.trivia,'Position',[0.015 0.03 0.21 0.16])
			xlabel('t'),	ylabel('y')
			title(['Steady Response of First-Order System to ', ...
				'F=A sin(2\pift)'])
			p.dimline=line('Parent',gca,'XData',0,'YData',p.y0, ...
				'Color','k','Tag','dimline');
			th=get(gca,'Title');
			font=get(th,'FontName');	fsz=get(th,'FontSize')+3;
			p.text1=text('Parent',gca,'Position',[0 0],'String','', ...
				'FontName',font,'FontSize',fsz);
			p.text2=text('Parent',gca,'Position',[0 0],'String','', ...
				'FontName',font,'FontSize',fsz, ...
				'HorizontalAlignment','right');
			doit(p)
		end
	case {'m','c','h','As'}
	   eval(['p.',action,'=str2num(get(gcbo,''String''));'])
	   p.tau=p.m*p.c/p.h/p.As;
	   doit(p)
	case {'Ti','Tf','tau','K','y0','a','b','A','f'}
	   eval(['p.',action,'=str2num(get(gcbo,''String''));'])
	   doit(p)	
	case 'Response'
		str=sprintf('%0.15f',gcbo);	% relocate window, or build a new one
		bfh=findobj('tag',str);		% the button handle becomes a unique
		if ishandle(bfh)			% tag to identify the response figure
			figure(bfh)
		else
			bfh=figure('tag',str);
			mph=uimenu('Parent',bfh,'Label','disown');
			uimenu('Parent',mph,'Label','disown', ...
				'Callback',['set(gcbf,''Tag'',''''),' ...
					'ans=[0;0];ans(1)=get(gcbo,''Parent'');' ...
					'ans(2)=gcbo;delete(ans)']);
		end
		fc=1/(2*pi*p.tau);		% center plots around cutoff frequency
		f=[exp(linspace(log(fc/10),log(fc*10),201)),p.f];
		M=1./sqrt(1+(2*pi*f*p.tau).^2);
		phi=-180/pi*atan(2*pi*f*p.tau);
		subplot(2,1,1)
		loglog(f(1:201),M(1:201),'b',f(202),M(202),'ko')
		xlabel('Frequency'),ylabel('Magnitude Ratio M')
		title('Frequency Response')
		set(gca,'XLim',[min(f) max(f)])
		subplot(2,1,2)
		semilogx(f(1:201),phi(1:201),'b',f(202),phi(202),'ok')
		xlabel('Frequency'),ylabel('Phase Shift \phi (\circ)')
		set(gca,'XLim',[min(f) max(f)],'YLim',[-90 0])
	case 'return'
		choice=questdlg('Select Output Time Basis','Numeric Output', ...
			'Display','Resample','Display');
		switch (choice)
		case 'Display'	% get output already calculated for line using
			t=get(p.yline,'XData');	% ode45 internal time basis
			y=get(p.yline,'YData');
		case 'Resample'	% recalculate output using new sample rate
			T=p.t(end)-p.t(1);
			fs=(length(p.t)-1)/T;
			fstr=inputdlg('Enter Sampling Frequency','fs',1,{num2str(fs)});
			fs=str2num(fstr{1});
			trange=p.t(1):1/fs:p.t(end);
			parms=struct('tau',p.tau,'K',p.K,'t',p.t,'F',p.F);
			[t,y]=ode45('firstord',trange,p.y0,[],parms);
		end
		assignin('base','tout',t(:))
		assignin('base','yout',y(:))
		disp('Returned tout and yout to workspace')
	case 'p'
		out=p;
	case 'help'
    	helpwin('aboutfirstord')
	otherwise
		disp(['callback ',action,' unknown'])
	end
	set(gcbf,'UserData',p)
	
elseif (nargin==0 | (nargin==2 & isnumeric(arg1) & isnumeric(arg2)))
	if nargin==0					% easy invocation
		t=[-5:49]';
		F=[ones(1,10),2*sin(1:25),12-(1:20)]';
		mode='step';	popval=1;
	else						% user data invocation
		t=arg1(:);	F=arg2(:);	% make sure they are 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
		mode='general';	popval=3;
	end
	p=struct('m',8.855e-4, ...
	   'c',0.1393, ...
	   'h',0.245, ...
	   'As',7.854e-5, ...
	   'Ti',20, ...
	   'Tf',37, ...
	   'tau',6.4104, ...
	   'K',1, ...
	   'y0',2, ...
	   'a',1, ...
	   'b',-1, ...
	   'A',1, ...
	   'f',0.1, ...
	   'mode',mode, ...
	   't',t, ...
	   'F',F, ...
	   'Fline',0,'yline',0,'asympline',0,'dimline',0, ...
	   'text1',0,'text2',0,'trivia',0,'tauedit',0, ...
	   'stepclump',zeros(13,1), ...
	   'rampclump',zeros(10,1), ...
	   'genclump',zeros(7,1), ...
	   'startclump',zeros(11,1), ...
	   'steadyclump',zeros(9,1));
    
    a=figure('Units','normalized','Position',[0.2 0.2 0.75 0.65], ...
      'Numbertitle','off','Name','FirstOrd','Tag','FirstOrdBR549');
	set(0,'ShowHiddenHandles','on')
	wh=findobj('Tag','FirstOrdBR549');
	s=length(wh);
	if (s>1)
		set(a,'Name',['FirstOrd #',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');
	b=axes('Parent',a,'Units','normalized','XGrid','on','YGrid','on', ...
		'Position',[0.33 0.12 0.6 0.80],'Box','on');
   c=line('Parent',b,'XData',[],'YData',[],'Color','k','Tag','Fline', ...
      'LineWidth',1.5);
	p.Fline=c;
   c=line('Parent',b,'XData',[],'YData',[],'Color','b','Tag','yline', ...
      'LineWidth',1.5);
	p.yline=c;
	c=line('Parent',b,'XData',[],'YData',[],'Color','b', ...
		'LineStyle',':','Tag','asympline');
	p.asympline=c;
	% dimline will be created only when needed
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.94 0.195 0.05],'BackgroundColor',color, ...
	   'String','Input Function');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'Position',[0.015 0.88 0.195 0.05],'BackgroundColor',color, ...
	   'String',['Thermal Step|Ramp|Line Segments|', ...
	   		'Start Harmonic|Steady Harmonic'], ...
	   'Callback','firstord popup','Value',popval);
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.767 0.195 0.05],'BackgroundColor',color, ...
	   'String','Parameters');
	% edit boxes for step follow
	medit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.707 0.15 0.06],'String',num2str(p.m), ...
	   'callback','firstord m');
	mlabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.707 0.03 0.05],'BackgroundColor',color, ...
	   'String','m');
	cedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.610 0.15 0.06],'String',num2str(p.c), ...
	   'callback','firstord c');
	clabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.610 0.03 0.05],'BackgroundColor',color, ...
	   'String','c');
	hedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.513 0.15 0.06],'String',num2str(p.h), ...
	   'callback','firstord h');
	hlabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.513 0.03 0.05],'BackgroundColor',color, ...
	   'String','h');
	Asedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.415 0.15 0.06],'String',num2str(p.As), ...
	   'callback','firstord As');
	Aslabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.415 0.03 0.05],'BackgroundColor',color, ...
	   'String','As');
	Tiedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.318 0.15 0.06],'String',num2str(p.Ti), ...
	   'callback','firstord Ti');
	Tilabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.318 0.03 0.05],'BackgroundColor',color, ...
	   'String','Ti');
	Tfedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.221 0.15 0.06],'String',num2str(p.Tf), ...
	   'callback','firstord Tf');
	Tflabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.221 0.03 0.05],'BackgroundColor',color, ...
	   'String','Tf');
	trivia=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.015 0.03 0.15 0.16],'BackgroundColor',color, ...
	   'String',{['tau = ',num2str(p.tau)];'K=1'});
	p.trivia=trivia;
	p.stepclump=[medit;mlabel;cedit;clabel;hedit;hlabel;
			Asedit;Aslabel;Tiedit;Tilabel;Tfedit;Tflabel;trivia];
	% edit boxes for ramp follow
	tauedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.707 0.15 0.06],'String',num2str(p.tau), ...
	   'callback','firstord tau','Visible','off');
	p.tauedit=tauedit;
	taulabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.707 0.04 0.05],'BackgroundColor',color, ...
	   'String','tau','Visible','off');
	Kedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.610 0.15 0.06],'String',num2str(p.K), ...
	   'callback','firstord K','Visible','off');
	Klabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.610 0.03 0.05],'BackgroundColor',color, ...
	   'String','K','Visible','off');
	y0edit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.513 0.15 0.06],'String',num2str(p.y0), ...
	   'callback','firstord y0','Visible','off');
	y0label=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.513 0.03 0.05],'BackgroundColor',color, ...
	   'String','y0','Visible','off');
	aedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.415 0.15 0.06],'String',num2str(p.a), ...
	   'callback','firstord a','Visible','off');
	alabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.415 0.03 0.05],'BackgroundColor',color, ...
	   'String','a','Visible','off');
	bedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.318 0.15 0.06],'String',num2str(p.b), ...
	   'callback','firstord b','Visible','off');
	blabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.318 0.03 0.05],'BackgroundColor',color, ...
	   'String','b','Visible','off');
	p.rampclump=[tauedit;taulabel;Kedit;Klabel;y0edit;y0label;
		aedit;alabel;bedit;blabel];
	% return button for general mode
	rb=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
	   'Position',[0.015 0.015 0.15 0.06],'Visible','off', ...
	   'String','Return','Callback','firstord return');
	p.genclump=[tauedit;taulabel;Kedit;Klabel;y0edit;y0label;rb];
	% edit boxes for sinusoidal modes
	Aedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.415 0.15 0.06],'String',num2str(p.A), ...
	   'callback','firstord A','Visible','off');
	Alabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.415 0.03 0.05],'BackgroundColor',color, ...
	   'String','A','Visible','off');
	fedit=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.318 0.15 0.06],'String',num2str(p.f), ...
	   'callback','firstord f','Visible','off');
	flabel=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'Position',[0.18 0.318 0.03 0.05],'BackgroundColor',color, ...
	   'String','f','Visible','off');
	frb=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
		'Position',[0.015 0.221 0.15 0.06],'Visible','off', ...
		'String','Freq. Response','Callback','firstord Response');
	p.startclump=[tauedit;taulabel;Kedit;Klabel;y0edit;y0label;
		Aedit;Alabel;fedit;flabel;trivia];
	p.steadyclump=[tauedit;taulabel;Kedit;Klabel;
		Aedit;Alabel;fedit;flabel;frb;trivia];
	b = uimenu('Parent',a,'Label','Firstord','Tag','uimenu1');
	c1 = uimenu('Parent',b,'Label','Help','Enable','on','Tag','helpmenu', ...
		'Callback','firstord help');
	c2 = uimenu('Parent',b,'Label','Visible','Enable','on','Tag','onmenu', ...
		'Callback',['set(gcf,''Handlevisibility'',''on'');', ...
			'set(get(gcbo,''UserData''),''Enable'',''on'');', ...
			'set(gcbo,''Enable'',''off'')']);
	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.18 0.015 0.05 0.06],'String','p', ...
% 	   'Callback','p=firstord(''p'')');
	% start program
	switch (p.mode)
	case 'step'	
		xlabel('time (sec)'),	ylabel('Temperature (øC)')
		title('Response of Thermal System to Step Input')
		doit(p)
	case 'general'
		set(p.stepclump,'Visible','off')
		set(p.genclump,'Visible','on')
		xlabel('t'),	ylabel('y')
		title('Response of First-Order System to Line Segments F(t)')
		doit(p)
	end
	set(gcf,'HandleVisibility','callback')
	set(a,'UserData',p)

else
	error('Argument Error!')
end


function doit(p)
switch p.mode
case 'step'
	t=-p.tau:p.tau*6/100:5*p.tau;
	T=ones(1,101)*p.Ti;
	T(t>0)=p.Tf-(p.Tf-p.Ti)*exp(-t(t>0)/p.tau);
	set(p.Fline,'XData',[-p.tau 0 0 5*p.tau], ...
		'YData',[p.Ti p.Ti p.Tf p.Tf])	% input
	set(p.yline,'XData',t,'YData',T)	% output
	set(p.asympline,'XData',[0 p.tau],'YData',[p.Ti p.Tf])
	set(p.trivia,'String',{['tau = ',num2str(p.tau)];['K = ',num2str(p.K)]})
	legend([p.Fline,p.yline,p.asympline], ...
		'Environment','Response','initial slope',4)
case 'ramp'
	A=p.K*p.a;
	B=p.K*(p.b-p.a*p.tau);
	C=p.y0-B;
	t=[0:p.tau*5/100:5*p.tau]';
	y=A*t+B+C*exp(-t/p.tau);
	set(p.Fline,'XData',[0 5*p.tau],'YData',p.a*[0 5*p.tau]+p.b)
	set(p.yline,'XData',t,'YData',y)
	set(p.asympline,'XData',[0 5*p.tau],'YData',A*[0 5*p.tau]+B)
	legend([p.Fline,p.yline,p.asympline], ...
		'F','Response','asymptote',4)
case 'general'
	set(p.Fline,'XData',p.t,'YData',p.F)
	parms=struct('tau',p.tau,'K',p.K,'t',p.t,'F',p.F);
	trange=[p.t(1) p.t(end)];	% plot the line using ode45 internal 
	[t,y]=ode45('firstord',trange,p.y0,[],parms);	% time basis
	set(p.yline,'XData',t,'YData',y)
	set(p.asympline,'XData',p.t(1),'YData',p.F(1))
	legend([p.Fline,p.yline],'F','Response',1)
case 'start'
	w=2*pi*p.f;
	B=p.K*p.A/sqrt(1+(w*p.tau)^2);
	C=p.y0+p.K*p.A*w*p.tau/(1+(w*p.tau)^2);
	phi=-atan(w*p.tau);
	xb=1/p.f-phi/w;
	t=[0:1/200:1]*max([5*p.tau, xb]);
	F=p.A*sin(w*t);
	y=C*exp(-t/p.tau)+B*sin(w*t+phi);
	ys=B*sin(w*t+phi);
	set(p.Fline,'XData',t,'YData',F)
	set(p.yline,'XData',t,'YData',y)
	set(p.asympline,'XData',t,'YData',ys)
	xd=[0 1.1 NaN 0 1.1 NaN 1 1 NaN, ...
		0.95 1 1.05 NaN 0.95 1 1.05]*p.tau;
	y0=p.y0;
	yd=[y0 y0 NaN [1 1]*(y0-C) NaN y0-C y0 NaN, ...
		y0-C/20 y0 y0-C/20 NaN y0-C*0.95 y0-C y0-C*0.95];
	set(p.dimline,'XData',xd,'YData',yd);
	set(p.text1,'Position',[p.tau*1.1 (p.y0*2-C)/2],'String','C');
	legend([p.Fline,p.yline,p.asympline], ...
		'F','Response','asymptote',1)
	set(p.trivia,'String',['C = ',num2str(C)])
case 'steady'
	w=2*pi*p.f;
	B=p.K*p.A/sqrt(1+(w*p.tau)^2);
	phi=-atan(w*p.tau);
	delay=phi/w;
	xb=1/p.f-delay;
	t=[0:1/200:1]*xb;
% 	t=[0:1/200:1]*max([5*p.tau, xb]);
	F=p.A*sin(w*t);
	ys=B*sin(w*t+phi);
	set(p.Fline,'XData',t,'YData',F)
	set(p.yline,'XData',t,'YData',ys)
	set(p.asympline,'XData',0,'YData',0)
	h=max(abs([p.A B]))*sign(B);
% 	xd=[1/4/p.f-delay+[0 p.tau/3.5] NaN, ...
	xd=[-delay+[1 1.3]*1/4/p.f NaN, ...
		 [1 1]/2/p.f NaN [1 1]/2/p.f-delay 1/2/p.f+[-delay -0.1/p.f] NaN, ...
		 1/2/p.f-[1 0 1]*delay/3.5 NaN 1/2/p.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];
	set(p.dimline,'XData',xd,'YData',yd);
	set(p.text1,'Position',[-delay+1/4/p.f*1.35 B],'String','B')
	set(p.text2,'String','\beta_1', ...
		'Position',[0.4/p.f -0.6*h])
	legend([p.Fline,p.yline],'F','Response',1)
	set(p.trivia,'String',{['B = ',num2str(B)];['M = ',num2str(B/p.K/p.A)];['delay = ',num2str(delay)]})
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 10/3/99
