Per la cronaca: ho scritto il codice per risolvere l'equazione:
f = @(z,y) [y(4:6); (q/m)* cross( cross( E, 1/y(4:6) ), 1/y(4:6) )+ (q/m)*cross( B, 1/y(4:6) ) ];
[z,y] = ode23t(f,zspan,y0);
fa cinque minuti di conti e poi fa una minchiata, quindi penso sia sbagliata l'equazione...
Poi ho provato così:
eV=1.6021765e-19; %[J], fattore di conversione eV-->J
c= 2.99792458e8; %[m/s], velocità della luce
%KB=(8.617332*10^(-5))*eV; %[eV/K], costante di Boltzmann
uma1=931.494028; %[MeV/c^2], fattore di conversione uma-->MeV/c^2 (massa)
KB=8.617e-5; %[eV*K^-1], costante di Boltzmann
KBT=1076*KB;
E=0.03; %Eneriga 30KeV
K_ion=E;
Etot=K_ion+((m_H*uma1)/uma); % ho tutto in MeV
betasquare_H=1-(((m_H*uma1)/(uma*Etot))^2);
v_H=sqrt(betasquare_H*c^2);
v0 = [0 0 v_H]'; %initial velocity column vector
r0 = [0 0 0]'; %initial position of particle column vector
y0 = [r0; v0]; %Concatena i due vettori colonna r0 e v0 e crea il vettore colonna delle condizioni iniziali
% m = 2; % mass
% q = 1; %charge on particle
q_over_m=q/m_H;
dt= 0.01; %intervallo di tempo
tspan = [0 dt];
%%
%Distanze in metri
Deriva1= 0.05;
LB=0.15;
inizioE=0.13;
fineB= 0.20;
LE=0.07;
fineE=0.21;
Dist_Spettr_riv=0.40;
z=0; %posizione iniziale
traiettoria= zeros(11,6);
while z<Deriva1
B=[0 0 0]';
E=[0 0 0]';
f = @(t,y) [y(4:6); (q_over_m).*cross(y(4:6),B)+(q_over_m).*E]; % the expression [y(4:6); (q/m)*cross(y(4:6),B)]
% combines two vectors of length 3 to make a vector of
% length 6(as long as y is a column vector).
[t,y] = ode23t(f,tspan,y0);
traiettoria= [traiettoria;y];
tspan = dt+tspan;
y0=y(end,:);
z=y(end,3);
end
figure
plot3(traiettoria(:,3),traiettoria(:,2),traiettoria(:,1))
grid on
xlabel ('direzione del fascio');
zlabel ('deflessione elettrica');
ylabel ('deflessione magnetica');
%%
while (z>=Deriva1 || z<inizioE)
B=[1 0 0]';
E=[0 0 0]';
f = @(t,y) [y(4:6); (q_over_m).*cross(y(4:6),B)+(q_over_m).*E]; % the expression [y(4:6); (q/m)*cross(y(4:6),B)]
% combines two vectors of length 3 to make a vector of
% length 6(as long as y is a column vector).
[t,y] = ode23t(f,tspan,y0);
traiettoria= [traiettoria;y];
tspan = dt+tspan;
y0=y(end,:);
z=y(end,3);
end
figure
plot3 ( traiettoria(:,3), traiettoria(:,2), traiettoria(:,1))
grid on
xlabel ('direzione del fascio');
zlabel ('deflessione elettrica');
ylabel ('deflessione magnetica');
La mia idea era quella di fargli risolvere l'equazione in funzione del tempo ma con valori di B ed E differenti a seconda della zona in cui si trova la particella.
Per fare questo volevo fargli risolvere l'equazione per in un ciclo while (con un if-else si pianta) in cui verifica che la coordinata z = y(end,3) sia più grande o più piccola di un certo valore.
Nella mia testa fino a che è verificata la condizione del while dovrebbe rientrare nel ciclo, risolvere l'equazione incrementare il tempo aggiornare le condizioni iniziali e la z e verificare la condizione. Se è vera ripete se è falsa passa all'altro while in cui dovrebbe rifare le stesse cose ma con un B diverso... poi ci andrebbe un altro while che fa la stessa cosa dove B ed E sono diversi da zero e così via.
Il punto è che il primo ciclo sembra lo faccia correttamente ma se entra nel secondo while mi sa che non riesce più ad uscirne perchè ho aspettato un'ora e passa ma mi continua a comparire la scritta "busy" vicino allo "start"....