%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% three function files 
function [ dy ] = harmonic( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1))/m;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ dy ] = anharmonicA( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^2)/m;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ dy ] = anharmonicB( t,y,m,k,l,c )
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=-(c*y(2)+k*y(1)+l*y(1)^3)/m;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
addpath('plot2svg');

m=1;
k=1;
l=0.2;
c=0;
x0=0.2;
x01=2.4;
x02=0;
p02=2;
p0=0;

PSX0range=0:(0.4):(2.4);

PSRange=[-4.5 4.5];
 
 
 
dy=@(t,y) anharmonicA(t,y,m,k,l,c);
dyB=@(t,y) anharmonicB(t,y,m,k,l,c);
dyh=@(t,y) harmonic(t,y,m,k,l,c);
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
 
 
h=figure(1);
subplot(3,2,1);
[T,Y] = ode45(dy,[0 40],[x0, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x0, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x0, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
ylim([-1.5 1.5]*x0);
xlabel('Zeit t');
ylabel('Ort x(t)');
title(['\bf Trajektorie für kleine Anfangsauslenkung (x_0=' num2str(x0) ', p_0=0)']);
 
subplot(3,2,2);
for x0=PSX0range
    [Th,Yh] = ode45(dyh,[0 10],[x0, p0],options);
    plot(Yh(:,1),Yh(:,2), 'b-');
    if x0==PSX0range(1)
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
daspect([1 1 1]);
xlabel('position x');
ylabel('impulse p');
title(['\bf Phasenraum-Plot, harmonischer Oszillator']);
 
 
subplot(3,2,3);
[T,Y] = ode45(dy,[0 40],[x01, p0],options);
[Th,Yh] = ode45(dyh,[0 40],[x01, p0],options);
[TB,YB] = ode45(dyB,[0 40],[x01, p0],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-');
xlabel('Zeit t');
ylabel('Ort x(t)');
title(['\bf Trajektorie für große Anfangsauslenkung (x_0=' num2str(x01) ', p_0=0)']);
ylim([-2 2]*x01);
subplot(3,2,4);
for x0=PSX0range
    [T,Y] = ode45(dy,[0 10],[x0, p0],options);
    plot(Y(:,1),Y(:,2), 'r-');
    if x0==PSX0range(1)
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
daspect([1 1 1]);
xlabel('Ort x');
ylabel('Impuls p');
title(['\bf Phasenraum-Plot, asymmetrisches Kraftgesetz']);
 
 
subplot(3,2,5);
[T,Y] = ode45(dy,[0 40],[x02, p02],options);
[Th,Yh] = ode45(dyh,[0 40],[x02, p02],options);
[TB,YB] = ode45(dyB,[0 40],[x02, p02],options);
plot(Th,Yh(:,1), 'b-', T,Y(:,1),'r-', TB,YB(:,1),'g-', [0 40], [2 2], 'b:', [0 40], -[2 2], 'b:');
xlabel('Zeit t');
ylabel('Ort x(t)');
title(['\bf Trajektorie bei gleicher Energie (x_0=' num2str(x02) ', p_0=' num2str(p02) ')' ]);
ylim([-2.5 2.5]);
subplot(3,2,6);
for x0=PSX0range
    [TB,YB] = ode45(dyB,[0 10],[x0, p0],options);
    plot(YB(:,1),YB(:,2), 'g-');
    if x0==PSX0range(1)
        hold on
    end
end
hold off
xlim(PSRange);
ylim(PSRange);
daspect([1 1 1]);
xlabel('Ort x');
ylabel('Impuls p');
title(['\bf Phasenraum-Plot, symmetrisches Kraftgesetz']);

% uses the plot2svg tool from MatlabCentral by Juerg Schwizer 
plot2svg('anharmonic_oscillators_solutions.svg', h)