Montecarlo method for particle simulation

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?

Answer 1

[Image Analyst]

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.

Answer 2

[Francesco]

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?

[Image Analyst]

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.

Answer 3

[Francesco]

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.

[Image Analyst]

I don't know. It's hard for me to follow since I don't understand Italian.

[Francesco]

Sorry, I forgot to translate my comments... I'm fixing it

Answer 4

[Francesco]

 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.