Figure 8.8:

RTD p(\theta ) versus \theta for n CSTRs in series, \tau =2.

Code for Figure 8.8

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
# Converted from ncstrspdens.m - RTD density for n CSTRs in series
import numpy as np
from scipy.special import gammaln
from misc import octave_save

x   = np.linspace(0., 5., 200)
tau = 2.
ns  = [1, 2, 3, 5, 10, 25, 100]
p   = np.zeros((len(x), len(ns)))

for i, n in enumerate(ns):
    logp    = -gammaln(n) + n*np.log(n/tau) + (n-1)*np.log(np.maximum(x, 1e-300)) - n*x/tau
    p[:,i]  = np.exp(logp)

ptable = np.column_stack([x, p])
octave_save('ncstrspdens.dat', ('ptable', ptable))