function [] = explore_poincare_map(name, E, xlim, ylim)

% exhaustive_map(SD, V, E)... 
%    this really needs a blurb!

global pmArgs           % arguments from previous Poincare' plot
global pmFigure 	% handles for figures for the Poincare' plot
global pmCurrent	% which figure is current?
global pmLimits		% current axis limits for q2, p2 (yes, in that order)
global pmDefaultLimits	% default axis limits for q2, p2 (yes, in that order)
global m                % particle mass

nx = 5;  ny = 5;
T = 600;
m = 1;
%format short g

if isunix
  tmpdir = '/tmp/';
else
  tmpdir = '\windows\temp\';
end

%%%% initial setup

V = name;
sysder  = [ 'd_' name '_f' ];
reverse = [ 'd_' name '_b' ];

% V = fcnchk(V, 'x', 'y');
plots = setup_for_map(sysder, V, E);

% set up axes
ax = { gcfa(pmFigure(1)), gcfa(pmFigure(2)) };

% ODE integrator options
opts = odeset('Events','on', 'AbsTol',[1e-7 1e-7 1e-7 1e-7], 'RelTol',1e-8);

state = [];		% current state vector (empty => none)

col_bgd  = { [ 0  0  0], [ 0  0  0] };	% plot background
col_axis = { [ 1  1  1], [ 1  1  1] };	% plot axes
col_new  = { [ 1  1  1], [ 1  1  1] };	% most recent points
col_old  = { [.5 .5  1], [ 1 .5 .5] };	% old points
col_edge = { [.6 .4  0], [.5 .5 .5] };	% energy boundary

done = 0;		% 1 if we want to end the loop

if isempty(pmLimits)
  pmLimits = pmDefaultLimits;
end

% reset axis limits
figure(pmFigure(pmCurrent));  axes(ax{pmCurrent});
axis(pmLimits{pmCurrent});

% reset colors
cx = 1 + (pmCurrent ~= 1);
set(ax{1}, 'Color', col_bgd{cx}, ...
    'XColor', col_axis{cx}, 'YColor', col_axis{cx});
set(plots{1}(1),        'Color',     col_new{cx});
set(plots{1}(2),        'Color',     col_old{cx});
set(plots{1}(3),        'EdgeColor', col_edge{cx});

cx = 1 + (pmCurrent ~= 2);
set(ax{2}, 'Color', col_bgd{cx}, ...
    'XColor', col_axis{cx}, 'YColor', col_axis{cx});
set(plots{2}(1),        'Color',     col_new{cx});
set(plots{2}(2),        'Color',     col_old{cx});
set(plots{2}(3),        'EdgeColor', col_edge{cx});

%%% set plot limits

if (nargin < 3) | isempty(xlim)
  sx = pmLimits{pmCurrent}(2)-pmLimits{pmCurrent}(1);
  xlim = pmLimits{pmCurrent}(1):(sx/nx):pmLimits{pmCurrent}(2);
end
if (nargin < 4) | isempty(ylim)
  sy = pmLimits{pmCurrent}(4)-pmLimits{pmCurrent}(3);
  ylim = pmLimits{pmCurrent}(3):(sy/ny):pmLimits{pmCurrent}(4);
end

%%%% main loop


for gX = xlim,
  for gY = ylim,

    % clear other relevant data
    times = [];  states = [];  Ptimes = [];  Pstates = [];  state = [];

    % set coordinates (compute the state vector)
    if (pmCurrent == 1)
      p1sqr = 2*m*(E - feval(V, 0, gX)) - gY.^2;
      % vector will have q2,p2 from click coords; q1=0, p1 calculated from E.
      if (p1sqr >= 0), state = [0, gX, sqrt(p1sqr), gY]; end
    else
      p2sqr = 2*m*(E - feval(V, gX, 0)) - gY.^2;
      % vector will have q1,p1 from click coords; q2=0, p2 calculated from E.
      if (p2sqr >= 0), state = [gX, 0, gY, sqrt(p2sqr)]; end
    end

    if isempty(state),
      disp('Outside of allowed energy region!');
    else
      disp(sprintf('[%f %f %f %f]', state));
      istate = state;

      % have starting state: plot!

      % crunch data
      disp('integrating...');
      [times,states,Ptimes,Pstates,Pwhy] = ode45(sysder, [0;T], state, opts);

      % add data to plots
      disp('plotting...');
      if ~isempty(Pstates)
	Xs = find(Pwhy == 1);
	Ys = find(Pwhy == 2);
	if ~isempty(Xs)
	  set(plots{1}(1), ...
	      'XData', [ get(plots{1}(1), 'XData'), Pstates(Xs,2)' ], ...
	      'YData', [ get(plots{1}(1), 'YData'), Pstates(Xs,4)' ]);
	end
	if ~isempty(Ys)
	  set(plots{2}(1), ...
	      'XData', [ get(plots{2}(1), 'XData'), Pstates(Ys,1)' ], ...
	      'YData', [ get(plots{2}(1), 'YData'), Pstates(Ys,3)' ...
	      ]);
	end
	drawnow;  % force graphics update
    end

      % set new start state to state at the end of integration
      if isempty(Pstates)
	sidx = size(states);
	state = states(sidx(1), :);
      else
	sidx = size(Pstates);
	state = Pstates(sidx(1), :);
	% regurgitate energy
	if 0
	disp(sprintf('[%f %f %f %f]  E0=%f', state, ...
	    feval(V, state(1), state(2)) + (state(3).^2 + state(4).^2)/(2*m)));
	if (Pwhy(sidx(1)) == 1)
	  state = [0, state(2), ...
		sqrt(2*m*(E-feval(V,0,state(2))) - state(4).^2), state(4)];
	else
	  state = [state(1), 0, ...
		state(3), sqrt(2*m*(E-feval(V,state(1),0)) - state(3).^2)];
	end
	end
      end

      disp(sprintf('[%f %f %f %f]   E=%f', state, ...
	  feval(V, state(1), state(2)) + (state(3).^2 + state(4).^2)/(2*m)));
    end
  end
end

hold off
disp('Bye bye.');

%%%% helper functions

function [axis_handle] = gcfa(fig)
% like gca, but for a specified figure
saved = gcf;
figure(fig);
axis_handle = gca;
figure(saved);

function [axis_tag] = xytag()
% return 'x' or 'y' depending on pmCurrent
global pmCurrent;
axis_tag = char(122-pmCurrent);

function [exch] = other(curr)
% select the other plot
exch = 3 - curr;
