data.py

import numpy as np

class BBdata:
    def __init__(self, version=1):
        self.version = version

    def get_bicep2021(self, symmetric=True):
        if self.version == 1:
            bicepdat = np.loadtxt('Data/v1/bicep.txt')
            ℓ, D_ℓ, σ_D_ℓ = bicepdat[:, 1], bicepdat[:, 6], bicepdat[:, 7]
            return ℓ, D_ℓ, σ_D_ℓ
        elif self.version == 2:
            bicepdat = np.loadtxt('Data/v2/bicep.txt')
            ℓ, D_ℓ, Db_low, Db_high = bicepdat[:, 1], bicepdat[:, 3], bicepdat[:, 4], bicepdat[:, 5]
            low, high = D_ℓ - Db_low, Db_high - D_ℓ
            if symmetric:
                σ_D_ℓ = np.sqrt(low**2 + high**2)
            else:
                σ_D_ℓ = np.stack([low, high])
            return ℓ, D_ℓ, σ_D_ℓ
        else:
            raise ValueError("Invalid version. Choose 1 or 2.")

    def get_polarbear2021(self, remove_neg=False):
        polarbeardat = np.loadtxt('Data/v2/polar.txt')
        ℓ, D_ℓ, σ_D_ℓ = polarbeardat[:, 2], polarbeardat[:, 3], polarbeardat[:, -1]
        if remove_neg:
            sel = np.where(D_ℓ > 0)[0]
            ℓ, D_ℓ, σ_D_ℓ = ℓ[sel], D_ℓ[sel], σ_D_ℓ[sel]
        return ℓ, D_ℓ, σ_D_ℓ

    def get_polarbear2017(self):
        polarbeardat = np.loadtxt('Data/v1/polar.txt')
        ℓ, D_ℓ, σ_D_ℓ = polarbeardat[:, 0], polarbeardat[:, 1], polarbeardat[:, 2]
        return ℓ, D_ℓ, σ_D_ℓ
    
    def get_spt(self):
        sptdat = np.loadtxt('Data/v1/spt.txt')
        ℓ, D_ℓ, σ_D_ℓ = sptdat[:, 0], sptdat[:, -2], sptdat[:, -1]
        return ℓ, D_ℓ, σ_D_ℓ
    def get_sptpol(self):
            """Load SPTpol 500d BB data with bandpowers, covariance, window functions, and beam errors."""
            base_path = 'data/SPTpol/data/sptpol_500d_bb/'
            
            # Bandpowers (3 frequencies: 95x95, 95x150, 150x150)
            _, _, _, bb_95x95, bb_95x150, bb_150x150 = np.loadtxt(
                base_path + 'sptpol_500d_BB_bandpowers.txt', unpack=True, usecols=(0, 1, 2, 3, 4, 5)
            )
            bb_data = np.concatenate((bb_150x150, bb_95x150, bb_95x95))  # Order matches original code
            
            # Covariance matrix
            with open(base_path + 'sptpol_500d_BB_covariance.bin', 'rb') as f:
                cov = np.reshape(np.fromfile(f, dtype=np.float64), (len(bb_data), len(bb_data)))
            
            # Window functions
            with open(base_path + 'sptpol_500d_BB_bpwfs.bin', 'rb') as f:
                l0, l1 = np.fromfile(f, dtype=np.int32, count=2)
                window = np.reshape(np.fromfile(f, dtype=np.float64), (len(bb_data), l1 - l0 + 1))
            ell = np.arange(l0, l1 + 1)
            
            # Beam errors (workaround for 148 elements)
            with open(base_path + 'sptpol_500d_BB_beamerror.bin', 'rb') as f:
                beam_data = np.fromfile(f, dtype=np.float64)
                n_bandpowers = len(bb_data)  # 21
                n_beam_modes = 7  # Assuming 7 modes as intended
                expected_size = n_bandpowers * n_beam_modes  # 147
                if beam_data.size == 148:
                    beam_data = beam_data[:147]  # Trim to 147, assuming extra value is at end
                    print("Warning: Trimmed beam error data from 148 to 147 elements.")
                elif beam_data.size != expected_size:
                    raise ValueError(f"Beam error data size ({beam_data.size}) doesn’t match expected ({expected_size})")
                beam_err = np.reshape(beam_data, (n_bandpowers, n_beam_modes))
            
            return ell, bb_data, cov, window, beam_err

    def get_act2025(self):
        actdat = np.loadtxt('Data/v1/act.txt')
        ℓ, D_ℓ, σ_D_ℓ = actdat[:, 0], actdat[:, 1], actdat[:, 2]
        return ℓ, D_ℓ, σ_D_ℓ