function map = poincare(plane, xyz, derivative, derivative_wrt)
%
% Computes a Poincare map of a trajectory based on an arbitrary plane
%
% arguments:
%  plane: the three points that define the plane, i.e. 
%   [[x1 y1 z1];[x2 y2 z2];[x3 y3 z3]]
%  xyz: trajectory (3 dimensions)
%  derivative: either 'positive' or 'negative'
%  derivative_wrt: derivative with respect to what variable (1=x, 2=y, 3=z)
%
% return:
%  map: a 3xm matrix containing intersection points, i.e. the Poincare map
%

% get the length of the time series
%
l = length(xyz(:,1));

% initialize the poincare map
%
map=[];

% loop over time series
%
for i=1:l-1
   
    % determine the sign of the derivative between these two points
    %
    delta = xyz(i,derivative_wrt)-xyz(i+1,derivative_wrt);
    
    % only continue if sign agrees with 'derivative' argument
    %
    if (((delta < 0) & (derivative == 'negative')) |...
        ((delta > 0) & (derivative == 'positive')))
        
        % determine whether the line formed by the two points intersects
        % the plane defined by the user
        %
        line = [xyz(i,:);xyz(i+1,:)];

        den = [[1 1 1 0];
               [plane(1,1),plane(2,1),plane(3,1),line(2,1)-line(1,1)];
               [plane(1,2),plane(2,2),plane(3,2),line(2,2)-line(1,2)];
               [plane(1,3),plane(2,3),plane(3,3),line(2,3)-line(1,3)]];
   
        num = [[1 1 1 1];
               [plane(1,1),plane(2,1),plane(3,1),line(1,1)];
               [plane(1,2),plane(2,2),plane(3,2),line(1,2)];
               [plane(1,3),plane(2,3),plane(3,3),line(1,3)]];

        t = -det(num)/det(den);

        x = line(1,1) + (line(2,1)-line(1,1))*t;
        y = line(1,2) + (line(2,2)-line(1,2))*t;
        z = line(1,3) + (line(2,3)-line(1,3))*t;

        % the point at which the line intersects the user defined plane
        %
        point = [x,y,z];
        
        % determine whether or not this point of intersection lies between
        % the current two points
        %
        if (((derivative == 'negative') &... 
             (point(derivative_wrt)-xyz(i,derivative_wrt) > 0) &... 
             (point(derivative_wrt)-xyz(i+1,derivative_wrt) < 0)) |...
            ((derivative == 'positive') &...
             (point(derivative_wrt)-xyz(i,derivative_wrt) < 0) &... 
             (point(derivative_wrt)-xyz(i+1,derivative_wrt) > 0)))
         
            % if so, add this point to the Poincare map
            %
            map = [map;[x,y,z]];
        end
    end
end

