Se scrivo:
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=0.03;
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);
%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);
v_z=v_ion;
%calcolo la componente x
v_x= unifrnd(-0.7875,0.7875, n,1); %angolo compreso da -45 e +45
v_y= unifrnd(-0.7875,0.7875, n,1);
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];
Pinhole1=1;
DiametroPinhole1=0.002;
r_pin1=DiametroPinhole1/2;
O=[0 0 0]';
%%%Risolvo il moto
f = @(t,ys) [ys(4:6); O];
% the expression [y(4:6); (q/m)*cross(y(4:6),B)]
options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1));
[t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
T(:,:,i)=ys;
end
va bene?
T(:,:,i)=ys è la matrice che mi serve?