function ReturnVal=ODEPlot(fn,arg)

%    ODEPlot: Differential equation plotter w/Direction and Vector fields
%    Copyright (C) 1996  Brad North

%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.

%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

%    Problems, Comments, and Suggestions can be sent to 
%            northbh@nextwork.rose-hulman.edu

if nargin == 0
	fn = 'init';
end;

Figure_Name='ODE Plotter';

%Activates the figure and makes sure it has been properly initialized

if ~strcmp(fn,'init')
	h_MainFig = gcf;

%If the current figure is not named Figure_Name then find 
%the figure that is and put its handle in h_MainFig. This will 
%be found by checking the children of root.

	if ~strcmp(get(h_MainFig,'name'),Figure_Name)
		h_figs=get(0,'children');
		h_MainFig = findobj(h_figs,'flat',...
				    'name',Figure_Name);

	%If the figure does not exist, initialize it and then
	%run the function that was originally intended

		if length(h_MainFig)==0
			odeplot('init');
			odeplot(fn);
			return;
		end;
	end;
end;
%____________________________________________________________________________

%Function: PlotDField

if strcmp(fn,'PlotDField')

%This function reads in minx, maxx, step, mint, and maxt
%which are the min and max x and t values and the step values.


	h_data=get(h_MainFig,'User');
	h_tMin=h_data(2);
	h_tMax=h_data(3);
	h_xMin=h_data(4);
	h_xMax=h_data(5);
	h_Step=h_data(6);
	h_EqnString=h_data(7);
	h_axis=h_data(12);

	mint=str2num(get(h_tMin,'string'));
	maxt=str2num(get(h_tMax,'string'));

	minx=str2num(get(h_xMin,'string'));
	maxx=str2num(get(h_xMax,'string'));

	step=str2num(get(h_Step,'string'));

	T=mint:step:maxt;
	stepx=(maxx-minx)/(maxt-mint)*step;
	X=minx:stepx:maxx;

	[t,x]=meshgrid(T,X);

	dx=eval(get(h_EqnString,'string'));
	ind=find(imag(dx)~=0);
	dx(ind)=NaN*ones(size(ind));

	m=dx(:)';
	theta=atan(m.*stepx./step);
	Arrows=[zeros(1,length(theta));...
		step*cos(theta)+j*stepx*sin(theta);...
		NaN*ones(1,length(theta));...
		.9*step*cos(theta-pi/6)+j*.8*stepx*sin(theta-pi/24);...
		step*cos(theta)+j*stepx*sin(theta);...
		.9*step*cos(theta+pi/6)+j*.8*stepx*sin(theta+pi/24)];
	Grid=t-j*x;
	Grid=Grid(:)';
	PlotData=.75.*Arrows+ones(6,1)*Grid;
	PlotSize=size(PlotData);
	PlotData=[PlotData; ones(1,PlotSize(2))*NaN];
	PlotData=PlotData(:);
	axes(h_axis);

	set(h_data(10),'xdata',real(PlotData),...
		       'ydata',imag(PlotData),...
		       'color',[1 1 1]);
	
	xlabel('t');
	ylabel('x');
	title(['dx/dt =' get(h_EqnString,'string')]);
	set(h_axis,'xlim',[mint maxt],...
		   'ylim',[minx maxx]);
	set(h_MainFig,'User',h_data);

end;
%____________________________________________________________________________

%Function: VectorField

if strcmp(fn,'VectorField');

	h_data=get(h_MainFig,'User');
	h_tMin=h_data(2);
	h_tMax=h_data(3);
	h_xMin=h_data(4);
	h_xMax=h_data(5);
	h_Step=h_data(6);
	h_EqnString=h_data(7);
	h_axis=h_data(12);

	mint=str2num(get(h_tMin,'string'));
	maxt=str2num(get(h_tMax,'string'));

	minx=str2num(get(h_xMin,'string'));
	maxx=str2num(get(h_xMax,'string'));

	step=str2num(get(h_Step,'string'));

	T=mint:step:maxt;
	X=minx:step:maxx;

	[t,x]=meshgrid(T,X);

	dx=eval(get(h_EqnString,'string'));

	arrow=[0;1;NaN;.8+j*.08;1;.8-j*.08];
	Slopes=dx(:)';
	ind=find(imag(Slopes)~=0);
	Scalar=.9*step/(max(abs(Slopes)));
	Slopes(ind)=NaN*ones(size(ind));

	Grid=t-j*x;
	Grid=Grid(:)';
	PlotData=Scalar*arrow*Slopes+ones(6,1)*Grid;
	PlotSize=size(PlotData);
	PlotData=[PlotData; ones(1,PlotSize(2))*NaN];
	PlotData=PlotData(:);
	axes(h_axis);
	h_data(10)=plot(real(PlotData),imag(PlotData),'w');

	xlabel('t');
	ylabel('x');
	title(['dx/dt =' get(h_EqnString,'string')]);
	set(h_axis,'xlim',[mint maxt],...
		   'ylim',[minx maxx]);
	set(h_MainFig,'user',h_data);
end;
%____________________________________________________________________________

%Function: SetEqnLim

if strcmp(fn,'SetEqnLim')

	h_data=get(h_MainFig,'User');
	
	h_DataFig=h_data(1);
	Handle=findobj('name','First Order Differential Equation');
	if Handle~=h_DataFig

		odeplot('SetupDataFig');
		h_data=get(h_MainFig,'User');
		set(h_data(1),'visible','on');

	else
		set(h_DataFig,'Visible','on');

	end;
	
end;
%____________________________________________________________________________

%Function: PlotButton

if strcmp(fn,'PlotButton')

	h_data=get(h_MainFig,'User');

	h_DataFig=h_data(1);
	h_DirectionType=h_data(8);

	set(h_DataFig,'visible','off');
	Checked=get(h_DirectionType,'Checked');

	if strcmp(Checked,'on')

		odeplot('PlotDField');

	else

		odeplot('VectorField');

	end;
	

end;

%____________________________________________________________________________

%Function: SetupDataFig

if strcmp(fn,'SetupDataFig')

	h_DataFig=figure('menu','none',...
			 'name','First Order Differential Equation',...
			 'number','off',...
			 'pos',[100 100 300 150],...
			 'visible','off',...
			 'color',[.73 .73 .73]); 

	%UIControls for h_DataFig

	uicontrol(h_DataFig,...
	'pos',[.1 .75 .8 .1],...
	'back',[.73 .73 .73],...
	'string','Enter differential equation using x and t',...
	'style','text',...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.01 .6 .14 .1],...
	'back',[.73 .73 .73],...
	'string','dx/dt =',...
	'style','text',...
	'units','normal');

	h_EqnString=uicontrol(h_DataFig,...
		    'pos',[.17 .6 .78 .1],...
		    'back',[1 1 1],...
		    'style','edit',...
		    'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.16 .58 .8 .13],...
	'back',[1 1 1],...
	'fore',[.5 .5 .5],...
	'style','frame',...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.16 .4 .1 .1],...
	'string','t min',...
	'style','text',...
	'units','normal',...
	'back',[.73 .73 .73]);

	uicontrol(h_DataFig,...
	'pos',[.28 .38 .14 .14],...
	'style','frame',...
	'fore',[.5 .5 .5],...
	'back',[1 1 1],...
	'units','normal');

	h_tMin=uicontrol(h_DataFig,...
	'pos',[.3 .4 .1 .1],...
	'style','edit',...
	'back',[1 1 1],...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.5 .4 .1 .1],...
	'string','t max',...
	'style','text',...
	'units','normal',...
	'back',[.73 .73 .73]);

	uicontrol(h_DataFig,...
	'pos',[.62 .38 .14 .14],...
	'style','frame',...
	'fore',[.5 .5 .5],...
	'back',[1 1 1],...
	'units','normal');

	h_tMax=uicontrol(h_DataFig,...
	'pos',[.64 .4 .1 .1],...
	'style','edit',...
	'back',[1 1 1],...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.16 .22 .1 .1],...
	'string','x min',...
	'style','text',...
	'units','normal',...
	'back',[.73 .73 .73]);

	uicontrol(h_DataFig,...
	'pos',[.28 .2 .14 .14],...
	'style','frame',...
	'fore',[.5 .5 .5],...
	'back',[1 1 1],...
	'units','normal');

	h_xMin=uicontrol(h_DataFig,...
	'pos',[.3 .22 .1 .1],...
	'style','edit',...
	'back',[1 1 1],...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.5 .22 .1 .1],...
	'string','x max',...
	'style','text',...
	'units','normal',...
	'back',[.73 .73 .73]);

	uicontrol(h_DataFig,...
	'pos',[.62 .2 .14 .14],...
	'style','frame',...
	'fore',[.5 .5 .5],...
	'back',[1 1 1],...
	'units','normal');

	h_xMax=uicontrol(h_DataFig,...
	'pos',[.64 .22 .1 .1],...
	'style','edit',...
	'back',[1 1 1],...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.3 .06 .18 .1],...
	'style','text',...
	'back',[.73 .73 .73],...
	'units','normal',...
	'string','Step Size');

	h_Step=uicontrol(h_DataFig,...
	'pos',[.51 .06 .1 .1],...
	'style','edit',...
	'back',[1 1 1],...
	'units','normal');

	uicontrol(h_DataFig,...
	'pos',[.49 .04 .14 .14],...
	'style','frame',...
	'back',[1 1 1],...
	'fore',[.5 .5 .5],...
	'units','normal');
	
	uicontrol(h_DataFig,...
	'pos',[.78 .1 .2 .16],...
	'units','normal',....
	'string','Plot',...
	'callback','odeplot(''PlotButton'');');


	h_data=get(h_MainFig,'User');

	h_data(1)=h_DataFig;
	h_data(2)=h_tMin;
	h_data(3)=h_tMax;
	h_data(4)=h_xMin;
	h_data(5)=h_xMax;
	h_data(6)=h_Step;
	h_data(7)=h_EqnString;

	set(h_MainFig,'User',h_data);
end;
%____________________________________________________________________________

%Function: SwitchChecks

if strcmp(fn,'SwitchChecks')

	h_data=get(h_MainFig,'user');
	h_DirectionType=h_data(8);
	h_VectorType=h_data(9);
	h_HideGraph=h_data(11);
	
	if strcmp(arg,'Direction')

		set(h_DirectionType,'check','on');
		set(h_VectorType,'Check','off');
		Checked=get(h_HideGraph,'Check');
		set(h_HideGraph,'label','Hide Direction Field');
		set(h_HideGraph,'check',Checked);
	else

		set(h_DirectionType,'check','off');
		set(h_VectorType,'Check','on');
		Checked=get(h_HideGraph,'Check');
		set(h_HideGraph,'label','Hide Vector Field');
		set(h_HideGraph,'check',Checked);

	end;


end;

%____________________________________________________________________________

%Function: HideGraph

if strcmp(fn,'HideGraph')

	h_data=get(h_MainFig,'User');
	h_VPlot=h_data(10);
	h_HideMenu=h_data(11);
	h_axis=h_data(12);

	xlim=get(h_axis,'xlim');
	ylim=get(h_axis,'ylim');
	
	Checked=get(h_HideMenu,'checked');

	if strcmp(Checked,'on')


		set(h_HideMenu,'checked','off');
		set(h_VPlot,'visible','on');

	else

		set(h_HideMenu,'checked','on');
		set(h_VPlot,'visible','off');
	end;
	
	set(h_axis,'xlim',xlim,...
		   'ylim',ylim);

end;

%____________________________________________________________________________

%Function: PlotCurve

if strcmp(fn,'PlotCurve')

	h_data=get(h_MainFig,'User');
	h_tmin=h_data(2);
	h_tmax=h_data(3);
	h_EqnString=h_data(7);

	set(h_MainFig,'pointer','watch');
	Header='function dx=odetemp(t,x)\ndx=';
	Function=[get(h_EqnString,'string') ';'];
	fid=fopen('c:\odetemp.m','w');
	path(path,'c:\');
	fprintf(fid,[Header Function]);
	fclose(fid);
	t0=1;
	x0=1;
	tf=str2num(get(h_tmax,'string'));
	[t,x]=ode23('odetemp',t0,tf,x0);
	set(h_MainFig,'pointer','arrow');
	figure(h_MainFig);
	axes(h_data(12));
	hold on;
	plot(t,x)
	!erase c:\odetemp.m
	end;



%____________________________________________________________________________

%Function: init

if strcmp(fn,'init')

	delete(findobj('name',Figure_Name));
	delete(findobj('name','First Order Differential Equation'));
	
	h_MainFig=figure('menu','none',...
			 'name',Figure_Name,...
			 'number','off',...
			 'pos',[30 30 500 350],...
			 'PaperPosition',[.25 .5 8 10]);

	h_axis=axes;

	h_RootFile=uimenu(h_MainFig,...
			  'label','&File');


	uimenu(h_RootFile,...
	       'label','Printer Se&tup',...
	       'CallBack','print -dsetup');

	uimenu(h_RootFile,...
	       'label','&Print',...
	       'CallBack','Print');

	h_RootEdit=uimenu(h_MainFig,...
			  'label','&Edit');

	uimenu(h_RootEdit,...
	       'label','&Copy',...
	       'CallBack','print -dbitmap');

	h_RootGraph=uimenu(h_MainFig,...
			   'label','&Graph');

	uimenu(h_RootGraph,...
	      'label','&First Order DE',...
	      'callback','odeplot(''SetEqnLim'');')

	h_HideGraph=uimenu(h_RootGraph,...
	      'label','&Hide Direction Field',...
	      'CallBack','odeplot(''HideGraph'');');


	h_GraphType=uimenu(h_RootGraph,...
			   'label','&Type',...
			   'separator','on');

	h_DirectionType=uimenu(h_GraphType,...
	      'label','&Direction Field',...
	      'Checked','on',...
	      'CallBack','odeplot(''SwitchChecks'',''Direction'');');

	h_VectorType=uimenu(h_GraphType,...
	      'label','&Vector Field',...
	      'CallBack','odeplot(''SwitchChecks'',''Vector'');');


	uimenu(h_RootGraph,...
		'separator','on',...
		'label','Plot Integral Curve',...
		'CallBack','odeplot(''PlotCurve'');');
%---------------------------------------------------------
	odeplot('SetupDataFig');
%---------------------------------------------------------

	h_data=get(h_MainFig,'user');
	figure(h_MainFig);
	h_data(8)=h_DirectionType;
	h_data(9)=h_VectorType;
	h_data(10)=plot(0,0,'k');     %Reserved for the vector field
	h_data(11)=h_HideGraph;
	h_data(12)=h_axis;

	set(h_MainFig,'User',h_data);
end;
