# 
# Typische Bahnkurve in der Penningfalle und elektr. Potential
#
# Alex Weisse 5/6/2008
# 

# Bewegungsgleichung und num. Loesung
function doty = penning(y,t)
  b = 10;
  doty = [y(4), y(5), y(6), y(1)+b*y(5), y(2)-b*y(4), -2*y(3)];
endfunction

t = 0:0.01:100;
y0 = [1, 0, 0.5, 0, 2, 0];

y = lsode("penning",y0,t);

# Plot der Bahn
vtk_blackonwhite();
vtk_line3(y(:,1),y(:,2),y(:,3));

# Contour plots des el. Potentials 
[a,b] = meshgrid(linspace(-1.5,1.5,30),linspace(-1.5,1.5,30));
c = 0.5*(a.^2 + b.^2);

vtk_surf(a,b,zeros(30),c,"Opacity",0.7);

[a,b] = meshgrid(linspace(-1.5,1.5,30),linspace(-0.6,0.6,30));
c = 0.5*(a.^2 - 2*b.^2);

vtk_surf(a,zeros(30),b,c,"Opacity",0.7);

