%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 2]);

[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 (x(0)=' num2str(x01) ', p(0)=0)']);
ylim([-2 2]*x01);

 

subplot(3,2,[3 4 5 6]);


FFTST=0.02;
FFTSF=1/FFTST;
FFTT=FFTST:FFTST:2000;
NFFT = 2^nextpow2(length(FFTT));
FFTF = FFTSF/2*linspace(0,1,NFFT/2+1);
[T,Y] = ode45(dy,FFTT,[x01, p0],options);
[Th,Yh] = ode45(dyh,FFTT,[x01, p0],options);
[TB,YB] = ode45(dyB,FFTT,[x01, p0],options);
FFTY=fft(Y,NFFT)/length(FFTT);
FFTYh=fft(Yh,NFFT)/length(FFTT);
FFTYB=fft(YB,NFFT)/length(FFTT);

windowSize = 5;

xdata=FFTF(1:windowSize:end);
ydata1=2*abs(FFTYh(1:NFFT/2+1));
ydata2=2*abs(FFTY(1:NFFT/2+1));
ydata3=2*abs(FFTYB(1:NFFT/2+1));

ydata1=filter(ones(1,windowSize)/windowSize,1,ydata1);
ydata1=ydata1(1:windowSize:end);
ydata2=filter(ones(1,windowSize)/windowSize,1,ydata2);
ydata2=ydata2(1:windowSize:end);
ydata3=filter(ones(1,windowSize)/windowSize,1,ydata3);
ydata3=ydata3(1:windowSize:end);


plot(xdata, ydata1, 'b-', xdata, ydata2, 'r-', xdata, ydata3, 'g-')
xlim([0 0.3])
xlabel('Schwingungsfrequenz f');
ylabel('Fourier-Anteil der Frequenz');
title('\bf Fourier-Spektrum der Schwingung')
plot2svg('anharmonic_oscillators_solutions_fft.svg', h)