%% Visualisierung der Potentiale einer Penningfalle

% konstantes elektrisches Quadrupolfeld und magnetisches Feld

% B0 Magnetfeld in z-Richtung in Tesla
% Ur = U/ro²  U angelegte Spannung r0 Radius der Ringelektrode  
% qm = q/m in As/kg q Ladung; m Masse

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

% Vers. 1.2 11/2014

%% 
U = -0.75; %V
r0 = 3.5e-03;
vorf = U/(2*r0^2);
%% z = 0;
x = linspace(-r0,r0,25);
y = x;
[X,Y] = meshgrid(x,y);
Phiz = vorf*(X.^2 + Y.^2);
subplot(1,2,1)
surf(X*1e03,Y*1e03,Phiz),shg
xlabel('x'),ylabel('y')
%figure
subplot(1,2,2)
contour(X*1e03, Y*1e03, Phiz, 20, 'LineWidth', 2)
hold on        % schuetzt die Abbildung vor loeschen
xlabel('x')
ylabel('y')

% Berechnen und Plotten des elektrischen Feld

[Exz, Eyz] = gradient(Phiz, 0.1);        % "0.1" Ortsaufloesung
quiver(X*1e03, Y*1e03, -Exz, -Eyz, 1, 'k'), shg     % "1" ist ein Skalierungsfaktor, 
axis equal                                      % damit die Pfeile nicht ueberlappen
                                      % 'k': schwarz
                                      

%% y = 0
z = x;
[X,Z] = meshgrid(x,z);
Phiy = vorf*(X.^2 - 2*Z.^2);
figure, subplot(1,2,1)
surf(X*1e03,Z*1e03,Phiy),shg
xlabel('x'), ylabel('z')
%figure
subplot(1,2,2)
contour(X*1e03, Z*1e03, Phiy, 20, 'LineWidth', 2)
axis equal
hold on        % schuetzt die Abbildung vor loeschen
xlabel('x')
ylabel('z')

% Berechnen und Plotten des elektrischen Feld

[Exy, Ezy] = gradient(Phiy, 0.1);        % "0.1" Ortsaufloesung
quiver(X*1e03, Z*1e03, -Exy, -Ezy, 1, 'k'), shg     % "1" ist ein Skalierungsfaktor, 
                                      % damit die Pfeile nicht ueberlappen
                                      % 'k': schwarz
                                      



