La storia era che mi servono le componenti della velocità lungo i tre assi, non una velocità lungo un asse e due angoli per le altre due componenti, che è una cosa senza senso....
Allora, 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=100;
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);
%%%%%%%%%%%%
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.
%T=zeros(2,6,n); % matrice tridimensionale dove immagazino le traiettore
%calcolo la componente z della velocità
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);
%spread angolare
A = -45; B = 45;
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(fi)*sind(teta);
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];
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);
%struttura
Particella(i).traiettoria=ys(:,1:3);
Particella(i).velocita=ys(:,4:6);
end
per ammazzare le particelle gli faccio fare un altro ciclo:
for i=1:n
if ( abs(Particella(i).traiettoria(end,1))>= r_pin1 || abs(Particella(i).traiettoria(end,1))>= r_pin1)
Particella(i)=[];
n=n-1;
i=n-1;
end
end
e un poco le ammazza perchè la struttura Particella (senti che termine da uno che ne capisce...) esce dal primo ciclo che è 1x100
ma ad un certo punto ho un errore nel secondo ciclo:
??? Index exceeds matrix dimensions.
Error in ==> Prova at 99
if ( abs(Particella(i).traiettoria(end,1))>= r_pin1 ||
abs(Particella(i).traiettoria(end,1))>= r_pin1)
e la struttura si è ridotta a 1x50 ma lui si incazza (come sempre) perchè l'indice diventa più grosso della matrice... io però ho copiato il codice tuo....(e questo dimostra che non ci capisco nulla)
Mi dici anche se con le strutture faccio giusto?