giug ha scritto:
E porti fuori anche tutta la parte sotto fino all'ode.
Sì, per gli indici hai ragione, (i) l'avevo dimenticato, (1,n) non lo so perchè l'ho messo...
Per il resto aspetta un attimo. Il codice ora è cambiato, dentro al ciclo interno c'è tutta la storia degli angoli...
for s=1:stati
%Ein_s=Ein;
Ein_s=s*Ein;
m_ions=m_ion-(s*m_e);
K_ion=normrnd(Ein_s,Ein_s/10);
Etot_ion=K_ion+((m_ions*uma1)/uma);
betasquare_ion=1-(((m_ions*uma1)/(uma.*Etot_ion)).^2);
v_ion=sqrt(betasquare_ion.*c^2);
for i=1:n
v0=feval(@God_Speed,n,ang_spread,v_ion);
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);
%oppure y_sorgente
%decido quali particelle passano il pimo diaframma
if (abs(ys(end,1)) <=r_pin1 && abs(ys(end,2))<=r_pin1) %%????? LA CONDIZIONE DEVE ESSERE && o || ??????
%T(:,:,k)=ys;
Particella(k).ionizzazione=s;
Particella(k).stato_di_carica=s*q;
Particella(k).massa=m_ions;
Particella(k).q_over_m= s*q/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