%% Henon–Heiles: Poincaré section (x, p_x) at y=0 with p_y>0
% H = 1/2(px^2+py^2+x^2+y^2) + x^2*y - (1/3)*y^3
% Section: y(t)=0 crossings with p_y(t)>0, plot (x, p_x).

clear; clc;
clf;

%% ---- User settings ----
E       = 1/12 + 0.03;        % total energy (e.g. 0.125; below escape is ~1/6)
Nsec    = 5000;       % number of section points to record
Tmax    = 2e5;        % max integration time (increase if needed)

% Choose one-parameter family of initial conditions on the section:
% We start exactly on the section: y(0)=0, p_y(0)>0
% Pick x(0) values, set p_x(0)=0, and determine p_y(0) from energy.
x_list  = linspace(0, 0.29, 16);   % several trajectories (adjust range)
px0     = 0.0;                      % common initial p_x
y0      = 0.0;

%% ---- Integrator options ----
opts = odeset( ...
    'RelTol',1e-10,'AbsTol',1e-12, ...
    'Events',@event_y0_pypos);

figure(1); hold on; box on;
xlabel('x','FontSize',14); ylabel('p_x','FontSize',14);
set(gca,'FontSize',14,'LineWidth',2);
title(sprintf('Henon–Heiles Poincare section: y=0, p_y>0,  E=%.6g',E));

%% ---- Loop over initial conditions ----
for k = 1:numel(x_list)
    x0 = x_list(k);

    % Compute p_y(0) from the energy constraint at y=0:
    % E = 1/2(px^2+py^2+x^2)  (since y=0 term and x^2*y - y^3/3 vanish)
    rad = 2*E - (px0^2 + x0^2);
    if rad <= 0
        warning('Skipping x0=%.4g: not allowed for this E with px0=%.4g',x0,px0);
        continue;
    end
    py0 = sqrt(rad);  % choose positive branch to start with p_y>0

    z0 = [x0; y0; px0; py0];

    % Storage for section points
    P = zeros(Nsec,2);
    count = 0;

    t0 = 0;
    zstart = z0;

    while count < Nsec && t0 < Tmax
        % Integrate until next event (y=0 with p_y>0)
        [t,z,te,ze,ie] = ode45(@henon_heiles, [t0 Tmax], zstart, opts);

        if isempty(te)
            break; % no more events found before Tmax
        end

        % Each event corresponds to y=0; we requested direction +1 and p_y>0
        % Record (x, p_x) at the event
        for j = 1:size(ze,1)
            count = count + 1;
            P(count,:) = [ze(j,1), ze(j,3)]; % x, p_x
            if count >= Nsec, break; end
        end

        % Restart just after the last event to avoid detecting the same one again
        t0 = te(end) + 1e-9;
        zstart = ze(end,:).';
    end

    P = P(1:count,:);
    plot(P(:,1), P(:,2), '.', 'MarkerSize', 3);
    drawnow
end

axis tight; grid on; hold off;

%% ======== Dynamics ========
function dz = henon_heiles(~, z)
    % z = [x; y; px; py]
    x  = z(1);  y  = z(2);
    px = z(3);  py = z(4);

    % Hamilton's equations:
    % xdot = dH/dpx = px
    % ydot = dH/dpy = py
    % pxdot = -dH/dx = -(x + 2xy)
    % pydot = -dH/dy = -(y + x^2 - y^2)
    dz = zeros(4,1);
    dz(1) = px;
    dz(2) = py;
    dz(3) = -(x + 2*x*y);
    dz(4) = -(y + x^2 - y^2);
end

%% ======== Event: y=0 crossings with p_y>0, crossing upward ========
function [value, isterminal, direction] = event_y0_pypos(~, z)
    y  = z(2);
    py = z(4);

    % Trigger when y=0. Keep only upward crossings (direction=+1).
    value = y;
    direction = +1;     % y increasing through 0 => upward crossing

    % Stop at the event so we can record the point and restart.
    isterminal = 1;

    % Additionally reject crossings with py<=0 by making event "inactive" there.
    % A simple trick: if py<=0, shift value away from zero so it won't trigger.
    if py <= 0
        value = 1;      % prevents an event root near y=0
    end
end