clear all;

global om0 gamma

gamma=0.1;
om0=1;

% initial y(0) and v(0)
y0=1;
v0=2;

% to check
% [ydot,vdot ]=twoD(1,2)

%how many periods?
nper=5;

T=linspace(0,nper*2*pi,nper*40);
d= T(2)-T(1);

Y=zeros(1,length(T));
V=zeros(1,length(T));


Y(1)=y0;
V(1)=v0;


for k=2:length(T);
    
    y=Y(k-1);
    v=V(k-1);
    
    [ky1,kv1] = twoD( y         , v         );
    [ky2,kv2] = twoD( y+d*ky1/2 , v+d*kv1/2 );
    [ky3,kv3] = twoD( y+d*ky2/2 , v+d*kv2/2 );
    [ky4,kv4] = twoD( y+d*ky3   , v+d*kv3   );
    
    Y(k) = y + d/6*(ky1+2*ky2+2*ky3+ky4);
    V(k) = v + d/6*(kv1+2*kv2+2*kv3+kv4);
end


% i theoritiki lysi
omg = sqrt(om0^2-gamma^2);
A = -1i*(v0+y0*(gamma+1i*omg))/(2*omg);
B =  1i*(v0+y0*(gamma-1i*omg))/(2*omg);
Yth = A*(exp(-gamma*T + 1i*omg*T))+B*(exp(-gamma*T - 1i*omg*T));


% i troxia y(t)
figure(14);clf;
plot(T,Y,'b',T,Yth,'--r','LineWidth',2)
hold on;
plot(T,0*T,'--k');
hold off;
xlabel('t','FontSize',18);
ylabel('y(t)','FontSize',18);
title(['damped oscillator with \omega_0=' num2str(om0) ' and \gamma=' num2str(gamma)],'FontSize',18);
xlim([T(1) T(end)]);
h=legend('numerical','theoretical');

maxd = sqrt(y0^2 + om0^2*v0^2);

figure(15);clf;
plot(Y,V,'k')
hold on;
plot(Y(1),V(1),'or','MarkerSize',8)
hold off;
axis square;
axis([-1.2 1.2 -1.2 1.2]*maxd)
xlabel('y','FontSize',18);
ylabel('v','FontSize',18);
title('phase space','FontSize',18);