%% Relativistische Bahn: Sextupolmagnet

% 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);
%tmax = 3e-04;%-07;%-05;%-06;
c = 2.998e08; % m/s
v0= 0.999*c;
qe = 1.602e-19; % C Elementarladung des Elektrons 
me = 9.1096e-31; % Masse des Elektrons
%mC = 1836*me*12; % Masse C-12 (5-fach ionisiert)
mp = 1836*me;     % Masse des Protons
m = mp/sqrt(1-v0.^2/c.^2);
B0f = m*v0 / (qe*radius)
tmax = Umfang/v0;

x10 = 0.00;
x20 = 0.001;
y0 = [x10;x20;0;0;0;v0]; % x-y-z vx-vy-vz
qe = 1.602e-19; % C Elementarladung des Elektrons 
qm = -qe/m;

B0 = 10;0;%0.4;%0.1;%1;%0; %T 
% führt für Dipolfeld zu einer Kreisbewegung mit Radius
RDipol = m*v0/(qe*B0)
qB0m = B0*qm
qB0fm = B0f*qm

%% Sextopol
tic
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvent(td,yd,qB0m,qB0fm,v0));
[ts1,ys1,tse,yse,ie] = ode113(@(ts1,ys1) MagnetoptikSextopolDgl(ts1,ys1,qB0m,qB0fm,v0),[0,tmax],y0,options);
%
qB0m=0;
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[ts2,ys2] = ode113(@(ts2,ys2) MagnetoptikSextopolDgl(ts2,ys2,qB0m,qB0fm),[tse,tmax],yse,options);
ts = [ts1;ts2(2:end,:)];
ys = [ys1;ys2(2:end,:)];
toc

%% Visualisierung
figure,plot3(ys(:,1),ys(:,2),ys(:,3)),shg
figure,plot(ys(:,1),ys(:,2)),xlabel('x'),ylabel('y'),grid on
figure,plot(ys(:,1),ys(:,3)),xlabel('x'),ylabel('z'), grid on
figure,plot(ys(:,2),ys(:,3)),xlabel('y'),ylabel('z'), grid on
%% Fehler
hoehe = 0.5*9.81*tmax^2