Monte Carlo electron scattering problem

Hi;

I'm trying to do a Monte Carlo electron scattering problem. Briefly what I have to do is start with a number of electrons (say 100), test the probability that each one scatters into one state or another (by testing a random number) and then place the electrons in the proper states. I have to do this for as many time intervals (of 1 picosecond) as it takes to get the initial number to get to zero.

So, No is the initial number of electrons, Nf is one scattering state based on the conditions in the code, and Ng is the other state. Here's the code I have so far:

%Input the lifetimes and time interval

tau1 = 30*(10^(-12));

tau2 = 60*(10^(-12));

deltatau = 1*(10^(-12));

%set the initial number of electrons 

No = 100;

%initialize Ni,Nf and Ng

Nf = 0;
Ng = 0;
Ni = No;

%calculate probabilities for the first time interval

P1 = (deltat/tau1);

P2 = (deltat/tau2);

Psum = P1+P2;

%generate a random number
r = rand(10,1)

%for the first time interval test if the first electron scatters

for i=1

      if (r(i)>0 & r(i)<P1)

     %if r is between 0 and P1 add 1 to Nf and subtract 1 from Ni

         Ni(i) = Ni-1

         Nf(i) = Nf+1

      end

     %if r is between P1 and P1+P2, add 1 to Ng and subtract 1 from Ni

     if (r(i)>P1 & r(i)<Psum)

         Ni(i) = Ni-1

         Ng(i) = Ng+1

     end

end

The thing I'm struggling with is how to loop through the values (hence why I only have i=1 in the for loop). So, if the first random number meets the condition, then I need to add 1 to either Nf or Ng, subtract 1 from Ni (so I'll have 99 Ni's left) and then go back to the beginning, look at the next random number and if that one meets the condition, add another value to Nf or Ng and subtract another value from Ni, until the value of Ni reaches zero.

It's a pretty convoluted problem, but it's basically a for loop within a for loop (I think) and I just can't set it up (my Matlab skills aren't that good). Any help is very much appreciated! Thanks!

Answer 1

[Nick]

Try this:

Ni=100;
Ng=0;
Nf=0;
Ni_mc=[];
Ng_mc=[];
Nf_mc=[];
i=0;

while Ni>0

    r=rand;
    i=i+1;

    if r>0 && r<P1
        Ni=Ni-1;
        Ni_mc(i)=Ni;
        Nf=Nf+1;
        Nf_mc(i)=Nf;

        Ng_mc(i)=Ng;
    elseif r>=P1 && r<Psum
        Ni=Ni-1;
        Ni_mc(i)=Ni;
        Ng=Ng+1;
        Ng_mc(i)=Ng;

        Nf_mc(i)=Nf;
    else
        Ni_mc(end+1)=Ni;
        Nf_mc(end+1)=Nf;
        Ng_mc(end+1)=Ng;
    end

end

Im not sure what you want to do if r=P1, I just set the electron to the Ng state if that happens, you can change it if you want.

The final value of "i" should tell you how many tries it took to get Ni down to 0.

[Nick]

Wait, before you run this script, I just realized that given how low your probabilities are, its going to take a long time before Ni gets to 0. Does the problem deal with more than 100 electrons? If the initial electron number is really high then we might have to do something else.