Update `position` for how solar system masses will deflect its light.

    Given the ICRS `position` [x,y,z] of an object (au) that is being
    viewed from the `observer` also expressed as [x,y,z], and given an
    ephemeris that can be used to determine solar system body positions,
    and given the time `t` and Boolean `apply_earth` indicating whether
    to worry about the effect of Earth's mass, and a `count` of how many
    major solar system bodies to worry about, this function updates
    `position` in-place to show how the masses in the solar system will
    deflect its image.

def add_deflection(position, observer, ephemeris, t,
                   include_earth_deflection, count=3):
    """Update `position` for how solar system masses will deflect its light.

    Given the ICRS `position` [x,y,z] of an object (au) that is being
    viewed from the `observer` also expressed as [x,y,z], and given an
    ephemeris that can be used to determine solar system body positions,
    and given the time `t` and Boolean `apply_earth` indicating whether
    to worry about the effect of Earth's mass, and a `count` of how many
    major solar system bodies to worry about, this function updates
    `position` in-place to show how the masses in the solar system will
    deflect its image.

    """
    # Compute light-time to observed object.

    tlt = length_of(position) / C_AUDAY

    # Cycle through gravitating bodies.

    jd_tdb = t.tdb
    ts = t.ts
    for name in deflectors[:count]:
        try:
            deflector = ephemeris[name]
        except KeyError:
            deflector = ephemeris[name + ' barycenter']

        # Get position of gravitating body wrt ss barycenter at time 't_tdb'.

        bposition = deflector.at(ts.tdb(jd=jd_tdb)).position.au  # TODO

        # Get position of gravitating body wrt observer at time 'jd_tdb'.

        gpv = bposition - observer

        # Compute light-time from point on incoming light ray that is closest
        # to gravitating body.

        dlt = light_time_difference(position, gpv)

        # Get position of gravitating body wrt ss barycenter at time when
        # incoming photons were closest to it.

        tclose = jd_tdb

        # if dlt > 0.0:
        #     tclose = jd - dlt

        tclose = where(dlt > 0.0, jd_tdb - dlt, tclose)
        tclose = where(tlt < dlt, jd_tdb - tlt, tclose)

        # if tlt < dlt:
        #     tclose = jd - tlt

        bposition = deflector.at(ts.tdb(jd=tclose)).position.au  # TODO
        rmass = rmasses[name]
        _add_deflection(position, observer, bposition, rmass)

    # If observer is not at geocenter, add in deflection due to Earth.

    if include_earth_deflection.any():
        deflector = ephemeris['earth']
        bposition = deflector.at(ts.tdb(jd=tclose)).position.au  # TODO
        rmass = rmasses['earth']
        # TODO: Make the following code less messy, maybe by having
        # _add_deflection() return a new vector instead of modifying the
        # old one in-place.
        deflected_position = position.copy()
        _add_deflection(deflected_position, observer, bposition, rmass)
        if include_earth_deflection.shape:
            position[:,include_earth_deflection] = (
                deflected_position[:,include_earth_deflection])
        else:
            position[:] = deflected_position[:]