%% Forced, Damped, Haronic Motion
% In this worksheet we study the equation \ddot{x} +b\;\dot{x} +\omega_0^2 x=A\;\cos \;\omega \;t.
%
% Simple harmonic motion case


clear
omega = 1;
shm = @(t,y) [y(2); -omega^2*y(1)];
tspan = [0 30];
y0 = [1 0]';
[t,y] = ode45(shm,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Simple Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('SHM Trajectory')
axis equal

%
% Damped Oscillator

clear
omega = 1;
b=0.1;
dhm = @(t,y) [y(2); -b*y(2)-omega^2*y(1)];
tspan = [0 30];
y0 = [1 0]';
[t,y] = ode45(dhm,tspan,y0);
plot(t,y)
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Damped Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('DHM Trajectory')
axis equal

%
% Forced, damped oscillations as an autonous system. Note the overlapping trajectory in the phase spae.

clear
omega0 = 1;
omega=2;
b=0.1;
A=1;
fdhm = @(t,y) [y(2); -b*y(2)-omega0^2*y(1)+A*cos(y(3)); omega];
tspan = [0 30]; % tspan = 0:.1:30;  % Makes solution smoother
y0 = [1 0 0]';
[t,y] = ode45(fdhm,tspan,y0);
plot(t,y(:,1),t,y(:,2))
legend({'y','dy/dt'})
xlabel('t')
ylabel('y')
title('Forced, Damped Harmonic Motion')

figure
plot(y(:,1),y(:,2))
xlabel('x')
ylabel('y')
title('FDHM Trajectory')
axis equal

%
% Create a Poincaré Surface of Section plot.

period=2*pi/omega;
dt=period/10;
tspan = 0:dt:2000;  
y0 = [1,0,0]';
[t,y] = ode45(fdhm,tspan,y0);


P_x = y(mod(t, period) < 1e-6, 1);
P_y = y(mod(t, period) < 1e-6, 2);

figure
scatter(P_x,P_y,'.')
xlabel('Position');
ylabel('Velocity');
title('Poincaré Surface of Section');

%
% Here we add more initial conditions. They seem to converge to the red point.

% Add more ICs
figure
scatter(P_x,P_y,'.')
xlabel('Position');
ylabel('Velocity');
title('Poincaré Surface of Section');
hold on
for j=1:50
    y0 = [j/50,0,0]';
    [t,y] = ode45(fdhm,tspan,y0);
    P_x = y(mod(t, period) < 1e-6, 1);
    P_y = y(mod(t, period) < 1e-6, 2);
    L=length(P_y);
    Npts=180;
   % scatter(P_x,P_y,'.r')
    scatter(P_x(L-Npts:L),P_y(L-Npts:L),'.r')    
end
hold off