% Struktur des Programms Pendel.m
Tend=10;h=0.01;N=floor(Tend/h);
% Anfangswerten
q=zeros(3,1);p=zeros(3,1);
q(1)=sin(1.3);q(2)=0;q(3)=cos(1.3);
p(1)=3*cos(1.3);p(2)=6.5;p(3)=-3*sin(1.3);
ham0= ... ; % Anfang gesamte Energie 
...
spher0=norm(q); % um zu sehen ob man auf der Sphäre bleibt
...
for j=1:N
  % Bestimmung von lambda_{n+1} so, dass
  % g(q_{n+1})=0 (bei Hand, explizit!)
  ...
  % Berechnung von phut_{n+1} und q_{n+1}
  ...
  % Bestimmung von mu_{n+1} so, dass
  % G(q_{n+1})H_p(p_{n+1},q_{n+1})=0 (bei Hand, explizit!)
  ...
  % Berechnung von p_n+1
  ...
  % Berechnung Energie+Sphäre
  % Einspeichern für die Plots
  ...
end  
% Plot Energie+Sphäre
...
% Schöne Visualisierung
fig=figure;set(fig,'doublebuffer','on')
xlim([-1 1]);ylim([-1 1]);zlim([-1 1])

p=plot3(q1plot(1),q2plot(1),q3plot(1),'ro');
for j=2:N
  set(p,'xdata',q1plot(1:j),'ydata',q2plot(1:j),'zdata',q3plot(1:j));
  drawnow;
end  
