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

%clear all;
%close all;

% Parameters:
Xinit= -2;                 % initial value of x        %
Yinit= 2;                  % initial value of y        %
tinit=0;                   % initial time              % These parameters can be changed
Max_Time= 150;             % end of the simulation     % What happens in longer simulations?
b=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)

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          Reminder          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%0.1. basic plotting
%{
figure                                  % new figure window
plot(x,sin(x));                         % plot (many other parameters)   'r' = 'red', 'r*'
waitforbuttonpress();
plot3(x,sin(x),cos(x))                  % 3D plot, the 3 parameters need to have the same size
waitforbuttonpress();
figure
subplot(2,2,1:2)
text(0.5,0.5,'Because we like plots! subplot(2,2,1:2)',...
    'FontSize',14,'HorizontalAlignment','center')
subplot(2,2,3)
text(.5,.5,'see the help menu - subplot(2,2,3)',...
    'FontSize',12,'HorizontalAlignment','center')
subplot(2,2,4)
text(.5,.5,'for the details - subplot(2,2,4)',...
    'FontSize',12,'HorizontalAlignment','center')
% hold on/off, xlabel('something'), 
%}


% 0.2. understanding meshgrid:
%{
% [X,Y] = meshgrid(1:3,10:14)
% X =
%      1     2     3
%      1     2     3
%      1     2     3
%      1     2     3
%      1     2     3
% Y =
%     10    10    10
%     11    11    11
%     12    12    12
%     13    13    13
%     14    14    14
figure 
[X,Y] = meshgrid(-2:.2:2, -2:.2:2);                                
Z = X .* exp(-X.^2 - Y.^2); %x*e^(-x^2-y^2)                                        
surf(X,Y,Z)
%}

% 0.3. contour megértése:
%{
load penny;                 % load image (P(128*128))
figure;
x = linspace(1,15,128);     % generates 128 values between 1 and 15 evenly distributed
y = linspace(1,15,128);
contour(x,y,flipud(P),5);  % flipud - mirror matrix elements to the central axis 
                            % (different indexing for images and plots)
axis square;                % equal axis division
colormap('Bone');           % other colormaps: Jet, pink, winter, ....
%}
%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          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=; 

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

%w=;                                         % Parameter (arbitrary)

%equ_4=; 

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

%equ_5=; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         equations 6         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%equ_6=;                                                    

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

%equ_7=;   

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

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

%equ_8=;                                                    

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

%equ_9=; 

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

%w=0.9;
%b=0.05;
%equ_10=;


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

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             ODE45              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[t, xy]=ode45(equ_2, [tinit Max_Time], [Xinit; Yinit]); %calculate differential equation 

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         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(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');            
xlabel('x');
ylabel('y');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             x-t           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
subplot(2,2,3);
plot(xy(:,1), t, 'b');
xlabel('x');
ylabel('t');

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

subplot(2,2,2);
plot(t, xy(:,2), 'b');
xlabel('t');
ylabel('y');

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

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