function plot_trajectory(s0,T)

if (nargin < 2) | isempty(T),  T = 500;  end

global energy_contour

oldcwd = pwd;
cd /mit/bert/Thesis/matlab/poincare;

[times,states] = ode45('d_delfi_f', [0;T], s0, ...
    odeset('AbsTol',[1e-4 1e-4 1e-4 1e-4], 'RelTol',1e-5));

delta = states - (ones(length(states),1) * s0(:)');
d = sum((delta.^2)');

plot(states(:,1),states(:,2));

set_contour(delfi(s0(1),s0(2)) + (s0(3).^2 + s0(4).^2)./2);
line(energy_contour(1,:),energy_contour(2,:), 'Color', [0 0 1]);
set(gca, 'DataAspectRatio',[1 1 1], 'Box','on');

title(sprintf('period = ~ %.2f', times( find_minimum(d) )));

cd(oldcwd);



function [index] = find_minimum(data)

% find indices with low values
low = find(data < .1);
% eliminate contiguous indices near the beginning
low = low( find(low - low(1) + 1 - (1:length(low))) );

if isempty(low), index=1; return; end

% extract next contiguous indices
band = low( find( ~(low - low(1) + 1 - (1:length(low))) ) );

% find the lowest value; return first occurrence
index = find(data == min( data(band) ));
index = index(1);
disp(sprintf('residual = %g', data(index)));