%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         Practice 1          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear Dynamical Systems %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2019.09.18.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear all;
%close all;

% Parameters:
Xinit= -5;%-2;                 % initial value of x        %
Yinit= 3;%5;%2;                  % initial value of y        %
tinit=0;                   % initial time              % These parameters can be changed
Max_Time= 1500;             % end of the simulation     % What happens in longer simulations?
b=0.1;%0.5;%1.98;%2.02;%0.01;%3;% 2;%           % damping                   % Change the damping value! What happens?
 
x =-10:0.05:10;            % not necessary to modify
y =-4:0.05:4;              % not necessary to modify
[X,Y] = meshgrid(x,y);     % makes grid based on x and y vectors (for calculating the energies in each point)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          Exercises          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% How do the systems behave? 
% What are the effecs of the parameters? 
%Try other initial condition and plot the into the same plot!
% Are there stable points in the systems? Zoom into the fix points!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          equation 1         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_1=@(t, xy) [xy(2); -xy(1)]; % 

energy_string=@(X, Y) 0.5*Y.^2+(X.^2)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 2          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_2=@(t, xy) [xy(2); -b*xy(2)-xy(1)]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 3          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_3=@(t, xy) [xy(2); -xy(1)+cos(t)]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 4          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w=0.5;%1.5;%                                         % Parameter (arbitrary)

equ_4=@(t, xy) [xy(2); -xy(1)+cos(w*t)]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 5          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_5=@(t, xy) [xy(2); -2*xy(1)+cos(t)]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 6          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_6=@(t, xy) [xy(2); -b*xy(2)-xy(1)+cos(t)]; %
%equ_6=@(t, xy) [xy(2); -b*xy(2)-xy(1)+cos(0.1*t)]; % for moving fix point

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 7          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_7=@(t, xy) [xy(2); -sin(xy(1))]; %                                                  

energy_pendulum= @(X,Y) 0.5*Y.^2-cos(X)+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 8          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_8=@(t, xy) [xy(2); -b*xy(2)-sin(xy(1))]; %                                                 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 9          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

equ_9=@(t, xy) [xy(2); -sin(xy(1))-b*xy(2)+cos(t)]; %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equation 10         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w=0.9;
b=0.05;
equ_10=@(t, xy) [xy(2); -sin(xy(1))-b*xy(2)+cos(w*t)]; %


Z =energy_pendulum(X,Y); % Energy (string or pendulum)

% 2. initial point
Xinit2= -5;%3;%-2.05;%-0.2007;
Yinit2= 2.99;%2.05;%0.3;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             ODE45              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t, xy]=ode45(equ_9, [tinit Max_Time], [Xinit; Yinit]); %calculate differential equation 
%opt=odeset('RelTol', 1e-3, 'AbsTol', 1e-24); %default RelTol=1e-3; AbsTol=1e-6;
%[t2, xy2]=ode45(equ_9, [tinit Max_Time], [Xinit; Yinit], opt); %calculate differential equation 
[t2, xy2]=ode45(equ_9, [tinit Max_Time], [Xinit2; Yinit2]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         PLOTTING          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        x-y , energy       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure;
subplot(2,2,1);                                            % subplot positions in the figure, calculated from left to rightm from up to down
hold on;
energy_levels=[(linspace(sqrt(0),sqrt(2),5)).^2, (linspace(sqrt(2),sqrt(8),5)).^2];%linspace(0,8,9);                             % energy levels to calculate and plot
%Choose better energy lines to visualise (linear scale, visualise fix points)!
contour(X,Y,Z, energy_levels);                             % plot chosen energy levels
plot(xy(:,1), xy(:,2), 'b');            
plot(xy2(:,1), xy2(:,2), 'r');
xlabel('x');
ylabel('y');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             x-t           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
subplot(2,2,3);
hold on;
plot(xy(:,1), t, 'b');
plot(xy2(:,1), t2, 'r');
xlabel('x');
ylabel('t');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             t-y          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,2);
hold on;
plot(t, xy(:,2), 'b');
plot(t2, xy2(:,2), 'r');
xlabel('t');
ylabel('y');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       x-y-t 3D plot     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,4);
hold on;
plot3(xy(:,1), xy(:,2), t, 'b');
plot3(xy2(:,1), xy2(:,2), t2, 'r');
xlabel('x');
ylabel('y');
zlabel('t');
%}
