Function br_heating_integral

pub fn br_heating_integral(
    x_grid: &[f64],
    delta_n: &[f64],
    theta_z: f64,
    theta_e: f64,
    n_h: f64,
    n_he: f64,
    n_e: f64,
    x_e_frac: f64,
    y_he_ii: f64,
    y_he_i: f64,
) -> f64

BR emission rate integrated over frequency (for electron temperature equation).

H_BR = (1/(4 G₃ θ_z)) ∫ [1 − n(e^{x_e} − 1)] K_BR(x)/x³ · x³ dx = (1/(4 G₃ θ_z)) ∫ [1 − n(e^{x_e} − 1)] K_BR(x) dx

This enters the electron temperature evolution as a cooling term, analogous to dc_heating_integral() in double_compton.rs.