source (or a group of radiation sources), $\alpha$, and $i = {\ensuremath{\mbox{{H$\,${\sc i}}}}}$, ${\mbox{{He$\,${\sc i}}}}$, and ${\mbox{{He$\,${\sc ii}}}}$. The contribution from a point-like radiation source, $\alpha$, is represented by $$\Gamma_{i, \alpha} = \frac{1}{4 \pi h r_\alpha^2} \int_{\nu_i}^\infty 
  \frac{{{\rm d}}\nu}{\nu} \sigma_i(\nu) L_\alpha(\nu) 
  \exp\left[-\sum_j {\cal N}_{j, \alpha} \sigma_j (\nu)\right], 
\label{eq:photoionization}$$ where $\nu_{i}$ is the threshold frequency for the $i$-th species, $\sigma_i(\nu)$ is the cross-section of the $i$-th species, and $r_\alpha$, $L_\alpha(\nu)$, and ${\cal N}_{i, \alpha}$ are respectively the distance between the luminosity centre and the target mesh, the intrinsic luminosity of the radiation source (or the group of the sources), and the column density of $i$-th species. The sum in the exponent runs over all three chemical species. When all sources have the same spectral shape, i.e. $L_\alpha(\nu) = C_\alpha f(\nu)$, we generate a look-up table for each species as a function of column densities: $$g_i({\cal N}_k) = \int_{\nu_i}^{\infty}   
\frac{{{\rm d}}\nu}{\nu} \sigma_i(\nu) f(\nu) 
  \exp\left[-\sum_j {\cal N}_j \sigma_j (\nu)\right]. 
\label{eq:tables}$$ In our case, a look-up table for each chemical species becomes three-dimensional table. We have confirmed that 20 logarithmic bins for each column density is sufficient. By using the look-up tables, the RT calculation is reduced to evaluating the column densities. Following @anninos97, 