close;

xmax=1.0;

x = linspace(-xmax,xmax);
y = linspace(-xmax,xmax);
[X,Y] = meshgrid(x,y);

V = 0.5*(X.^2+Y.^2) + X.*X.*Y-(1/3)*Y.*Y.*Y;

figure(1)
[C, h]=contour(X,Y,V,16);
axis equal

xlabel 'Position x'
ylabel 'Position y'
clabel(C,h)
h.LineWidth = 2;

%

% sympletic Euler integrator 

clear
close 

Tmax = 10000;
dt   = 0.01;
Tsteps = round(Tmax/dt) +1;

t=zeros(Tsteps,1);

x=zeros(Tsteps,1);
y=zeros(Tsteps,1);
px=zeros(Tsteps,1);
py=zeros(Tsteps,1);

x0 = 0;

%e0  = 0.5*1/8;
e0  = 1/8;
y0  = 0.1;
py0 = 0.14;
px0 = sqrt(2*e0-py0^2-y0^2+(2/3)*y0^3 );

t(1)=0;

x(1)=x0;
y(1)=y0;

px(1)=px0;
py(1)=py0;

s=0;

pc1=zeros(Tmax,1) + y0 ;
pc2=zeros(Tmax,1) + py0 ;

for n=1:Tsteps-1

    t(n+1)   = t(n) + dt;

    px(n+1) = px(n) + dt*( -x(n) - 2*x(n)*y(n) );
    py(n+1) = py(n) + dt*( -y(n) + y(n)^2-x(n)^2);
 
    x(n+1)  = x(n) + dt* px(n+1);
    y(n+1)  = y(n) + dt* py(n+1);

    if x(n+1) * x(n) < 0

        s=s+1;
        pc1(s) = y(n);
        pc2(s) = py(n);

    end
end

figure(2)
T10= round(Tsteps/20);
scatter(t(1:T10), x(1: T10 ) , '.')
xlabel 'Time t' 
ylabel 'Position x' 


ham=0.5*(px.^2+py.^2) +0.5*(x.^2+y.^2) + x.*x.*y-(1/3)*y.*y.*y;

figure(3)
scatter(t(1:100:end), ham(1:100:end) , 'ro');
hold on;
plot([0,Tmax], [e0,e0] , 'k-', LineWidth=2);
xlim([0,Tmax])
ylim([0,2*e0])

xlabel 'Time t' 
ylabel 'Energy H' 

hold off;

figure(4)
scatter(pc1, pc2,'.')
hold on;
scatter(pc1(Tmax), pc2(Tmax),'o')
title(['Poincare section H=', num2str(e0)])
xlabel 'Position y'
ylabel 'Momentum p_y'
axis equal
xlim([-0.7,0.7])
ylim([-0.7,0.7])
hold on;
if e0==1/8
yb=[-0.439694:0.00001:0.673648];
pbp=(sqrt(6*e0 - 3*yb.^2 + 2*yb.^3)/sqrt(3));
pbm=-(sqrt(6*e0 - 3*yb.^2 + 2*yb.^3)/sqrt(3));
plot(yb,pbm, 'k', LineWidth=2);
plot(yb,pbp, 'k', LineWidth=2);
end



%

% NON-sympletic Euler integrator does not conserve energy

close 

Tmax = 67;
dt   = 0.01;
Tsteps = Tmax/dt;

t=zeros(Tsteps,1);

x=zeros(Tsteps,1);
y=zeros(Tsteps,1);
px=zeros(Tsteps,1);
py=zeros(Tsteps,1);

x0 = 0;

e0  = 1/8;
y0  = 0.1;
py0 = 0.14;
px0 = sqrt(2*e0-py0^2-y0^2+(2/3)*y0^3 );

t(1)=0;

x(1)=x0;
y(1)=y0;

px(1)=px0;
py(1)=py0;

for n=1:Tsteps-1

    t(n+1)  = t(n) + dt;

    px(n+1) = px(n) + dt *( -x(n) - 2*x(n)*y(n) );
    py(n+1) = py(n) + dt *( -y(n) + y(n)^2-x(n)^2);
 
    x(n+1)  = x(n) + dt* px(n); % px(n) instead of px(n+1)
    y(n+1)  = y(n) + dt* py(n);

end

figure(2)
scatter(t,x, 'o')
xlabel 'Time t' 
ylabel 'Position x' 




ham=0.5*(px.^2+py.^2) +0.5*(x.^2+y.^2) + x.*x.*y-(1/3)*y.*y.*y;

figure(3)
scatter(t, ham , 'ro');
hold on;
plot([0,Tmax], [e0,e0] , 'k-', LineWidth=2);
ylim([0,2*e0])

xlabel 'Time t' 
ylabel 'Energy H' 

hold off;

figure(4)
scatter(y,py,'.')



