%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% df is a tool for plotting a phase plane and solutions of a
%% first order differential system of the form
%%
%%     x' = f(x,y)
%%     y' = g(x,y)
%%
%%
%%    df( F, x0, x1, y0, y1, vals)
%%
%%    F: a function handle describing the right-hand side [f(x,y),g(x,y)]
%%  x0,x1: minimum and maximum values of the x coordinate for the plot window
%%  y0,y1: minimum and maximum values of the x coordinate for the plot window
%%   vals: A (n x 2) matrix of initial data for solutions to be plotted, of the form
%%         [ X0 , Y0; X1 , Y1; X2, Y2, ... ]
%%

function pp( F, x0, x1, y0, y1, vals)
  % Number of grid points for the direction field
	Nx=21; Ny=21;

	% Create the grid for the direction field
	x=linspace(x0,x1,Nx);
	y=linspace(y0,y1,Ny);
	[xx,yy]=meshgrid(x,y);

  % The vectors to plot at each grid point	
        s1 = @(x) x(1);
        s2 = @(x) x(2);
        f = @(x,y) s1(F(0,x,y));
        g = @(x,y) s2(F(0,x,y));
        
        dx = arrayfun(f,xx,yy);
        dy = arrayfun(g,xx,yy);
	
	Lx = max(max(abs(dx)));
	Ly = max(max(abs(dy)));
        DX = abs(x1-x0)/Nx; DY = abs(y1-y0)/Ny;
        lam = max(DX/Lx,DY/Ly);
	
  % DT = abs(t1-t0)/Nt; DU = abs(u1-u0)/Nu;
  % lam = 3*max(abs(dt)/DT,abs(du)/DU);
	dt = dx.*lam;
	du = dy.*lam;
	
	% Undo any scaling that Octave/Matlab might prefer
  S = 0;

  % Quiver plots vectors with their base at each vertex.  To get centered line segments,
  % use two vectors, each with a base at the vertex, pointing in opposite directions.
  % Plot the first half of each line segment
	h=quiver(xx,yy,dt,du,S);
	set(h,'linewidth',1.4)
        % set(h,'showarrowhead', 'off');
	
        % hold on % (keep drawing on the same figure)
        % % Plot the second half of each line segment
        % h=quiver(tt,uu,-dt,-du,S);    
        % set(h,'showarrowhead', 'off');
	

        % % For each initial data, plot the curve corresponding to that solution
        % if( exist('vals') )
        %       nv=size(vals,1);        
        %       for(k=1:nv)
        %               tv=vals(k,1);
        %               uv=vals(k,2);
        %               dt=(t1-t0)/NT;
        % 
        %               % Compute the solution curve forward in time
        %               [U,T,istate,msg] = odewalk( ode, tv, uv, dt, t0, t1, u0, u1);
        %               h=plot(T,U,'color','black');
        %               set(h,'linewidth',1.5);
        % 
        %               bwode = @(t,u) -ode(-t,u); % ODE for running backwards in time
        %               [U,T,istate,msg] = odewalk( bwode, -tv, uv, dt, -t1, -t0, u0, u1 );
        %               h=plot(-T,U,'color','black');
        %               set(h,'linewidth',1.5);
        % 
        %         plot(tv,uv,'rs')
        %       end
        % end
	% Constrain the plot window
	axis([x0,x1,y0,y1]);
	% Stop plotting everything in the same window
	hold off
end

function [U,T,istate,msg]=odewalk(ode, tv, uv, dt, t0, t1, u0, u1)

	% Detect if we are running octave or matlab
  Octave = exist('OCTAVE_VERSION') ~= 0;

	% Call the appropriate code
	if( Octave )
		[U,T,istate,msg]=odewalkOctave(ode, tv, uv, dt, t0, t1, u0, u1);
	else
		[U,T,istate,msg]=odewalkMatlab(ode, tv, uv, dt, t0, t1, u0, u1);
	end
end

function [U,T,istate,msg]=odewalkOctave(ode, tv, uv, dt, t0, t1, u0, u1)
	U=[uv];
	T=[tv];
	istate = 0;
	msg = 'initial value out of bounds';

	% Octave uses a different convention for argument order!
	odeOctave = @(u,t) ode(t,u);

	% Step the solution forward until we leave the view window
	% or encounter an error from lsode
%  N=2*(t1-t0)/dt;
  k = 0;
	while( (tv<t1) )
    k = k + 1;
		[UU,istate,msg]=lsode(odeOctave,uv,[tv,tv+dt]);
		if(istate==2)
			U=[U, UU(2)];
			T=[T, tv+dt];
			if( UU(2)<u0 || UU(2)>u1)
				break
			end
	  else
			break
		end
	  tv = tv +dt;
	  uv = UU(2);
	end
end

function [U,T,istate,msg]=odewalkMatlab(ode, tv, uv, dt, t0, t1, u0, u1)
	U=[uv];
	T=[tv];
	% Emulate lsode return codes and pretend like we always succeed
	istate = 2;
	msg = '';
	% Step the solution forward until we leave the view window
	% or encounter an error from lsode
	
	while( tv<t1)
		[TT,UU] = ode45(odeMatlab,[tv,tv+dt],uv);
		U=[U, UU(end)];
		T=[T, TT(end)];
		if( UU(end)<u0 || UU(end)>u1)
				break
		end
		tv = tv+dt;
		uv = UU(end);
	end
end
