%% Relativistische Bahn: Quadrupolmagnet
% Teilchen starten laengs einer Kreisbahn
% Abbildung der Kreisbahn nach Durchlaufen des Quadrupolmagnets --> Abb.:
% 2.7
% Ausgangspunkt ist ein kreisfoermiges Profil mit einem Radius von 100 mu
% Darstellung dieses Profils nach Durchlaufen eines 6m langen
% Quadrupolmagneten

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

% Vers. 2.1 01/2015

%clear,clc
%% 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;


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
tq = [];
yq = [];
yq0 = [];
yort0 = []
rs = 0.0001;
phil = linspace(0,2*pi,100);
for phi=phil
ort0 = [rs * sin(phi), rs * cos(phi)];
yort0 =[yort0;ort0];
y0 = [ort0';0;0;0;v0]; % x-y-z vx-vy-vz
%% Quadrupol
%tic
options = odeset('RelTol',1e-12,'AbsTol',1e-12,'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;
% options = odeset('RelTol',1e-10,'AbsTol',1e-10);
% [tq2,yq2] = ode113(@(tq2,yq2) MagnetoptikQuadrupolDgl(tq2,yq2,qB0m,qB0fm),[tqe,tmax],yqe,options);
yq = [yq;yqe];
%toc
end
%% Quadrupol
% tic
% [tq,yq] = ode15s(@(tq,yq) MagnetoptikQuadrupolTDgl(tq,yq,qB0m),[0,tmax],y0);
% toc

%% Visualisierung
%figure,plot3(yq(:,1),yq(:,2),yq(:,3)),shg
korr=(min(yq(:,1))+max(yq(:,1)))/2
figure,plot(yort0(:,1),yort0(:,2)),shg
hold on
plot(yq(:,1)-korr,yq(:,2)),shg
legend('Start','nach 6m')
%axis equal