"""Low-frequency (M omega << 1) amplification factor of a massless scalar on Kerr.

Brito, Cardoso and Pani, Lect. Notes Phys. 906 (2015), arXiv:1501.06570v8, Eq. (4.104) with s = 0,
in units G = c = M = 1:
    Z = -8 r_+ (omega - m Omega_H) omega^(2l+1) (r_+ - r_-)^(2l) [(l!)^2 / ((2l)! (2l+1)!!)]^2
        * prod_{k=1}^{l} [1 + ((omega - m Omega_H) / (k pi r_+ T_H))^2 ],   T_H = (r_+ - r_-) / (4 pi r_+^2).
The sign is set by -(omega - m Omega_H): the same superradiance condition as the exact solution.
"""
import math


def double_factorial(n):
    return math.prod(range(n, 0, -2)) if n > 0 else 1


def low_frequency_amplification(omega, a, l=1, m=1):
    root = math.sqrt(1 - a * a)
    r_plus, r_minus = 1 + root, 1 - root
    omega_h = a / (2 * r_plus)
    temperature = (r_plus - r_minus) / (4 * math.pi * r_plus ** 2)
    k_h = omega - m * omega_h
    prefactor = (math.factorial(l) ** 2 / (math.factorial(2 * l) * double_factorial(2 * l + 1))) ** 2
    product = 1.0
    for k in range(1, l + 1):
        product *= 1 + (k_h / (k * math.pi * r_plus * temperature)) ** 2
    return -8 * r_plus * k_h * omega ** (2 * l + 1) * (r_plus - r_minus) ** (2 * l) * prefactor * product