Simulating Two-Level Electron Scattering

[quark.antiquark]

Homework Statement

I have to simulate a two-level electron scattering process. 1000 electrons start in state i and over time, end up in state f. I need to plot the number of electrons in each state as a function of time.

Homework Equations

None.

The Attempt at a Solution

The probability of going from state i to f is given. I had previously generated an array of random numbers between 0 and 1 (called xArray in the code below). At each time interval, I pick a random number and compare it to the probability. If the number is less than the probability, the electron goes to state f. Otherwise, it stays in state i. I do this for each electron in state i. Then, in the next time interval, I repeat the process, and so on. This all looks good on paper, but I'm not exactly the world's best programmer. My program seems to have a bug, and I've been trying to find it all evening. I end up with a negative number of electrons in the initial state, which doesn't make any sense. I've attached the relevant parts of my code (java); extraneous stuff like the random number generator have been left out for brevity. Hopefully someone will be able to spot the bug!

final double pf = 1./30.; // Probability of going to state f
int[] nfArray = new int[100000]; // Number in state f as a function of t
nfArray[0] = 0; // Start with 0 electrons in state f
int[] niArray = new int[100000]; // Number in state i as a function of t
niArray[0] = 1000; // Start with 1000 electrons in state i
int ni = 1000; // Number of electrons in state i at time t
int[] tArray = new int[100000]; // Elapsed time (for plotting purposes)
tArray[0] = 0;
int t = 1; // Time in ps
int j = 0; // Index of xArray
int k = 0; // # of electrons that have decayed (or not) at
// time t
while(ni > 0) // Stop when no electrons remain in state i
{
k = 0;

while(k < ni) // Loop through all electrons in state i
{
if(xArray[j] < pf) // If true, the electron goes to state f
{
nfArray[t] = nfArray[t] + 1; // One electron goes to state f
niArray[t] = niArray[t] - 1; // Also decrease the number in state i
}

else
{
; // The electron stays in state i
}

j++; // Pick the next random number
k++; // Go to the next electron in the system
}

ni = niArray[t]; // Number of electrons remaining in state i after time t
tArray[t] = t;
t++; // Next time interval
}

System.out.println("Time elapsed: " + tArray[t-1]);
System.out.println("Number in state i: " + niArray[t-1]);
System.out.println("Number in state f: " + nfArray[t-1]);

 

[dotman]

Hello,

I don't see where you're initializing the values of the ni and nf arrays, other than for index = 0, and I think this may be your problem. Either loop through them after initialization:

int[] nfArray = new int[100000]; // Number in state f as a function of t
for (int i = 0; i < 100000; i++)
nfArray = 0;

Or just initialize each one in your while loop as you come to it:

while(ni > 0) // Stop when no electrons remain in state i
{
k = 0;
niArray[t] = ni;
if (t != 0) nfArray[t] = nfArray[t-1];

Or something. I think this is your problem. If not, maybe try adding some System.out.printlns throughout the code, looking for funky values-- this usually helps me.

Hope this helps!

 

[quark.antiquark]

I already checked that out... the array elements initialize to 0 by default.

 

[dotman]

Hello,

Well, its still good form to initialize them-- that way you -know- they're 0.

Another thing is I don't see where you're initializing nfArray and niArray after each round-- ie, if 30 flip in round 1, you'll have:

nfArray[0] is 30
niArray[0] is 970

but

nfArray[1] is 0
niArray[1] is 0

So while you're loop will operate right (because it runs on ni, not niArray), the array value is off.

Could this be it? Try the printlns :-)

 

[quark.antiquark]

Haha, yes, I just realized that as soon as I replied to your post. I've modified my code (below) but I still have the same problem. I'll keep tinkering...

while(ni > 0) // Stop when no electrons remain in state i
{
k = 0;

while(k < ni) // Loop through all electrons in state i
{
if(xArray[j] < pf) // If true, the electron goes to state f
{
nfArray[t] = nfArray[t] + 1; // One electron goes to state f
niArray[t] = niArray[t] - 1; // Also decrease the number in state i
}

else
{
; // The electron stays in state i
}

j++; // Pick the next random number
k++; // Go to the next electron in the system
}

ni = niArray[t]; // Number of electrons remaining in state i after time t
tArray[t] = t;
t++; // Next time interval
niArray[t] = ni; // # of electrons starting out in state i at time t
nfArray[t] = nfArray[t-1]; // # starting out in state f at time t
}

 

[quark.antiquark]

If it's any help, my simulation ends at t = 1 which definitely should not be happening

 

[dotman]

Well, if its ending at t=1, that's telling you that for some reason ni is less than or equal to zero after the initial iteration. Not quite seeing where exactly this is happening, based on what you're showing. Try some printlns would be my suggestion, especially inside the loops.

I also wanted to note that, if you've got the time and you're up for some simplification, you can simplify this business a bit with something like

niArray[0] = 1000;
nfArray[0] = 0;
int t = 0;
int flipcount;

while (niArray[t] > 0)
{
t++;
niArray[t] = niArray[t-1]; // Initialize
nfArray[t] = nfArray[t-1]; // Initialize

flipcount = 0;
randnum = getRandNum();
for (int i = 0; i < niArray[t]; i++)
{
if (randnum < pf)
{
flipcount++;
}
}

niArray[t] -= flipcount;
nfArray[t] += flipcount;
}

Or something, there may be a bug or two in there. Basically, you have more indicies that your going over (ni, k, etc) than you really need, if you want to simplify.

 

[quark.antiquark]

You, sir (or madam), are my hero. I don't know what was wrong with my code, but when I used your simplified version and plotted the results, I got the most beautiful pair of exponential curves I have ever seen!

Now I'll be awake for the rest of the night trying to figure out what was wrong with my code...

Thanks!