def plot_consistent_Av(interpolation_kind = "linear"):
    filename = os.path.dirname(os.path.abspath(__file__))\
             + "/Data/" + "taurus.dat"
    taurus = np.genfromtxt(filename, skip_header=1,
                           usecols=(3,7,8,9,10,11,12,1), dtype=None)
    SpT_taurus = [x[0].decode("utf-8") for x in taurus]
    photometry = {"V": [x[1] for x in taurus],
                  "Rc": [x[2] for x in taurus],
                  "Ic": [x[3] for x in taurus],
                  "J":[x[4]+x[5]+x[6] for x in taurus],
                  "H": [x[5]+x[6] for x in taurus],
                  "Ks": [x[6] for x in taurus]}
    star_names = [x[7] for x in taurus]

    #Instantiate a photometric system.
    ps = PhotometricSystem()
    k = list(it.combinations(photometry.keys(), 2))
    #print (sorted(k, key=ps))

    # Create Whitney-Mathis extinction law interpolator:
    Whitney = WhitneyInterpolator(interpolation_kind=interpolation_kind)
    Mathis = MathisInterpolator(interpolation_kind=interpolation_kind)
    WhitneyMathis = MergedFallbackInterpolator([Whitney, Mathis])

    # Create Pecaut & Mamajek pre-MS with extended dwarf color interpolator:
    PM = {}
    C_Av = {}

    #Assumed color in V-Ic is YSO and Extended Dwarf:
    PMY_V_Ic =\
    PecautMamajekYoungInterpolator("SpT", "V-Ic",
                                   interpolation_kind=interpolation_kind)
    PMED_V_Ic =\
    PecautMamajekExtendedDwarfInterpolator("SpT", "V-Ic",
                                           interpolation_kind=\
                                           interpolation_kind)
    PM["V-Ic"] =\
    MergedFallbackInterpolator([PMY_V_Ic,
                                PMED_V_Ic])

    C_Av["V-Ic"] =\
    ColorAvCalculator("V-Ic", PhotometricSystem=ps, ColorInterpolators=[
            PecautMamajekYoungInterpolator,
            PecautMamajekDwarfInterpolator,
            PecautMamajekExtendedDwarfInterpolator
        ],
                      ExtinctionInterpolator=WhitneyMathis,
                      MergingInterpolator=MergedFallbackInterpolator,
                      interpolation_kind="linear")

    Av = {}
    Av_C = {}
    for color in PM:
        Av[color] = calc_Av(PM[color](SpT_taurus),
        np.array(photometry[color.split("-")[0]]) -\
        np.array(photometry[color.split("-")[1]]),
        WhitneyMathis(ps[color.split("-")[0]]) -\
        WhitneyMathis(ps[color.split("-")[1]]))

        Av_C[color] = C_Av[color].evaluate(\
                                np.array(photometry[color.split("-")[0]]) -\
        np.array(photometry[color.split("-")[1]]),
        SpT_taurus)

        plt.plot(Av[color], Av_C[color])
    #plt.plot(0., np.max([np.max(Av[color]) for color in PM]))

    plt.show()