%% Vergleich relativistische Bewegung - nicht-relativistisch in einem Dipolfeld

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

% Vers. 1.2 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.999*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

qm = -qe/m;

B0 = 0;% muss hier immer Null sein 
% führt für Dipolfeld zu einer Kreisbewegung mit Radius
RDipol = m*v0/(qe*B0)


%% Dipol nichtrelativistisch
B0f = m*v0 / (qe*radius)
tmax = Umfang/v0
qB0m = B0*qm
qB0fm = B0f*qm
%  function dy = MagnetoptikDipolDgl(t,y,qB0)
tic
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
qB0m = 0;
[tdn,ydn] = ode113(@(tdn,ydn) MagnetoptikDipolDgl(tdn,ydn,qB0m,qB0fm,v0),[0,tmax],y0,options);
toc
%% Dipol relativistisch
%tmax = tmax/10
% m,mr = m/sqrt(1-(v0/c).^2)
% B0f = mr*v0 / (qe*radius)
% tmax = Umfang/v0
% qB0m = B0*qm
% qB0fm = B0f*qm
radiusrel = m./sqrt(1-v0^2/c^2) *v0 /(qe*B0f),radius
tmaxr = 2*pi*radiusrel/v0, tmax
tic
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
qB0m = 0;
[tdr,ydr] = ode113(@(tdr,ydr) MagnetoptikDipolRelDgl(tdr,ydr,qB0m,qB0fm,v0,c),[0,tmaxr],y0,options);
toc
%% Visualiserung
figure, plot(ydn(:,1),ydn(:,3)), grid on, shg
hold on
plot(ydr(:,1),ydr(:,3),'m'), axis equal
%% Magnetfeld relativistisch - nichtrelativistisch
v0r = linspace(0,0.999*c);
mr = m./sqrt(1-(v0r/c).^2);
B0fr = mr.*v0r / (qe*radius);
B0fn = m.*v0r / (qe*radius);
figure,plot(v0r,B0fr,'k',v0r,B0fn,'k')
xlabel('v_0'),ylabel('B'),shg
hold on
gca = axes('Position',[0.2400 0.2200 0.60 0.60]);shg
plot(gca,ydn(:,1),ydn(:,3),'k'), grid on, shg
hold on
plot(gca,ydr(:,1),ydr(:,3),'k')
xlabel('x'),ylabel('z'),axis equal
