Montecarlo method for particle simulation

12 次查看（过去 30 天）
Francesco
Francesco 2012-2-8
Hi all,
I have to simulate a number of particle with a initial random velocity along z axis emitted by a point source and distributed inside a solid angle passing through a pinhole. I would like to "kill" particles that cannot pass through the pinhole becuase of the angular spread and mantain those that can pass.
Is there any example I can use to write this simulation?

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2012-2-8
It's really not that hard. 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. It's like 7-10 lines or so.

更多回答（3 个）

Francesco
Francesco 2012-2-9
I was thinking to do something similar. My concern is that I have two pinholes (the second one is 10cm further) Do I need to solve the motion equation from the source to the first pinhole and then do the same between the pinholes? Or is it useless?
  1 个评论
Image Analyst
Image Analyst 2012-2-9
You need to figure out what the distance is from the start, and where that particular particle would hit laterally the blocking screen/aperture plane, and decide if that location is within the pinhole locations or not.

请先登录，再进行评论。

Francesco
Francesco 2012-2-10
I think I'm making the whole thing too complicate... I give an initial energy Ein and n=1000 particles.
the for loop is:
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 = -20; B = 20;
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);
%loop to select particles
if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1)
T(:,:,i)=ys;
Particella(i).traiettoria=ys(:,1:3);
Particella(i).velocita=ys(:,4:6);
k=k+1;
end
end
But I'm not sure of the result.
  2 个评论
Image Analyst
Image Analyst 2012-2-10
I don't know. It's hard for me to follow since I don't understand Italian.
Francesco
Francesco 2012-2-10
Sorry, I forgot to translate my comments... I'm fixing it

请先登录，再进行评论。

Francesco
Francesco 2012-2-10
My beam is moving along z axis
for i=1:n
K_ion=normrnd(Ein,Ein/10); %energy distribution
Etot_ion(i)=K_ion+((m_ion*uma1)/uma); % Total Energy in MeV
betasquare_ion(i)=1-(((m_ion*uma1)/(uma*Etot_ion(i)))^2);
v_ion(i)=sqrt(betasquare_ion(i)*c^2); %Particle velocity
%%angular spread
A = -20; B = 20;
angolo = A + (B-A) * rand(1); %random number in [A,B]
teta=angolo; %angle between velocity zy plane
fi=angolo; % angle between velocity x
%%%%projection of v along z
v_z= v_ion*cosd(teta)*cosd(fi);
%%%%%projection of v along x
v_x= v_ion*sind(teta);
%%%%%projection of v along 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 conditions (position and velocity at the source)
tspan=[0,100];
O=[0 0 0]';
%%%solve the motion equation (drift)
f = @(t,ys) [ys(4:6); O];
options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1)); % options
%set error ---% Event_Stop_Sorgente stops ode solver when particles reaches the pinhole
[t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
%loop to select particles
if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1) %if x or y coordinates of a particle is smaller of the pinhole radius save particle, esle the particle is ignored
T(:,:,i)=ys; %matrix containing particles passing through the pinhole
Particella(i).traiettoria=ys(:,1:3); %structure containg particle
Particella(i).velocita=ys(:,4:6);
k=k+1; %particle counter
end
end
I'm using a matrix and a structer at the moment because I'm structures are really new for me and I use both to be safe, I'm going to delet the matrix.
I think you need the Event_Stop function:
function [value,isterminal,direction]=Event_Stop_Sorgente(t,ys, Pinhole1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Function stopping ode solver. % In the 'events' case, the ODE file returns to the solver the values % that it needs to perform event location. When the Events property is set % to on, the ODE solvers examine any elements of the event vector for % transitions to, from, or through zero. If the corresponding element of % the logical isterminal vector is set to 1, integration will halt when a % zero-crossing is detected. The elements of the direction vector are -1, % 1, or 0, specifying that the corresponding event must be decreasing, % increasing, or that any crossing is to be detected.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % When "value" is equal to zero, an event is triggered.
% % set isterminal to 1 to stop the solver at the first event, or 0 to
% % get all the events.
% % direction=0 if all zeros are to be computed (the default), +1 if
% % only zeros where the event function is increasing, and -1 if only
% % zeros where the event function is decreasing.
% value = ......; % when value = 0, an event is triggered
% isterminal = 1; % terminate after the first event
% direction = 0; % get all the zeros
%Pinhole1=1; distance between source and pinhole
z=ys(3); %particle position z
cond_z=z-Pinhole1; %%%zero passing condition
value=cond_z; % event function
isterminal=1; %%%Stop the integration when the particle reach the pinhole
direction=0; %%%Stop the integration when the zero crossing occurs from every direction.

类别

在 Help Center 和 File Exchange 中查找有关 General Applications 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by