giug ha scritto:
Continua a non essermi chiara la domanda. più che altro non ho capito, descrivendola come fosse una funzione, cosa vuoi che abbia in ingresso e cosa in uscita. Tu sai che si è spostata di un metro e in uscita vuoi che ti metta l'1 al terzo elemento del vettore? o deve trovare quelle che tu hai messo come punti interrogativi nell'y_diaframma?
In ingresso il vettore che uso come condizione iniziale, in uscita deve trovare le coordinate x e y (ossia i punti interrogativi) della particella una volta raggiunto il diaframma.
Praticamente fargli fare quello che fa l'ode solver senza odesolver.
giug ha scritto:
k è il numero di particelle sopravvissute (+1 dato che l'ultima volta che esce dal ciclo aumenta di 1 anche se ha finito). Dovrebbe essere esattamente uguale alla dimensione della struttura +1
Lo pensavo pure io ma non è così.
o probabilmente fa qualcosa che non capisco perchè se nell'if gli faccio fare pure una matrice
T(:,:,i)=ys;
è di dimensioni 2x6x794 e anche la struttura è 794 ma è tutta fatta di zeri...
Ti dò il codice:
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);
%spread angolare
A = -20; B = 20;
angolo(i) = A + (B-A) * rand(1); %random number in [A,B]
teta=angolo(i); %angolo nel piano zy
fi=angolo(i); %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];
% 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)
Particella(i).traiettoria=ys(:,1:3);
Particella(i).velocita=ys(:,4:6);
k=k+1;
end
end