%% Zur Visualisierung eines geladenen Teilchens in einer Penningfalle
% Berechnung der Werte
% konstantes elektrisches und magnetisches Feld

% Aufruf:  konstMagnquadEfeldDGL(t,y,qm,B0,Ur)
% 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

tmax = 5e-05;%6;%5;
y0 = [0.00;0;0.00;100;100;100]; % x-y-z vx-vy-vz
qe = 1.602e-19; % C Elementarladung des Elektrons 
me = 9.1096e-31; % Masse des Elektrons
m = 1836*me*12; % Masse C-12 (5-fach ionisiert)
qm = qe/m;%5*qe/m;
B0 = 3.7; %T
U = -0.75; %V
r0 = 3.5e-03 ; % 7 mm 
Ur = U/(r0^2);

tic
[t,y] = ode15s(@(t,y) konstMagnquadEfeldDGL(t,y,qm,B0,Ur),[0,tmax],y0);
toc
% Zusatzdaten
v_xy = sqrt(y0(4)^2 + y0(5)^2);
r_xy = v_xy*tmax*1000 % in mm  freies Teilchen radial
z = y0(5)*tmax*1000 % in mm    freies Teilchen z-Richtung
f = sqrt(-2*Ur*qm)/(2*pi)

%% Plotten
plot3(y(:,1),y(:,2),y(:,3)),shg
%% Ausdehnung Vergleich mit freien Teilchen
figure
plot3(y(:,1),y(:,2),y(:,3)),shg
hold on
pause(3)
xlim([-r_xy,r_xy]),ylim([-r_xy,r_xy]),zlim([0,z])
shg