%% Loesen der Bewegungsgleichungen verschiedener Magnetsysteme

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

% Vers. 2.1 12/2014

%% 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.9*c;
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)
m = 1836*me;     % Masse des Protons
B0f = m*v0 / (qe*radius)
tmax = Umfang/v0;

x10 = 0.000;
x20 = 0.0;
y0 = [x10;x20;0;0;0;v0]; % 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)
m = 1836*me;     % Masse des Protons
qm = -qe/m;

B0 = 10;%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

%% Dipol
%  function dy = MagnetoptikDipolDgl(t,y,qB0)
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvent(td,yd,qB0m,qB0fm,v0));
tic
[td1,yd1,te,ye,ie] = ode113(@(td,yd) MagnetoptikDipolDgl(td,yd,qB0m,qB0fm,v0),[0,tmax],y0,options);
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
qB0m = 0;
[td2,yd2] = ode113(@(td,yd) MagnetoptikDipolDgl(td,yd,qB0m,qB0fm,v0),[te,tmax],ye,options);
td = [td1;td2(2:end)];
yd = [yd1;yd2(2:end,:)];
toc
%% Quadrupol
tic

options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[tq,yq] = ode113(@(tq,yq) MagnetoptikQuadrupolDgl(tq,yq,qB0m,qB0fm,v0),[0,tmax],y0,options);
toc
%% Quadrupol
tic
[tq,yq] = ode15s(@(tq,yq) MagnetoptikQuadrupolTDgl(tq,yq,qB0m),[0,tmax],y0);
toc