Add script with correct theoretical covariance matrix calculation

[Bai-Chiang]

If estimated power spectra is average over all cross spectra across all splits, the covariance matrix calculation was not correct.

Some math:

If the estimated power spectrum of frequency $f_i$ polarization $p_i$ cross frequency $f_j$ polarization $p_j$ is average over all cross spectra

\hat{C}_l^{f_ip_i \times f_j p_j} = \frac{1}{\text{num of cross pairs}} \sum_{\text{all cross splits } (s_{\alpha}, s_{\beta})} \hat{C}_l^{s_{\alpha}f_ip_i \times s_{\beta}f_j p_j}

Then the theoretical prediction of the covariance matrix is (ignoring $f_{sky}$ and binning factor)

< ( \hat{C}_l^{f_ip_i \times f_j p_j} - C_l^{f_ip_i \times f_j p_j} ) ( \hat{C}_{l'}^{f_k p_k \times f_l p_l} - C_{l'}^{f_k p_k \times f_l p_l} ) >

= \frac{\delta_{ll'}}{2l + 1} \frac{1}{\text{num of } (s_{\alpha}, s_{\beta})} \frac{1}{\text{num of } (s_{\mu}, s_{\nu})} \sum_{(s_{\alpha}, s_{\beta})} \sum_{(s_{\mu}, s_{\nu})} \Big[
<\hat{C}_l^{s_{\alpha}f_ip_i \times s_{\mu}f_k p_k}> <\hat{C}_l^{s_{\beta}f_jp_j \times s_{\nu}f_l p_l}> + <\hat{C}_l^{s_{\alpha}f_ip_i \times s_{\nu}f_l p_l}> <\hat{C}_l^{s_{\beta}f_jp_j \times s_{\mu}f_k p_k}>
\Big]

where $(s_{\alpha},s_{\beta})$ are all cross splits of $\hat{C}_l^{s_{\alpha}f_ip_i \times s_{\beta}f_j p_j}$, and $(s_{\mu},s_{\mu})$ are all cross splits of $\hat{C}_l^{s_{\mu}f_kp_k \times s_{\nu}f_l p_l}$

The point is that some of the terms $<\hat{C}_l^{s_{\alpha}f_i p_i \times s_{\mu}f_k p_k}>$ would have noise power, while others only have signal. The old calculation overestimated the covariance, the new one should decrease std in final result, which is what I found when testing with main branch code.

This commit simply loop over two summation $\sum_{(s_{\alpha}, s_{\beta})} \sum_{(s_{\mu}, s_{\nu})}$ and accumulate each term, since I find it's hard to calculate how many terms have (both or one or none) noise spectrum.