import numpy as np

def campbell(ala):
    """
    Campbell's leaf angle distribution function (ellipsoidal).
    
    Parameters:
        ala : average leaf angle in degrees

    Returns:
        freq : relative frequency for 13 inclination bins
        litab: bin centers (degrees)
    """
    tx1 = np.array([10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88., 90.])
    tx2 = np.array([0., 10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88.])
    litab = (tx1 + tx2) / 2
    n = len(litab)

    tl1 = np.radians(tx1)
    tl2 = np.radians(tx2)

    excent = np.exp(-1.6184e-5 * ala**3 + 2.1145e-3 * ala**2 - 0.12390 * ala + 3.2491)

    freq = np.zeros(n)
    for i in range(n):
        x1 = excent / np.sqrt(1 + (excent**2) * np.tan(tl1[i])**2)
        x2 = excent / np.sqrt(1 + (excent**2) * np.tan(tl2[i])**2)

        if excent == 1:
            freq[i] = abs(np.cos(tl1[i]) - np.cos(tl2[i]))
        else:
            alpha = excent / np.sqrt(abs(1 - excent**2))
            alpha2 = alpha**2
            x1_sq = x1**2
            x2_sq = x2**2

            if excent > 1:
                alpx1 = np.sqrt(alpha2 + x1_sq)
                alpx2 = np.sqrt(alpha2 + x2_sq)
                dum1 = x1 * alpx1 + alpha2 * np.log(x1 + alpx1)
                dum2 = x2 * alpx2 + alpha2 * np.log(x2 + alpx2)
                freq[i] = abs(dum1 - dum2)
            else:
                almx1 = np.sqrt(alpha2 - x1_sq)
                almx2 = np.sqrt(alpha2 - x2_sq)
                dum1 = x1 * almx1 + alpha2 * np.arcsin(x1 / alpha)
                dum2 = x2 * almx2 + alpha2 * np.arcsin(x2 / alpha)
                freq[i] = abs(dum1 - dum2)

    freq_sum = np.sum(freq)
    freq_normalized = freq / freq_sum if freq_sum != 0 else freq

    return freq_normalized, litab