def P_trans(N, M, DeltaJ, DeltaM):
    """
    Calculates transition probability.

    Parameters
    ----------
    N : int
        Total rotational angular momentum quantum number
    M : int
        Magnetic quantum number
    DeltaJ : int
        Change in rotational quantum number
    DeltaM : int
        Change in magnetic quantum number

    Returns
    -------
    float
        Transition probability
    """
    # fmt: off
    if DeltaJ == 1: # Liebe 1981
        if DeltaM == 1:
            return 3 * (N + M + 1) * (N + M + 2) / (4 * (N + 1) * (2*N + 1) * (2*N + 3))
        elif DeltaM == 0:
            return 3 * ((N + 1)**2 - M**2) / ((N + 1) * (2 * N + 1) * (2 * N + 3))
        elif DeltaM == -1:
            return 3 * (N - M + 1) * (N - M + 2) / (4 * (N + 1) * (2*N + 1) * (2*N + 3))
    elif DeltaJ == -1:
        if DeltaM == 1:
            return 3 * (N + 1) * (N - M) * (N - M - 1) / (4 * N * (2 * N + 1) * (2 * N**2 + N - 1))
        elif DeltaM == 0:
            return 3 * (N + 1) * (N**2 - M**2) / (N * (2 * N + 1) * (2 * N**2 + N - 1))
        elif DeltaM == -1:
            return 3 * (N + 1) * (N + M) * (N + M - 1) / (4 * N * (2 * N + 1) * (2 * N**2 + N - 1))
    # fmt: on