Dear Honghu,

You should clarify the notation (I don’t see any x in your f(x),

parameters appear loosely and normalization is not clear), and

exactly what you want to sample. It looks like you want to sample

from an exponential distribution with parameter T:

```
f(x) = (1/T) exp(-x/T).
```

Forget the approach in your source.f (there’s an underlying confusion

between tabulating a probability density function and sampling from

it). Sampling from an exponential distribution is analytical, so there’s

no need to tabulate anything. This distribution is normalized as it

stands:

```
\int_0^\infty dx (1/T) exp(-x/T) = 1
```

And moreover you can solve the inverse sampling problem analytically.

Let xi be a uniformly distributed random number in [0,1) i.e. what you

get from FLRNDM(XDUM). You then need to solve:

```
\int_0^x dx' (1/T) exp(-x'/T) = xi
```

This yields (double check it)

```
x = - T * ln(1-xi)
```

Or, since 1-xi and xi are both distributed homogeneously, you may take

```
x = -T * ln (xi)
```

Your x appears to be an energy, so you may want to discard sampled

values outside of a certain domain [E1,E2]. Try along the following

lines (maybe more elegant with a while loop but it does not matter

now):

```
E1 = <insert your lower energy bound here>
E2 = <insert your upper energy bound here>
666 XI = FLRNDM(XDUM)
XE = - T * LOG(XI)
IF (XE.LT.E1 .OR. XE.GT.E2) GO TO 666
```

This should work, but double check the math. It is more

expedient/cleaner/efficient to do the little bit of algebra above

restricting x in [E1,E2] from the onset (which affects both the

normalization and the inverse sampling equation) so as not

to waste time rejecting sampled values (depending on what

your E1 and E2 are this can be critical).

With kind regards,

Cesc