%% Beispiel fuer die Nullstellenberechnung zu Gl. (1.95)
% MATLAB-Funktion Roots

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% Version 1.0  2016

%% Konstanten
% Traegheitmoment typisch 3/10 m r^2 Beispiel Vollkegel
%
g = 9.81;
m = 0.51;
r = 0.35;
h0 = 0.5;
%
I1 = 3/10*m*r^2;
I3 = 3/15*m*r^2;
% Start mit Rotation um die Figurenachse
dphi0 = 0;
dpsi0 = 100.;
dtheta0 = 0.;
%
theta0 = 60/180*pi;
phi0 = 0;
psi0 = 0;
%%
A = dphi0*(I1*sin(theta0)^2 + I3*cos(theta0)^2)+dpsi0*I3*cos(theta0);
B = I3*(dphi0*cos(theta0) + dpsi0);
Es = 1/2*I1*(dphi0^2*sin(theta0)^2 +dtheta0^2) + m*g*h0*cos(theta0);
%

%
c0 = 1/(2*I1) * A^2 - Es;
c1 = m*g*h0 - A*B/I1;
c2 = B^2/(2*I1) + Es;
c3 = -m*g*h0;

C = [c3,c2,c1,c0];
uNull = roots(C);       % Nullstellen des Potentials
%
uNull = sort(uNull);
ui = linspace(uNull(2),uNull(1),100);
fintegral = @(x) 1./sqrt(-2/I1*polyval(C,x));
%%
t(1) = 0;
for n = 1:length(ui) 
    t(n)=quadgk(fintegral,ui(n),uNull(2));
end
%%
theta =acos(ui);


