function [] = explore_poincare_map(name, E)

% interactive_map(name, E) takes the name of a two-dimensional potential,
%    and plots the Poincare' section for a fixed energy E.  It allows
%    interactive exploration of the behavior of the system via keyboard and
%    mouse.  For NAME to be a valid argument to interactive_map, NAME,
%    d_NAME_f and d_NAME_b must all be valid functions.  (NAME is the
%    potential function [which is assumed to depend only on coordinates],
%    and the other two are forward-going and backward-going system
%    derivatives.)
%
%    Since we're working with a two-dimensional system, the state
%    vector of the system is assumed to be [q1, q2, p1, p2], where q1
%    and q2 are generalized coordinates and p1 and p2 their respecive
%    conjugate momenta.
%
%    [BACKGROUND: The 'system derivative' is a function returning the time
%    derivative of the system state vector.  For a Hamiltonian system, the
%    behavior is specified by the Hamilton's equations:
%        Dq/Dt =  dH/dp
%        Dp/Dt = -dH/dq
%    or, written out into components for a 2D system,
%        Dq1/Dt =  dH/dp1
%        Dq2/Dt =  dH/dp2
%        Dp1/Dt = -dH/dq1
%        Dp2/Dt = -dH/dq2
%    where D_/D_ indicates a total derivative and d_/d_ a partial one.  The
%    system derivative function will return the column vector made up of
%    the four partial derivatives on the RHS.]
%
%    Given the required assumptions, the "forwards" system derivative
%    function should return a column vector looking something like this:
%
%        [state(3)/m ; ...
%         state(4)/m ; ...
%         - << expression for dV/dq1 >> ; ...
%         - << expression for dV/dq2 >>]
%
%    Furthermore, the system derivative must support the 'event' option
%    (described in ODEFILE documentation), and use it to extract crossings
%    of qi==0 where pi>0, perhaps by doing something like this:
%
%      case 'events'   % crossings info; used only if odeset('Events','on').
%        % look at intersections [ x=0 y=0 ]...
%        out1 = [state(1) state(2)];
%        % ...where [ x y ] is increasing, ie [ px>0 py>0 ].
%        out3 = [      1        1 ];  
%        % 0=don't stop, 1=stop
%        out2 = [      0        0 ];
%
%    The "backwards" system derivative should just reverse the signs of the
%    system derivative and out3 (the intersection direction).  The only
%    reason for its existence is that (1) I couldn't get odeXX to do the
%    right thing by itself, and (2) calculating it using the "forwards"
%    function would involve clunky, inefficient kludges.
%
%    This function always attempts to continue an existing plot, so
%    you may need to use "clf" to start a new one.  The figure is
%    cleared automatically when switching evergy levels.

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 pmLimitsStack	% a list of saved limits views
global m                % particle mass

%direction = +1;
no_tweak_energy = 0;
T = 100;
N = 1;
m = 1;
%format short g

if isunix,
  printdir = '/tmp/';
  rawdir   = '/afs/zone.mit.edu/user/bert/Thesis-poincare-raw/';
else
  printdir = '\windows\temp\';
  rawdir   = '\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)) };

%%%% spew instructions

disp(sprintf('\nSorry, help is unavailable until I update it.\n'));

% disp(sprintf('\nmapping/evolution'));
% disp(          '  button1: set current position');
% disp(          '  button2: show trajectory');
% disp(          '  button3: plot time N*T forward from pointer position');
% disp(          '  e:       show phase space evolution');
% disp(          '  [,]:     increase/decrease time interval');
% disp(          '  D:       delete current data set (ie, white points)');
% disp(sprintf('\nviews'));
% disp(          'backspace: switch active plot');
% disp(          '  z:       zoom in 2x around current position');
% disp(          '  Z:       zoom out to the entire image');
% disp(          '  x:       zoom in 2x in X coord around the image center');
% disp(          '  X:       zoom out 2x in X coord around the image center');
% disp(sprintf('\nrandom other functions'));
% disp(          '  p:       print to file, landscape');
% disp(          '  P:       print to file, tall/swapped');
% disp(          '  S:       save to file, binary format');
% disp(sprintf(  '  q:       quit\n'));

%%%% loop initializations

% stepQ = .1;  stepP = .4;  stepDIV = 1;

% 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)

lockX=[];		% if nonempty, this specifies a fixed X coordinate
lockY=[];		% if nonempty, this specifies a fixed Y coordinate

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

%%%% main loop

% loop until the user says to stop
while ~done  
  % 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});

  % get a key press or button click
  [gX, gY, btn] = ginput(1);

  % adjust them if necessary
  if ~isempty(lockX)
    gX=lockX;
  end
  if ~isempty(lockY)
    gX=lockY;
  end

  if isempty(btn)
    disp('No button?  I think I am confused.');

  elseif (btn == 1)
    % button 1 pressed

    % A. move plots(1) plot data to plots(2)
    set(plots{1}(2), ...
	'XData', [ get(plots{1}(2), 'XData'), get(plots{1}(1), 'XData') ], ...
	'YData', [ get(plots{1}(2), 'YData'), get(plots{1}(1), 'YData') ]);
    set(plots{1}(1), 'XData',[], 'YData',[]);    
    set(plots{2}(2), ...
	'XData', [ get(plots{2}(2), 'XData'), get(plots{2}(1), 'XData') ], ...
	'YData', [ get(plots{2}(2), 'YData'), get(plots{2}(1), 'YData') ]);
    set(plots{2}(1), 'XData',[], 'YData',[]);    

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

    % C. 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;
    end

  elseif ((btn == 3)|(btn == 13))
    % button 3 or Enter pressed

    if ~isempty(state)
      % have starting state: plot!

      % repeat in a loop (for some reason helps conserve memory)
      for i = 1:N,

	% 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
	end

	% set new start state to state at the end of integration
	if no_tweak_energy | isempty(Pstates),
	  sidx = size(states);
	  state = states(sidx(1), :);
	else
	  sidx = size(Pstates);
	  state = Pstates(sidx(1), :);
	  % adjust "invisible" momentum to retain energy conservation
	  %disp(sprintf('[%f %f %f %f]  E0=%f  %d/%d', state, ...
	  %    feval(V, state(1), state(2)) ...
	  %    + (state(3).^2 + state(4).^2)/(2*m), ...
	  %    i, N));
	  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

	disp(sprintf('[%f %f %f %f]   E=%f  %d/%d', state, ...
	    feval(V, state(1), state(2)) ...
	    + (state(3).^2 + state(4).^2)/(2*m), ...
	    i, N));
	drawnow;
      
      end % of the 1:N loop

      disp('done.');

    else
      % no starting state
      disp('I need a starting point, Dave.');
    end

  elseif (btn == 2)
    % button 2 pressed
    if exist('states') & ~isempty(states)
      figure;
      plot(states(:,1), states(:,2));
      xlabel x;  ylabel y;  set(gca, 'DataAspectRatio', [1 1 1]);
      text(0, -30, sprintf('E=%f     s_0=[%f %f %f %f]', E, istate), ...
	  'Units', 'Points', ...
	  'FontSize', 7.5, 'FontUnits', 'Points', ...
	  'VerticalAlignment', 'baseline', 'HorizontalAlignment', 'left' );

    else
      disp('I need some data before I can do that... try button3.');
    end

  elseif (btn == 'e') | (btn == 'E')
    % button 2 pressed
    if ~isempty(states)
      figure;
      plot(states(:,1), states(:,2));
      xlabel x;  ylabel y;

      figure;
      plot(states(:,3), states(:,4));
      xlabel p_x;  ylabel p_y;

      figure;
      plot(states(:,1), states(:,3));
      xlabel x;  ylabel p_x;

      figure;
      plot(states(:,2), states(:,4));
      xlabel y;  ylabel p_y;

    else
      disp('I need some data before I can do that... try button3.');
    end

  elseif (btn == 'D')
    % remove the current plot
    state = [];
    set(plots{1}(1), 'XData',[], 'YData',[]);    
    set(plots{2}(1), 'XData',[], 'YData',[]);    
    times = [];  states = [];  Ptimes = [];  Pstates = [];
    disp('Where was the kaboom? There was supposed to be a giant plot-shattering kaboom!');

  elseif (btn == '-')
    % swap plot style: '-' <-> '.'
    if strcmp(get(plots{pmCurrent}(1), 'LineStyle'), '-'),
      set(plots{pmCurrent}(1), 'LineStyle', 'none');
    else
      set(plots{pmCurrent}(1), 'LineStyle', '-');
    end
  elseif (btn == 19) % ^S
    % sparsify grid by a factor of 2 (to do: factor changes?)
    tmpX = get(plots{pmCurrent}(1), 'XData');
    tmpY = get(plots{pmCurrent}(1), 'YData');
    set(plots{pmCurrent}(1), 'XData', tmpX(1:2:length(tmpX)), ...
	                     'YData', tmpY(1:2:length(tmpY)));
    tmpX = [];
    tmpY = [];
    disp('Yum, tasty.');

  elseif (btn == '[')
    if (N > 1)
      N = N./2;
      disp(sprintf('N -> %d   (N*T -> %f)', N, N*T));
    else
      disp('N cannot be less than 1!');
    end
  elseif (btn == ']')
    N = N*2;
    disp(sprintf('N -> %d   (N*T -> %f)', N, N*T));
  elseif (btn == '{')
    T = T./2;
    disp(sprintf('T -> %f   (N*T -> %f)', T, N*T));
  elseif (btn == '}')
    T = T*2;
    disp(sprintf('T -> %f   (N*T -> %f)', T, N*T));
  elseif (btn == 8) | (btn == 127)   % backspace (^H) or delete
    % flip direction
    [sysder, reverse] = swap(sysder, reverse);
    disp(sysder);

  elseif (btn == 'p')
    fname = [printdir sprintf('poin%d', E*100) xytag '.ps' ];
    disp([ 'printing [landscape] to ' fname ]);
    set(ax{pmCurrent}, 'Color', [1 1 1], 'XColor', [0 0 0], 'YColor', [0 0 0]);
    set(plots{pmCurrent}(1), 'Color', [0 0 0]);
    set(plots{pmCurrent}(2), 'Color', [0 0 0]);
    set(plots{pmCurrent}(3), 'EdgeColor', [0 0 0]);
    orient landscape;  print('-deps', fname);
  
  elseif (btn == 'P')
    fname = [printdir sprintf('poin%d', E*100) xytag 't.ps' ];
    disp([ 'printing [swapped, tall] to ' fname ]);

    % swap X, Y
    set(plots{pmCurrent}(1), 'XData', get(plots{pmCurrent}(1), 'YData'), ...
	                     'YData', get(plots{pmCurrent}(1), 'XData'));
    set(plots{pmCurrent}(2), 'XData', get(plots{pmCurrent}(2), 'YData'), ...
                             'YData', get(plots{pmCurrent}(2), 'XData'));
    set(plots{pmCurrent}(3), 'XData', get(plots{pmCurrent}(3), 'YData'), ...
                             'YData', get(plots{pmCurrent}(3), 'XData'));
    axis(pmLimits{pmCurrent}([3 4 1 2]));
    xlabel([ 'p_' xytag ]); ylabel(xytag);

    set(ax{pmCurrent}, 'Color', [1 1 1], 'XColor', [0 0 0], 'YColor', [0 0 0]);
    set(plots{pmCurrent}(1), 'Color', [0 0 0]);
    set(plots{pmCurrent}(2), 'Color', [0 0 0]);
    set(plots{pmCurrent}(3), 'EdgeColor', [0 0 0]);
    orient tall;  print('-deps', fname);

    % unswap X, Y
    set(plots{pmCurrent}(1), 'XData', get(plots{pmCurrent}(1), 'YData'), ...
	                     'YData', get(plots{pmCurrent}(1), 'XData'));
    set(plots{pmCurrent}(2), 'XData', get(plots{pmCurrent}(2), 'YData'), ...
                             'YData', get(plots{pmCurrent}(2), 'XData'));
    set(plots{pmCurrent}(3), 'XData', get(plots{pmCurrent}(3), 'YData'), ...
                             'YData', get(plots{pmCurrent}(3), 'XData'));
    axis(pmLimits{pmCurrent});
    xlabel(xytag); ylabel([ 'p_' xytag ]);

  elseif (btn == 's')
    fname = sprintf('%sf%d_%d', rawdir, pmCurrent, E*100);
    disp([ 'saving to ' fname ' (.m/.mat)' ]);
    print('-dmfile', fname);
  elseif (btn == 'S')
    figure(pmFigure(1));
    fname = sprintf('%sf1_%d', rawdir, E*100);
    disp([ 'saving to ' fname ' (.m/.mat)' ]);
    print('-dmfile', fname);
    figure(pmFigure(2));
    fname = sprintf('%sf2_%d', rawdir, E*100);
    disp([ 'saving to ' fname ' (.m/.mat)' ]);
    print('-dmfile', fname);

  elseif (btn == 'z')
    % zoom in 2x around cursor
    wid = (pmLimits{pmCurrent}(2) - pmLimits{pmCurrent}(1))/4;
    hgt = (pmLimits{pmCurrent}(4) - pmLimits{pmCurrent}(3))/4;
    pmLimits{pmCurrent} = [ gX-wid gX+wid gY-hgt gY+hgt ];
  
  elseif (btn == 'Z')
    % zoom out 2x around cursor [Q: should we make it around center? stack?]
    wid = (pmLimits{pmCurrent}(2) - pmLimits{pmCurrent}(1));
    hgt = (pmLimits{pmCurrent}(4) - pmLimits{pmCurrent}(3));
    pmLimits{pmCurrent} = [ gX-wid gX+wid gY-hgt gY+hgt ];

  elseif (btn == 26)	% ^Z
    % reset pmLimits to defaults
    pmLimits{pmCurrent} = pmDefaultLimits{pmCurrent};

  elseif (btn == 'x')
    pmLimits{pmCurrent} = ...
	[(3*pmLimits{pmCurrent}(1)+pmLimits{pmCurrent}(2))/4 ...
	 (3*pmLimits{pmCurrent}(2)+pmLimits{pmCurrent}(1))/4 ...
	 pmLimits{pmCurrent}(3) ...
	 pmLimits{pmCurrent}(4)];

  elseif (btn == 'X')
    pmLimits{pmCurrent} = ...
	[(3*pmLimits{pmCurrent}(1)-pmLimits{pmCurrent}(2))/2 ...
	 (3*pmLimits{pmCurrent}(2)-pmLimits{pmCurrent}(1))/2 ...
	 pmLimits{pmCurrent}(3) ...
	 pmLimits{pmCurrent}(4)];

  elseif (btn == 'v')
    pmLimitsStack{pmCurrent} = ...
	cat(2, pmLimitsStack{pmCurrent}, {pmLimits{pmCurrent}});
    disp('view pushed');

  elseif (btn == 'V')
    if isempty(pmLimitsStack{pmCurrent}),
      disp([ 'no saved ' xytag '/p' xytag ' views.' ]);
    else
      pmLimits{pmCurrent} = pmLimitsStack{pmCurrent}{1};
      pmLimitsStack{pmCurrent} = ...
	  pmLimitsStack{pmCurrent}(2:length(pmLimitsStack{pmCurrent}));
      disp('view popped');
    end

  elseif (btn == 9)   % TAB
    pmCurrent = other(pmCurrent);

  elseif ((btn == 'q') | (btn == 'Q'))
    done = 1;

  else
    disp(sprintf('button/key %d is not supported', btn));
%    [beep]
  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.  (hack!)
global pmCurrent;
axis_tag = char(122-pmCurrent);

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

function [two, one] = swap(one, two)
% swap 2 args (syntactic sugar, doesn't use a variable)
