close all
clear all
% Derive DE for spring
syms x dx a da m l g
T   = 1/2 * m * (dx.^2 + (da*(x)).^2);%m*l^2*dx.^2/2;
V   =  m*g*(x)*sin(a) + 1/2 * 400 * (x-l)^2;%m*g*l*(1 - cos(x));


%T   = 1/2 * m * (dx.^2 + (da*(x+l)).^2);%m*l^2*dx.^2/2;
%V   =  -m*g*(x+l)*sin(a) + 1/2 * 30 * x^2;%m*g*l*(1 - cos(x));

L   = T - V;
E   = T + V;
X   = {x dx a da};
Q_i = {0 0}; Q_e = {0 0};
R   = 0;
par = {m g l};
% Solve Lagrange equations and save DE as .m file
VF  = EulerLagrange(L,X,Q_i,Q_e,R,par,'m','Pendulum_sys');

% Solve DE numerically using ode45
m_num   = 1;
g_num   = 9.81;
l_num   = 10;
tspan   = linspace(0, 2*pi, 200);
Y0      = [l_num-0.5 5 90*pi/180 1*pi/180];
options = odeset('RelTol',1e-8);
% Use created .m file to solve DE
[t, Y]  = ode45(@Pendulum_sys,tspan,Y0,options,m_num,g_num,l_num);

% Plot state vector and total energy
%plot(t,Y)
polar(Y(:,3), Y(:,1), 'r')
hold on
theta = linspace(0,2*pi);
rho = l_num*ones(size(theta));
polar(theta,rho, 'b')
title('Ressort')
xlabel('t')
ylabel('{x}(t), d{x}/dt(t)')
legend('Position', 'Vitesse')
% End of script