Improvements to the oblate model

[rodluger]

Several improvements can be made to the oblate transit model (starry.Map(oblate=True)):

• Implement the general quartic solver for arbitrary size occultors. The current algorithm can only model cases where there are up to two points of intersection between the limbs of the two bodies. If the star is very oblate or if the occultor is large (such as a companion star), there can be up to four points of intersection, with four separate curves comprising the boundary of the occulted region. This requires a bit of re-structuring of the algorithm on the C++ side, but it's not conceptually difficult.
• Prior to implementing the general case, we should add checks for the condition fproj >= 1 - sqrt(ro) (where fproj is the projected oblateness and ro is the planet-star radius ratio). When this condition is satisfied, occultations with four points of intersection are possible. Raise an error on the theano side when this happens.
• Secondary eclipses are currently not modeled correctly: if the star is oblate and a planet gets occulted by it, the present algorithm ignores the oblateness. It shouldn't be too hard to solve the problem correctly. In the meantime, raise a warning when this happens.