%% Schwerer symmetrischer Kreisel
% Praezessions- und Nutationsbewegung

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

%% Visualisierung des Potentials V(u) und 
% Bezeichnungen A, B, u vgl. Buch
theta0 = 60/180*pi;
theta1 = 30/180*pi;
u0 = cos(theta0);
u1 = cos(theta1);
%
I1 = 0.05;%0.01;%0.015%0.15%;0.2;  % unterschiedliche Werte erproben
B = 0.6;1.25;%0.6; %1.25;
A = B*u1*0.75;%0.625;
%
V = conv([1 -u0],[1 -u1]);   % V(u) berechnen
V = -conv(V,[1 -2]);
%
ulv = linspace(0,2.01);
plot(ulv,polyval(V,ulv)),grid on
title('Potential V(u)'), shg
%
ul = linspace(u0,u1,200);
%
dphiN = (1 - B/A.*ul)./(1 - ul.^2);
figure, plot(ul,dphiN), title('dphiN'),grid on, shg
%
phiInt = @(x) A/I1 * (1-B/A*x)./(1-x.^2) ./ sqrt(-2/I1 .* polyval(V,x));
figure,plot(ul,phiInt(ul)), title('u,dphi'),grid on, shg
%
ul0=ul(1);
ul(1)=ul(1)*1.001;       %Nullstelle von V
ul(end)=ul(end)/1.001;
for n = 1:length(ul) 
    phi(n)=quadgk(phiInt,ul(n),ul0)/pi*180;      % Integration
    theta(n) = acos(ul(n))/pi*180;
end
figure,plot(phi,theta),title('phi - theta'), grid on, shg
%% Kugel fuer Visualisierung plotten

x = sin(phi/180*pi).*sin(theta/180*pi);
y = cos(phi/180*pi).*sin(theta/180*pi);
z = cos(theta/180*pi);
figure, sphere, hold on
plot3(x,y,z),shg
axis equal
%% Symmetrische Loesung dazupacken und Visualisieren
% thetan = [theta,fliplr(theta)];
% deltaphi=fliplr(diff(phi));
% deltaphi = [phi(end),deltaphi];
% deltaphitoo = cumsum(deltaphi);
% phin=[phi,deltaphitoo];
% x = sin(phin/180*pi).*sin(thetan/180*pi);
% y = cos(phin/180*pi).*sin(thetan/180*pi);
% z = cos(thetan/180*pi);
% figure, sphere, hold on
% plot3(-x,-y,z),zlim([0,1]),shg
%%  Abb. 1.33 Ueberlagerung Kugel - Nutationsbewegung 
thetan2 = [theta,theta]+90; % nur ein Ast der Loesung wurde berechnet s. Abb. phi - theta
phin2=[phi,-phi];
x = sin(phin2/180*pi).*sin(thetan2/180*pi);
y = cos(phin2/180*pi).*sin(thetan2/180*pi);
z = -cos(thetan2/180*pi);
figure, sphere, hold on
plot3(-x,-y,z),zlim([0,1]),shg