pub fn beam_diffusion_ms(
    sigma_s: Float,
    sigma_a: Float,
    g: Float,
    eta: Float,
    r: Float
) -> Float