def fit_linear(mu, I):
    """
    Calculate the coefficients for the linear LD law.
    It assumes input intensities are normalized.  For a derivation of the
    least-squares problem solved, see Espinoza & Jordan (2015).

    INPUTS:
      mu:  Angles at which each intensity is calculated (numpy array).
      I:   Normalized intensities (i.e., I(mu)/I(1)) (numpy array).

    OUTPUTS:
      a:   Coefficient of the linear law.
    """
    alpha_1_1 = sum((1.0-mu)**2)
    beta_1 = sum((1.0-mu)*(1.0-I))
    a = beta_1/alpha_1_1
    return a