Allora, il primo pezzo l'ho fatto e sembra funzionare.
Ho fatto così:
clear all
close all
uma=1.66e-27; %[Kg], unità di massa atomica
q=1.602e-19; %[C], carica elementare
m_e=5.485799e-4*uma; %[kg], massa dell'elettrone
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;
n=500; %numero di particelle per stato di carica
Ein=0.3;
m_ion=1.00794*uma;
%q_over_m=q/m_ion;
stati= 3; %numero stati di carica
Etot_ion=zeros(n,1);
betasquare_ion=zeros(n,1);
v_ion=zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%
%COLLIMAZIONE FASCIO
%
%%
%%
%Passaggio dal primo pinhole
%%%%%%%%%%%%
Pinhole1=1;
DiametroPinhole1=0.001;
r_pin1=DiametroPinhole1/2;
%Have a for loop over a bunch of particles. For each particle, set its
%velocity (speed and angle) vector using the rand() function. Then use an
%"if" statement to decide whether the particle passed through or not and
%count those that do.
k=1; %indice di supporto per selezionare le particelle nel ciclo interno if
for s=1:stati
Ein_s=s*Ein;
m_ions=m_ion-(s*m_e);
for i=1:n
K_ion=normrnd(Ein_s,Ein_s/10);
Etot_ion(i)=K_ion+((m_ions*uma1)/uma); % ho tutto in MeV
betasquare_ion(i)=1-(((m_ions*uma1)/(uma*Etot_ion(i)))^2);
v_ion(i)=sqrt(betasquare_ion(i)*c^2);
%angular spread
%A = -0.0575; B = 0.0575;
A=-0.08; B=0.08;
%A=0; B=0;
angolo = A + (B-A) * rand(1); %random number in [A,B]
teta=angolo; %angolo nel piano zy
fi=angolo; %angolo con l'asse x
%%%%Componente z della velocità
v_z= v_ion*cosd(teta)*cosd(fi);
%calcolo la componente x
v_x= v_ion*sind(teta);
%calcolo la componente y
v_y=v_ion*cosd(teta)*sind(fi);
v0=[v_x, v_y, v_z]; %initial velocity
r0=zeros(n,3); %initial position
y_sorgente=[r0,v0]; %vector of initial condition
tspan=[0,100];
O=[0 0 0]';
%%%Risolvo il moto
f = @(t,ys) [ys(4:6); O];
% the expression [y(4:6); O]
% combines two vectors of length 3 to make a vector of
% length 6(as long as y is a column vector).
options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1)); % opzione per risoluzione equazione differenziale:
%set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
[t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
%decido quali particelle passano il pimo diaframma
if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1)
% T(:,:,k)=ys;
Particella(k).stato_di_carica=s;
Particella(k).massa=m_ions;
Particella(k).q_over_m= s/m_ions;
Particella(k).energia_iniziale=Ein_s;
Particella(k).traiettoria=ys(:,1:3);
Particella(k).velocita=ys(:,4:6);
k=k+1;
end
end
end
%%
%Passaggio dal secondo pinhole
Pinhole2=0.1;
Distanzap1p2=Pinhole1+Pinhole2;
DiametroPinhole2=0.0005;
r_pin2=DiametroPinhole2/2;
m=k-1; %numero di particelle che passano
y_pin2=zeros(m,6); %%%matrice condizione iniziale per il passaggio dal secondo pinhole
% for i=1:m
%
% %y_pin2(i,:)=[Particella(i).traiettoria(end,:),Particella(i).velocita(end,:)]; %vector of initial condition
% %y_pin2(i,:)=T(end,:,i);
% end
j=1; %
for i=1:m
y_pin2(i,:)=[Particella(i).traiettoria(end,:),Particella(i).velocita(end,:)]; %vector of initial condition
tspan=[0,100];
%O=[0 0 0]';
%%%Risolvo il moto
f = @(t,yp) [yp(4:6); O];
% the expression [y(4:6); O]
% combines two vectors of length 3 to make a vector of
% length 6(as long as y is a column vector).
options2=odeset('RelTol',1e-7,'Events',@(t,yp)Event_Stop_Coll(t,yp,Distanzap1p2)); % opzione per risoluzione equazione differenziale:
%set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
[t,yp] = ode23t(f,tspan,y_pin2(i,:),options2);
%decido quali particelle passano il pimo diaframma
if (abs(yp(end,1)) <=r_pin2 || abs(yp(end,2))<=r_pin2)
% T2(:,:,j)=yp;
Part_collimata(j).stato_di_carica =Particella(j).stato_di_carica;
Part_collimata(j).massa = Particella(j).massa;
Part_collimata(j).q_over_m= Particella(j).q_over_m;
Part_collimata(j).traiettoria=yp(:,1:3);
Part_collimata(j).velocita=yp(:,4:6);
j=j+1;
end
end
%number of particles entering the spectrometer
p=j-1;
figure('color', 'white')
hold all
for j=1:m
plot3(Particella(j).traiettoria(:,3),Particella(j).traiettoria(:,2), Particella(j).traiettoria(:,1))
%F(j)= getframe;
end
for i=1:p
plot3(Part_collimata(i).traiettoria(:,3),Part_collimata(i).traiettoria(:,2), Part_collimata(i).traiettoria(:,1))
Part_trasp(i) = Part_collimata(i); %particelle da trasportare
end
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%
%%TRASPORTO DENTRO LO SPETTROMETRO
%
%%
%%
%Distanze in metri
P2_B = 0.091; %distanza secondo pinhole campo magnetico (Deriva)
Deriva1=Distanzap1p2+P2_B; %distanza totale percorsa dalla particella prima del settore di deflessione
LB=0.15; %Lunghezza campo magnetico
P2_E = 0.191; %distanza secondo pinhole campo elettrico (solo campo magnetico)
inizioE=Distanzap1p2+P2_E; %distanza totale percorsa dalla particella prima dell'ingresso nel campo elettrico
P2_fineB=0.241; %distanza secondo pinhole fine campo magnetico (zona a due campi)
fineB= Distanzap1p2+P2_fineB; %distanza totale percorsa dalla particella fino all'uscita dal campo magnetico
LE=0.07; %Lunghezza campo elettrico
P2_fineE=0.261; %distanza secondo pinhole fine campo elettrico (solo campo elettrico)
fineE=Distanzap1p2+P2_fineE; %distanza totale percorsa dalla particella fino all'uscita dal campo elettrico
Deriva2=0.393; %Distanza spettrometro rivelatore
Dist_Spettr_riv=Distanzap1p2+Deriva2; %distanza totale percorsa dalla particella fino al rivelatore
%%
%Zona 1
%Calcolo traiettoria nel tratto di deriva tra collimatori e spettrometro
y0=zeros(p,6); %inizializzo la matrice delle condizioni iniziali
a=1;
for i=1:p
%Campi nulli
B=[0 0 0]';
E=[0 0 0]';
y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)]; %vettore delle condizini iniziali
%%%Risolvo il moto
f = @(t,y1) [y1(4:6); (Part_trasp(i).q_over_m).*cross(y1(4:6),B)+(Part_trasp(i).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).
options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_1(t,y,Deriva1)); % opzione per risoluzione equazione differenziale:
%set precisione calcolo---% Sintassi per utilizzare Event_Stop con Deriva1 come variabile in input
[t,y1] = ode23t(f,tspan,y0(i,:),options);
%traiettoria= y1; % aggiorno la traiettoria
Part_trasp(i).traiettoria(end+1:end+size(y1,1),1:3)=y1(:,1:3); %aggiungo alla traiettoria le nuove coordinate partendo dall'utlima riga
Part_trasp(i).velocita(end+1:end+size(y1,1),1:3)=y1(:,4:6); %aggiungo alla veloctà le nuove coordinate partendo dall'utlima riga
end
%PLOT
figure('color', 'white')
hold all
for i=1:p
plot3(Part_trasp(i).traiettoria(:,3),Part_trasp(i).traiettoria(:,2), Part_trasp(i).traiettoria(:,1))
%Part_trasp(i) = Part_collimata(i); %particelle da trasportare
end
grid on
Poi nella parte successiva si blocca: dice "busy" e non da segni di vita. pensavo si bloccasse nel ciclo ma anche se dico:
i=1;
B2=[0.003 0 0]';
E2=[0 0 0]';
y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)]; %vettore delle condizini iniziali
%%%Risolvo il moto
f = @(t,y2) [y2(4:6); (Part_trasp(i).q_over_m).*cross(y2(4:6),B2)+(Part_trasp(i).q_over_m).*E2];
% 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).
options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_2(t,y,inizioE)); % opzione per risoluzione equazione differenziale, set precisione calcolo
[t,y2] = ode23t(f,tspan,y0(i,:),options);
%traiettoria= [y1;y2]; % aggiorno la traiettoria
Part_trasp(i).traiettoria(end+1:end+size(y2,1),1:3)=y2(:,1:3); %aggiungo alla traiettoria le nuove coordinate partendo dall'utlima riga
Part_trasp(i).velocita(end+1:end+size(y2,1),1:3)=y2(:,4:6); %aggiungo alla veloctà le nuove coordinate partendo dall'utlima riga
lui non la smette più di fare conti, e non capisco perchè