%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         Practice 2          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear Dynamical Systems %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2019.09.25.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear all;
%close all;

% damped pendulum equation: ddx+b*dx+sin(x)=0

%Exercises: analyse the system:
% - Show the equilibrium points!
%   - What kinds of fix points they are?
%       - Try to visualise the stable (stable spiral, sink) and the instable (saddle point) fix points!
%       - Which fix points are easy to show?
%           - Determine the separation line (separatrix) between two stable fix points!
%   - What is difference in cases of different damping parameter (see the shapes of the trajectories)?
% - Bolzano shooting: approach a saddle point with 1.0e-6 precision! (Give the initial x,y coordinates of the trajectory!) 

% parameters:
Xinit= -5;                 % initial value of x        %
Yinit= 2;                  % initial value of y        %
tinit=0;                   % initial time              % 
Max_Time= 150;             % end of the simulation     % 
b=0.5;                     % damping                   % Change the damping value! What happens?

% 1. energy level curves:
x =-6*pi:0.05:6*pi;        % resolution of the grid
y =-6:0.05:6;              % resolution of the grid
[X,Y] = meshgrid(x,y);     % makes grid based on x and y vectors (for calculating the energies in each point)
energy_pendulum=; % energy function: E=-cox(x)+1+(y^2)/2
Z =energy_pendulum(X,Y); % Energy level curves of the pendulum
% plot the energy levels with the contour function

% 2. calculate the gradient of the vector field in each grid point:
x2 =-6*pi:pi/5:6*pi;        % resolution of the grid
y2 =-6:0.5:6;              % resolution of the grid
[X2,Y2] = meshgrid(x2,y2);     % makes grid based on x and y vectors (for calculating the energies in each point)
% vector field:
dx=; % change of x
dy=; % change of y
% direction field (vector field after normalisation):
%dx=; %normalisation
%dy=;
%plot the vector field (quiver function)
% Which is the more informative/useful? (vector field or direction field)

% equation:
equ_8=; %damped pendulum equation: ddx+b*dx+sin(x)=0

% visualisation
figure;
hold on;



[t, xy]=ode45(equ_8, [tinit Max_Time], [Xinit; Yinit]); %calculate differential equation
plot(xy(:,1), xy(:,2), 'r');

quiver(X2,Y2,dx, dy, 'k'); % calculate the change of x and y in every X,Y grid ponts, and draw a vector from it starting from the X,Y point
energy_levels=[(linspace(sqrt(0),sqrt(2),5)).^2, (linspace(sqrt(2),sqrt(8),5)).^2, (linspace(sqrt(8),sqrt(16),5)).^2];% energy levels to calculate and plot
contour(X,Y,Z, energy_levels); % plot chosen energy levels
xlabel('x');
xticks(-6*pi:pi:6*pi)
ylabel('y');
xlim([-6*pi, 6*pi]);
ylim([-6, 6]);
hold off;