MATLAB Answers

Montecarlo method for particle simulation

36 views (last 30 days)
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?

Accepted Answer

Image Analyst
Image Analyst on 8 Feb 2012
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.

More Answers (3)

Francesco
Francesco on 9 Feb 2012
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 Comment
Image Analyst
Image Analyst on 9 Feb 2012
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.

Sign in to comment.


Francesco
Francesco on 10 Feb 2012
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 Comments
Francesco
Francesco on 10 Feb 2012
Sorry, I forgot to translate my comments... I'm fixing it

Sign in to comment.


Francesco
Francesco on 10 Feb 2012
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.

Community Treasure Hunt

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

Start Hunting!

Translated by