%% Pendulum_and_dirfield.m [H4T8heiluri] 
%  [Tehdään mieluummin suoraan mukaillen ohjetta:]
%%
% <http://math.aalto.fi/opetus/MatOhjelmistot/2014kevat/L/ODEohje.html>
%
%% ODEfunktio "heiluri" voidaan määritellä myös funktiokahvana:
heiluri=@(t,y) [y(2);-sin(y(1))]
Tvali=[0 20]; y0=[1;1];        % Alkukulma=1, alkukulmanopeus=1. 
[T,Y]=ode45(heiluri,Tvali,y0); % ode45 havaitsee funktiokahvan, sille
                               % ei tarvitse sitä enää @-merkillä
                               % ilmoittaa. (Muistisääntö 1 @-merkki
                               % joko määrittelyssä tai kutsussa, ei molemmissa.)
size(T)  % Aikapistevektorin pisteiden lkm.
size(Y)  % Kaksisarakkeinen matriisi, ratkaisut y1 ja y2 aikapisteissä t.
plot(T,Y(:,1),T,Y(:,2))  % Aikakuva: heilahduskulma ja kulmanopeus ...
                         %              ajan funktiona
figure
plot(Y(:,1),Y(:,2))     % Faasikuva: Kulmanopeus kulman "funktiona" ...
                        % eli käyrä (Theta(t),omega(t)), 0< t < 20).
%% 3 eri trajektorityypppiä:
close all
ya0=[1;1]; % Pieni alkukulma/nopeus -> heilahtelee edestakaisin 
yb0plus=[-pi;0.01]; % Kriittinen alkutilanne (Heiluri pysytyasennossa+)
yb0miinus=[-pi;-0.01]; % Kriittinen alkutilanne (Heiluri pysytyasennossa-)
yc0=[-5;2];   % Suuret alkuarvot, heiluri pyörii ympäri
yd0=[5;-2];   %    ... toisinpäin
Tvali=[0 50]; 
[ta,Ya]=ode45(heiluri,Tvali,ya0);
[tbplus,Ybplus]=ode45(heiluri,Tvali,yb0plus);
[tbmiinus,Ybmiinus]=ode45(heiluri,Tvali,yb0miinus);
[tc,Yc]=ode45(heiluri,Tvali,yc0);
[td,Yd]=ode45(heiluri,Tvali,yd0);
hold on
plot(Ya(:,1),Ya(:,2),Ybplus(:,1),Ybplus(:,2),Ybmiinus(:,1),Ybmiinus(:,2),...
    Yc(:,1),Yc(:,2),Yd(:,1),Yd(:,2))
axis equal; axis([-6 6 -4 4])
xlabel y_1(t); ylabel y_2(t)
grid on
%% Hiiri avuksi, samalla suuntakentta:
[y1,y2]=meshgrid(-6:0.5:6,-6:0.5:6); % Koordinaattihilapisteistö
dy1=y2; dy2=-sin(y1);  % Heiluridiffyhtälöryhmän derivaatat hilapisteissä
quiver(y1,y2,dy1,dy2,2)  % Suuntanuolet hilapisteisiin derivaattojen suuntaan.
% Toden totta, näillä 3:lla koodirivillä saat suuntakentän!
%% Now, choose initial values pointing on graphics window.
% help ginput
yy0 = [];  % Gather chosen initial points
n = 0;
% Loop, picking up the points.
disp('Move cursor into graphics window')
disp('******************************')
disp(' ')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [y1i,y2i,but] = ginput(1);
[Tr,Yr]=ode45(heiluri,[0 4*pi],[y1i;y2i]); % t= ti...3 (up to right window border)
    plot(Yr(:,1),Yr(:,2))
[Tl,Yl]=ode45(heiluri,[0 -4*pi],[y1i;y2i]); 
    plot(Yl(:,1),Yl(:,2))
    n = n+1;
    yy(:,n) = [y1i;y2i];
end
%%
hold on
plot(yy(1,:),yy(2,:),'r*')  % Add chosen initial points to fig.
