%=====================================================
%
% OS82 - 2013
% 
% Minimisation d'un syteme masse + pendule rigide
% avec ressort de torsion par Newton-Raphson
%
%=====================================================

%=== nettoyage de matlab
clc;
clear all;
close all;

%=== parametres
L=2;                %longueur barre
m=1;                %masse ponctuelle
g=10;               %acceleration pesanteur
t0=pi/4;            %angle de repos du ressort de torsion
C=10;               %rigidite du ressort

n=100;              %nb points tracés E(x)
epsilon_cvg=1E-6;   %epsilon de convergence pour l'algo
n_iter_max=1000;    %nb d'iterations max toléré

%=== courbe E(theta)
table_t=zeros(n,1);
table_E=zeros(n,1);
for i=1:n
    ti=-pi+2*pi*(i-1)/(n-1);
    table_t(i,1)=ti;
    table_E(i,1)=0.5*C*(ti-t0)^2+m*g*L*sin(ti);
end
figure(1)
hold on
plot(table_t,table_E,'-b');
grid;
xlabel('angle \theta');
ylabel('energie E(\theta)');
hold off

%=== algo Newton-Raphson
ti=0;       %première estimation de l'angle d'équilibre stable
fprintf('* estimation initial -> t0=%f\n',ti);
compteur=1;
test_fin=0;

while(test_fin==0)
    
    fi=C*(ti-t0)+m*g*L*cos(ti);
    fpi=C-m*g*sin(ti);
    tip1=ti-fi/fpi;
    if(abs(tip1-ti)<epsilon_cvg)|(compteur>n_iter_max)
        test_fin=1;
    end
    ti=tip1;
    fprintf('* fin itération %d -> ti=%f\n',compteur,ti);
    
    compteur=compteur+1;
end