def bns_postmerger_frequency(label, kappa , mass, nu):
    """
        Calculate the merger and postmerger frequencies with the aligned-spin NR fit
        by Matteo Breschi, Sebastiano Bernuzzi, Francesco Zappa et al.
        [https://arxiv.org/abs/1908.11418]
        this versions is calibrated with BAM and THC data
        
        label : 1,2,3 or 'm'
        kappa : tidal polarizability
        mass  : total mass [solar masses]
        nu    : symmetric mass ratios
    """
    if label == 1:
        Q  = 0.052182
        n1 = 0.002843
        n2 = 0.0
        d1 = 0.012868
        d2 = 0.0
        c  = 5767.6
    elif label == 2:
        Q  = 7.6356
        n1 = 0.066645
        n2 = 0.000040146
        d1 = 10.949
        d2 = 0.040276
        c  = -52.655
    elif label == 3:
        Q  = 4.5722
        n1 = 0.060385
        n2 = 0.00010661
        d1 = 4.1506
        d2 = 0.027552
        c  = 1875.5
    elif label == 'm':
        Q  = 0.033184
        n1 = 0.0013067
        n2 = 0.0
        d1 = 0.0050064
        d2 = 0.0
        c  = 3199.75613
    else:
        raise ValueError("Unknown label for bns_postmerger_frequency.\n Please use 1, 2, 3 or 'm' (merger) respectively for f_{2-0} , f_2 , f_{2+0} and f_merg")
        exit()
    MTsun = 4.925491025543576e-06
    
    # f_2-0 fit is not realible for kappa > 500
    x = kappa + c * (1 - 4*nu)
    if kappa > 500 and label == 1:
        x0  = 500 + c * (1 - 4*nu)
        q0  = Q*( (1 + n1*x0 + n2*x0*x0) / (1 + d1*x0 + d2*x0*x0) )
        m0  = Q*( (n1 - d1 + 2*(n2-d2)*x0 + (d1*n2 - d2*n1)*x0*x0) / ((1 + d1*x0 + d2*x0*x0)*(1 + d1*x0 + d2*x0*x0)))
        f = q0 + m0 * (x - x0)
    else:
        f = Q * (1 + n1*x + n2*x*x) / (1 + d1*x + d2*x*x)

    return f/(MTsun*mass)