## 2 Calculation Model
Moreover, the diffuse reflected part is stored as the starting energy of the diffuse energy transfer process, which can be expressed as follows:  
\[E_{patch}=E_{received}\cdot(1-\alpha_{i})\cdot s_{i} \tag{4}\]  
With respect to the directional distribution of the scattering energy, it was assumed that the diffuse reflected sound energy follows the Lambert cosine law, which has been widely used in different models [15; 17; 19; 37]. With the Lambert cosine law, the sound energy transfer between two patches is independent of the incident angle, which could significantly improve the time efficiency. As shown in Figure 2, the sound energy transfer from \(Patch_{i}\) to \(Patch_{i}\), namely Form Factor, can be calculated as follows:  
\[E_{received_{j}}=\int\limits_{S_{i}}\int\limits_{S_{j}}\frac{cos\theta_{i}cos \theta_{j}}{\pi{d_{ij}}^{2}}dS_{i}dS_{j}\cdot E_{patch_{i}} \tag{5}\]  
where, \(E_{received_{j}}\) and \(E_{patch_{i}}\) are the energy received by \(Patch_{j}\) and the energy emitted from \(Patch_{i}\), respectively. The sound energy emitted from an image patch source can be calculated as follows:  
\[E_{image\ patch_{i,n}}=E_{patch_{i}}\prod_{k=1}^{n}(1-\alpha_{k})\cdot(1-s_{k}) \tag{6}\]  
Figure 2: Sound energy transfer from sources to patches. In the figure above, \(d\) and \(\theta\) represent the distance and the incident angle from the source to the patch, respectively. where, \(E_{image\ patch_{i,n}}\) and \(E_{patch_{i}}\) are the energy of the nth order image patch and that of the patch source, respectively; and \(\alpha_{k}\) and \(s_{k}\) are the absorption coefficient and scattering coefficient in the \(k\)th order reflection, respectively,  
As shown in Figure 2, the sound energy transfer from an image patch to the patches can be calculated as follows:  
\[E_{received_{i}}=\int\limits_{S_{r}}\int\limits_{S_{j}}\frac{cos\theta_{i^{ \prime}}\cdot cos\theta_{j^{\prime}}}{\pi d_{r^{\prime}}}dS_{i^{\prime}}dS_{i^{ \prime}}E_{image\ patch_{i,n}} \tag{7}\]  
All the sound energy received by the receivers is either from point sources or surface sources. For the point source, which includes the original sound source and image source, the sound intensity at a receiver position can be calculated using Equation (2). The sound intensity due to a patch source can be calculated as follows:  
Figure 1:0  
where, \(E_{patch}\) and \(S_{i}\) are the energy and the area of the patch source, respectively; and \(\theta\) and \(d\) are the incident angle and the distance between the receiver and a point on the surface source, respectively. Ideally, the accuracy will improve when higher-order reflections and image patches are calculated. However, with the increase in the iteration number, the sound energy calculated will rapidly decrease and the average propagation time will increase. Therefore, the simulation accuracy will not significantly improve by taking into account more reflections when the considered reflection order is very high. Meanwhile, the exponential growth of the number of energy transfer paths can be found with the increase in the number of iterations. Hence, a termination condition is necessary to obtain a balance of accuracy and time costs. In IMCRI, an auto-termination process was applied to stop creating new image patches and the whole sound energy transfer process. For each iteration, the sound pressure level (SPL) and reverberation time (RT) of every receiver position were calculated with the received impulse response. The calculation process continued until the change in the SPL and RT caused by an entire iteration was less than a termination threshold. The termination threshold (TT) was defined as the minimum value of acoustic parameters to continue the iteration process. TT = 0.01 means the calculation program will stop when the change of SPL and RT in the entire iteration is less than 0.01 dB and 0.01 s, respectively. As shown in Figure 3, the RT simulation was more sensitive to TT because it was mainly determined by the late reflected sound.