def top_down_SolarSystem(discoveries,
                         inner_limit=6,      # truncate at 8 AU to show that we don't have sensitivity in close
                         extent=83,          # extend to 83 AU as indicative only, sensitivity is to ~300 AU
                         plot_discoveries=None,
                         plot_colossos=False,
                         plot_blocks=None,
                         plot_galaxy=False,
                         feature_blocks=None,
                         plot_Ijiraq=False,  # detected Ijiraq at 9.80 AU in the 13AE block
                         label_blocks=True,
                         savefilename=None):
    """
    Plot the OSSOS discoveries on a top-down Solar System showing the position of Neptune and model TNOs.
    Discoveries are plotted each at their time of discovery according to the value in the Version Release.
    Coordinates should be polar to account for RA hours, radial axis is in AU.
    :return: a saved plot

    This is not a 'projection' of the particles into any common plane.  Each wedge is a different latitude above
    the plane, but inside each wedge it is heliocentric distance vs RA.
     That is fine.  It's just important to know that if someone asks we know this is NOT (for example) a projection of each
    object down into the ecliptic.
     A very minor future long-term improvement is that the galactic
    plane wedge is not (I don't think) symmetric around 18h and 6h.
    I think that lines of constant (say 15 degree galactic lat) are
    centered a little off (I think the wedge would 'twist' a little
    bit counter-clockwise).  I am not absolutely sure along the
    ecliptic where the b= +/- 15 degree lines are but I don't
    think they are symmetric?
    """
    fig = plt.figure(figsize=(6, 6))
    rect = [0.01, 0.01, 0.95, .95]  # the plot occupies not all the figspace

    ax1 = fig.add_axes(rect, polar=True, frameon=False)  # theta (RA) is zero at E, increases anticlockwise
    ax1.set_aspect('equal')

    ax1.set_rlim(0, extent)
    ax1.set_rgrids([20, 40, 60, 80], labels=["", "", '20 au', '40 au', '60 au', '80 au'], angle=190, alpha=0.45)  # angle = 308
    ax1.yaxis.set_major_locator(MultipleLocator(20))
    ax1.xaxis.set_major_locator(MultipleLocator(math.radians(15)))  # every 2 hours
    ax1.grid(axis='x', color='k', linestyle='--', alpha=0.2)
    ax1.set_xticklabels(['', '0h', "", '2h', "", '4h', "", '6h', "", '8h', "", '10h', "", '12h', "", '14h', "", '16h',
                         "", '18h', "", '20h', "", '22h', "", ],
                        # ""])  # str(r)+'h' for r in range(-1,24)],
                        #        ['', '0h', '2h', '4h', '6h', '8h', '10h', '12h', '14h', '16h', '18h', '20h', '22h'],
                        color='b', alpha=0.6)  # otherwise they get in the way


    if plot_galaxy:
        # plot exclusion zones due to Galactic plane: RAs indicate where bar starts, rather than its centre angle
        width = math.radians(3 * 15)
        plt.bar(math.radians(4.5 * 15), extent, width=width, color=plot_fanciness.ALMOST_BLACK, linewidth=0, alpha=0.2)
        plt.bar(math.radians(16.5 * 15), extent, width=width, color=plot_fanciness.ALMOST_BLACK, linewidth=0, alpha=0.2)
        ax1.annotate('galactic plane', (math.radians(6.9 * 15), extent - 15), size=10, color='k', alpha=0.45)
        ax1.annotate('  galactic plane\navoidance zone', (math.radians(16.9 * 15), extent - 12), size=10, color='k', alpha=0.45)

    # again RAs indicate where bar starts, so subtract half the block width from the block centrepoints
    # and approximate blocks as math.radians(7) degrees wide.

    for blockname, block in parameters.BLOCKS.items():
        if plot_blocks:
            if blockname in plot_blocks:
                if feature_blocks is not None and blockname in feature_blocks:
                    colour = 'm' #'#E47833'
                    alpha = 0.25
                else:
                    colour = 'b'
                    alpha = 0.1
                # if blockname.startswith('13') or blockname.startswith('14'):
                    width = math.radians(7) # the 3 x 7 field layout
                # else:
                #     width = math.radians(5) # the 4 x 5 field layout

                # cmap = plt.get_cmap('Blues')
                # if blockname.endswith('AO'):
                #     colour = '#E47833'
                #     alpha = 0.17
                #     cmap = plt.get_cmap('Oranges')  # colorbrewer.sequential.Oranges_3.mpl_colors

                ax1.bar(ephem.hours(block["RA"]) - math.radians(3.5), extent, linewidth=0.1,
                        width=width, bottom=inner_limit, zorder=0, color=colour, alpha=alpha)
                if label_blocks:
                    ax1.annotate(blockname[3], (ephem.hours(block["RA"]) + math.radians(0.36), extent - 3.), size=15,
                                 color='b')

    # No optional on these just yet
    plot_planets_plus_Pluto(ax1)
    ra, dist, hlat, Hmag = parsers.synthetic_model_kbos(kbotype='resonant', arrays=True, maglimit=24.7)
    # can't plot Hmag as marker size in current setup.
    ax1.scatter(ra, dist, marker='o', s=2, facecolor=plot_fanciness.ALMOST_BLACK,
                edgecolor=plot_fanciness.ALMOST_BLACK, linewidth=0.1, alpha=0.12, zorder=1)

    if plot_discoveries is not None:
        plot_ossos_discoveries(ax1, discoveries, plot_discoveries, plot_colossos=plot_colossos)

    # special detection in 13AE: Saturn's moon Ijiraq at 2013-04-09 shows inner limit of sensitivity.
    if plot_Ijiraq:
        # Position from Horizons as it's excluded from the list of detections
        ax1.scatter(ephem.hours('14 29 46.57'), 9.805,
                    marker='o', s=4, facecolor='b', edgecolor=plot_fanciness.ALMOST_BLACK, linewidth=0.15, alpha=0.8)

    plt.draw()
    if savefilename is not None:
        outfile = '{}.pdf'.format(savefilename)
    else:
        outfile = 'topdown_RA_d_OSSOS_v{}.pdf'.format(parameters.RELEASE_VERSION)
    plt.savefig(outfile, transparent=True, bbox_inches='tight')

    return