def I(r, limbdark):
  '''
  The standard quadratic limb darkening law.
  
  :param ndarray r: The radius vector
  :param limbdark: A :py:class:`pysyzygy.transit.LIMBDARK` instance containing 
                   the limb darkening law information
  
  :returns: The stellar intensity as a function of `r`
  
  '''
  
  if limbdark.ldmodel == QUADRATIC:
    u1 = limbdark.u1
    u2 = limbdark.u2
    return (1-u1*(1-np.sqrt(1-r**2))-u2*(1-np.sqrt(1-r**2))**2)/(1-u1/3-u2/6)/np.pi
  elif limbdark.ldmodel == KIPPING:
    a = np.sqrt(limbdark.q1)
    b = 2*limbdark.q2
    u1 = a*b
    u2 = a*(1 - b)
    return (1-u1*(1-np.sqrt(1-r**2))-u2*(1-np.sqrt(1-r**2))**2)/(1-u1/3-u2/6)/np.pi
  elif limbdark.ldmodel == NONLINEAR:
    raise Exception('Nonlinear model not yet implemented!')                           # TODO!
  else:
    raise Exception('Invalid limb darkening model.')
The standard quadratic limb darkening law.
  
  :param ndarray r: The radius vector
  :param limbdark: A :py:class:`pysyzygy.transit.LIMBDARK` instance containing 
                   the limb darkening law information
  
  :returns: The stellar intensity as a function of `r`