Pare che ho finito.
Gli faccio fare il passaggio dai due collimatori con un angolo ragionevole e perdo poco meno di 400 particelle.
Puoi dare un'occhiata al codice e vedi che te ne pare.
Ho utilizzato le strutture anche se mi sono fatto anche le matrici per sicurezza... quando sarò sicuro di non doverle usare le levo.
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=1000; %numero diparticelle
Ein=3;
m_ion=1.00794*uma;
q_over_m=q/m_ion;
Etot_ion=zeros(n,1);
betasquare_ion=zeros(n,1);
v_ion=zeros(n,1);
%%
%Passaggio dal primo pinhole
%%%%%%%%%%%%
Pinhole1=1;
DiametroPinhole1=0.002;
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 i=1:n
K_ion=normrnd(Ein,Ein/10);
Etot_ion(i)=K_ion+((m_ion*uma1)/uma); % ho tutto in MeV
betasquare_ion(i)=1-(((m_ion*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;
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 speed
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).traiettoria=ys(:,1:3);
Particella(k).velocita=ys(:,4:6);
k=k+1;
end
end
%%
%Passaggio dal secondo pinhole
Pinhole2=0.1;
Distanzap1p2=Pinhole1+Pinhole2;
DiametroPinhole2=0.002;
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
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)
a=1;
T2(:,:,j)=yp;
Part_collimata(j).traiettoria=yp(:,1:3);
Part_collimata(j).velocita=yp(:,4:6);
j=j+1;
end
end