def fishqq(lon=None, lat=None, di_block=None):
    """
    Test whether a distribution is Fisherian and make a corresponding Q-Q plot.
    The Q-Q plot shows the data plotted against the value expected from a
    Fisher distribution. The first plot is the uniform plot which is the
    Fisher model distribution in terms of longitude (declination). The second
    plot is the exponential plot which is the Fisher model distribution in terms
    of latitude (inclination). In addition to the plots, the test statistics Mu
    (uniform) and Me (exponential) are calculated and compared against the
    critical test values. If Mu or Me are too large in comparision to the test
    statistics, the hypothesis that the distribution is Fisherian is rejected
    (see Fisher et al., 1987).

    Parameters:
    -----------
    lon : longitude or declination of the data
    lat : latitude or inclination of the data
        or
    di_block: a nested list of [dec,inc]
        A di_block can be provided in which case it will be used instead of
        dec, inc lists.

    Output:
    -----------
    dictionary containing
    lon : mean longitude (or declination)
    lat : mean latitude (or inclination)
    N : number of vectors
    Mu : Mu test statistic value for the data
    Mu_critical : critical value for Mu
    Me : Me test statistic value for the data
    Me_critical : critical value for Me

    if the data has two modes with N >=10 (N and R)
    two of these dictionaries will be returned

    Examples
    --------
    In this example, directions are sampled from a Fisher distribution using
    ``ipmag.fishrot`` and then the ``ipmag.fishqq`` function is used to test
    whether that distribution is Fisherian:

    >>> directions = ipmag.fishrot(k=40, n=50, dec=200, inc=50)
    >>> ipmag.fishqq(di_block = directions)
    {'Dec': 199.73564290371894,
    'Inc': 49.017612342358298,
    'Me': 0.78330310031220352,
    'Me_critical': 1.094,
    'Mode': 'Mode 1',
    'Mu': 0.69915926146177099,
    'Mu_critical': 1.207,
    'N': 50,
    'Test_result': 'consistent with Fisherian model'}

    The above example passed a di_block to the function as an input. Lists of
    paired declination and inclination can also be used as inputs. Here the
    directions di_block is unpacked to separate declination and inclination
    lists using the ``ipmag.unpack_di_block`` functionwhich are then used as
    input to fishqq:

    >>> dec_list, inc_list = ipmag.unpack_di_block(directions)
    >>> ipmag.fishqq(lon=dec_list, lat=inc_list)

    """
    if di_block is None:
        all_dirs = make_di_block(lon, lat)
    else:
        all_dirs = di_block

    ppars = pmag.doprinc(all_dirs)  # get principal directions

    rDIs = []
    nDIs = []
    QQ_dict1 = {}
    QQ_dict2 = {}

    for rec in all_dirs:
        angle = pmag.angle([rec[0], rec[1]], [ppars['dec'], ppars['inc']])
        if angle > 90.:
            rDIs.append(rec)
        else:
            nDIs.append(rec)

    if len(rDIs) >= 10 or len(nDIs) >= 10:
        D1, I1 = [], []
        QQ = {'unf': 1, 'exp': 2}
        if len(nDIs) < 10:
            ppars = pmag.doprinc(rDIs)  # get principal directions
            Drbar, Irbar = ppars['dec'] - 180., -ppars['inc']
            Nr = len(rDIs)
            for di in rDIs:
                d, irot = pmag.dotilt(
                    di[0], di[1], Drbar - 180., 90. - Irbar)  # rotate to mean
                drot = d - 180.
                if drot < 0:
                    drot = drot + 360.
                D1.append(drot)
                I1.append(irot)
                Dtit = 'Mode 2 Declinations'
                Itit = 'Mode 2 Inclinations'
        else:
            ppars = pmag.doprinc(nDIs)  # get principal directions
            Dnbar, Inbar = ppars['dec'], ppars['inc']
            Nn = len(nDIs)
            for di in nDIs:
                d, irot = pmag.dotilt(
                    di[0], di[1], Dnbar - 180., 90. - Inbar)  # rotate to mean
                drot = d - 180.
                if drot < 0:
                    drot = drot + 360.
                D1.append(drot)
                I1.append(irot)
                Dtit = 'Mode 1 Declinations'
                Itit = 'Mode 1 Inclinations'
        plt.figure(figsize=(6, 3))
        Mu_n, Mu_ncr = pmagplotlib.plot_qq_unf(
            QQ['unf'], D1, Dtit, subplot=True)  # make plot
        Me_n, Me_ncr = pmagplotlib.plot_qq_exp(
            QQ['exp'], I1, Itit, subplot=True)  # make plot
        plt.tight_layout()
        if Mu_n <= Mu_ncr and Me_n <= Me_ncr:
            F_n = 'consistent with Fisherian model'
        else:
            F_n = 'Fisherian model rejected'
        QQ_dict1['Mode'] = 'Mode 1'
        QQ_dict1['Dec'] = Dnbar
        QQ_dict1['Inc'] = Inbar
        QQ_dict1['N'] = Nn
        QQ_dict1['Mu'] = Mu_n
        QQ_dict1['Mu_critical'] = Mu_ncr
        QQ_dict1['Me'] = Me_n
        QQ_dict1['Me_critical'] = Me_ncr
        QQ_dict1['Test_result'] = F_n

    if len(rDIs) > 10 and len(nDIs) > 10:
        D2, I2 = [], []
        ppars = pmag.doprinc(rDIs)  # get principal directions
        Drbar, Irbar = ppars['dec'] - 180., -ppars['inc']
        Nr = len(rDIs)
        for di in rDIs:
            d, irot = pmag.dotilt(
                di[0], di[1], Drbar - 180., 90. - Irbar)  # rotate to mean
            drot = d - 180.
            if drot < 0:
                drot = drot + 360.
            D2.append(drot)
            I2.append(irot)
            Dtit = 'Mode 2 Declinations'
            Itit = 'Mode 2 Inclinations'
        plt.figure(figsize=(6, 3))
        Mu_r, Mu_rcr = pmagplotlib.plot_qq_unf(
            QQ['unf'], D2, Dtit, subplot=True)  # make plot
        Me_r, Me_rcr = pmagplotlib.plot_qq_exp(
            QQ['exp'], I2, Itit, subplot=True)  # make plot
        plt.tight_layout()

        if Mu_r <= Mu_rcr and Me_r <= Me_rcr:
            F_r = 'consistent with Fisherian model'
        else:
            F_r = 'Fisherian model rejected'
        QQ_dict2['Mode'] = 'Mode 2'
        QQ_dict2['Dec'] = Drbar
        QQ_dict2['Inc'] = Irbar
        QQ_dict2['N'] = Nr
        QQ_dict2['Mu'] = Mu_r
        QQ_dict2['Mu_critical'] = Mu_rcr
        QQ_dict2['Me'] = Me_r
        QQ_dict2['Me_critical'] = Me_rcr
        QQ_dict2['Test_result'] = F_r

    if QQ_dict2:
        return QQ_dict1, QQ_dict2
    elif QQ_dict1:
        return QQ_dict1
    else:
        print('you need N> 10 for at least one mode')