The general form for this equation is:

vtrap = x/(exp(x/y) - 1)

however as can be easily seen if x/y is 0, or close to zero, then the denominator of the given equation is really small or zero, leading to and infinite or very large output. To get around this when x/y is really small the output is given as:

vtrap = y*(1 - x/y/2) or as I prefer y - x/2

hence the following code:

Code: Select all

```
FUNCTION vtrap(x,y) { :Traps for 0 in denominator of rate eqns.
if (fabs(x/y) < 1e-6) {
vtrap = y*(1 - x/y/2)
}else{
vtrap = x/(exp(x/y) - 1)
}
}
```

To start exp(x/y) can be approximated by (where u = x/y) 1 + u + u^2/2 + more terms which make the approximation more accurate which I will not use (second order approximation). Thus from the first equation we have:

x/(u+u^2/2) => the 1's canceled

Going back to x/y the equation can be simplified fairly easily to:

y^2/(y+x/2)

however this still contains a denominator, which I guess would be nice to get rid of (??) or perhaps the squared term for y isn't stable (??) (if we assume that y may be very large to cause x/y to be small)

Regardless, the following simplification can be made:

y^2 * (y-x/2)

-------------------------

(y+x/2) * (y-x/2)

This leads to:

y^3 - y^2*x/2

-----------------

y^2 - x^2/4

If x^2/4 is assumed to be small, or I guess small compared to y^2, the denominator can be approximated as y^2 (this is the second and final approximation).

This simplifies the original equation down to y - x/2