%% Relativistische Bahn: Quadrupolmagnet

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

% Vers. 2.1 01/2015

%% Parameter
% Basieren auf den Daten des HCL
Umfang = 26658.883;         %m
radius = Umfang/(2*pi);
c = 2.998e08;               % m/s  Lichtgeschwindigkeit
v0= 0.999*c;                % Teilchengeschwindigkeit
qe = 1.602e-19;             % Elementarladung 
me = 9.1096e-31;            % Masse des Elektrons
mp = 1836*me;               % Ruhemasse des Protons 
m = mp/sqrt(1-v0.^2/c.^2);  % relativistische Masse
B0f = m*v0 / (qe*radius);   % Führungsfeld Dipol -> Kreisbahn
tmax = Umfang/v0;           % Umlaufdauer

x10 = 0.001;
x20 = 0.001;
y0 = [x10;x20;0;0;0;v0];    % x-y-z vx-vy-vz  Anfangswerte
qm = -qe/m;

B0 = 10;                    % T Quadrupolfeld
qB0m = B0*qm;               % -> an Differetialgleichung             
qB0fm = B0f*qm;

%% Quadrupol
tic
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvent(td,yd,qB0m,qB0fm,v0));
[tq1,yq1,tqe,yqe,ie] = ode113(@(tq1,yq1) MagnetoptikQuadrupolDgl(tq1,yq1,qB0m,qB0fm,v0),[0,tmax],y0,options);
%
qB0m=0;                     % Bahn nach Quadrupol
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[tq2,yq2] = ode113(@(tq2,yq2) MagnetoptikQuadrupolDgl(tq2,yq2,qB0m,qB0fm),[tqe,tmax],yqe,options);
tq = [tq1;tq2(2:end,:)];
yq = [yq1;yq2(2:end,:)];
toc
%% Visualisierung
figure,plot3(yq(:,1),yq(:,2),yq(:,3)),shg
