%% Relativistische Bahn: Sextupolmagnet

% Berechnung der notwendigen Daten zur Visualisierung der
% Abbildungseigenschaften.
% Wird von MagnetoptikquadsextAbb2 verwandt

% 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);
c = 2.998e08;               % m/s
v0= 0.999*c;
qe = 1.602e-19;             % Elementarladung 
me = 9.1096e-31;            % Masse des Elektrons
mp = 1836*me;               % Masse des Protons
m = mp/sqrt(1-v0.^2/c.^2);  % Rel. Masse
B0f = m*v0 / (qe*radius)    % Fuehrungsmagnetfeld -> Kreis
tmax = Umfang/v0;

qm = -qe/m;

B0 = 7.3;                   % Magnetfeld 1. Sextupol
%RDipol = m*v0/(qe*B0);

winkel = linspace(0,2*pi);  % Anfangswerte
r = 0.0005;                 % Anfangswerte 
xa = [];
ya = [];
xe = [];
ye = [];
for alpha = winkel
x10 = r*sin(alpha);
x20 = r*cos(alpha);
xa = [xa;x10];              % Ausgangsprofil
ya = [ya;x20];
y0 = [x10;x20;0;0;0;v0];
qB0m = B0*qm;               % Magnetfeld Sextupol
qB0fm = B0f*qm;             % Magnetfeld Dipol

%% Sextupol
% Magnet 1
%tic
rpos=6;
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvents(td,yd,qB0m,qB0fm,v0,rpos));
[ts1,ys1,tse1,yse1,ie1] = ode113(@(ts1,ys1) MagnetoptikSextopolDgl(ts1,ys1,qB0m,qB0fm,v0),[0,tmax],y0,options);
%
rpos=v0*tmax/2;
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvents(td,yd,qB0m,qB0fm,v0,rpos));
[ts2,ys2,tse2,yse2,ie2] = ode113(@(ts2,ys2) MagnetoptikSextopolDgl(ts2,ys2,0,qB0fm,v0),[tse1,tmax],yse1,options);
rpos=rpos+6;
xcorr=ys2(end,1);  % Korrektur auf lokale Sextupolkoordinaten
zcorr=ys2(end,3);
% Magnet 2
qB0mc=1.3*qm;
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'Events',@(td,yd) MagnetoptikEvents(td,yd,qB0m,qB0fm,v0,rpos));
[ts3,ys3,tse3,yse3,ie] = ode113(@(ts3,ys3) MagnetoptikSextupolDglsys(ts3,ys3,qB0mc,qB0fm,v0,xcorr,zcorr),[tse2,tmax],yse2,options);
%
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[ts4,ys4] = ode113(@(ts4,ys4) MagnetoptikSextopolDgl(ts4,ys4,0,qB0fm),[tse3,tmax],yse3,options);
ts = [ts1;ts2(2:end,:);ts3(2:end);ts4(2:end)];
ys = [ys1;ys2(2:end,:);ys3(2:end,:);ys4(2:end,:)];
%toc
xe =[xe;ys(end,1)];     % Endprofil
ye =[ye;ys(end,2)];
end
%% Visualisierung
% figure,plot3(ys(:,1),ys(:,2),ys(:,3)),shg, grid on
% 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