def iplot_hys(fignum, B, M, s):
    """
    function to plot hysteresis data

    This function has been adapted from pmagplotlib.iplot_hys for specific use
    within a Jupyter notebook.

    Parameters
    -----------
    fignum : reference number for matplotlib figure being created
    B : list of B (flux density) values of hysteresis experiment
    M : list of M (magnetization) values of hysteresis experiment
    s : specimen name
    """
    if fignum != 0:
        plt.figure(num=fignum)
        plt.clf()
    hpars = {}
# close up loop
    Npts = len(M)
    B70 = 0.7 * B[0]  # 70 percent of maximum field
    for b in B:
        if b < B70:
            break
    Nint = B.index(b) - 1
    if Nint > 30:
        Nint = 30
    if Nint < 10:
        Nint = 10
    Bzero, Mzero, Mfix, Mnorm, Madj, MadjN = "", "", [], [], [], []
    Mazero = ""
    m_init = 0.5 * (M[0] + M[1])
    m_fin = 0.5 * (M[-1] + M[-2])
    diff = m_fin - m_init
    Bmin = 0.
    for k in range(Npts):
        frac = old_div(float(k), float(Npts - 1))
        Mfix.append((M[k] - diff * frac))
        if Bzero == "" and B[k] < 0:
            Bzero = k
        if B[k] < Bmin:
            Bmin = B[k]
            kmin = k
    # adjust slope with first 30 data points (throwing out first 3)
    Bslop = B[2:Nint + 2]
    Mslop = Mfix[2:Nint + 2]
    polyU = polyfit(Bslop, Mslop, 1)  # best fit line to high field points
    # adjust slope with first 30 points of ascending branch
    Bslop = B[kmin:kmin + (Nint + 1)]
    Mslop = Mfix[kmin:kmin + (Nint + 1)]
    polyL = polyfit(Bslop, Mslop, 1)  # best fit line to high field points
    xhf = 0.5 * (polyU[0] + polyL[0])  # mean of two slopes
    # convert B to A/m, high field slope in m^3
    hpars['hysteresis_xhf'] = '%8.2e' % (xhf * 4 * np.pi * 1e-7)
    meanint = 0.5 * (polyU[1] + polyL[1])  # mean of two intercepts
    Msat = 0.5 * (polyU[1] - polyL[1])  # mean of saturation remanence
    Moff = []
    for k in range(Npts):
        # take out linear slope and offset (makes symmetric about origin)
        Moff.append((Mfix[k] - xhf * B[k] - meanint))
        if Mzero == "" and Moff[k] < 0:
            Mzero = k
        if Mzero != "" and Mazero == "" and Moff[k] > 0:
            Mazero = k
    hpars['hysteresis_ms_moment'] = '%8.3e' % (Msat)  # Ms in Am^2
#
# split into upper and lower loops for splining
    Mupper, Bupper, Mlower, Blower = [], [], [], []
    deltaM, Bdm = [], []  # diff between upper and lower curves at Bdm
    for k in range(kmin - 2, 0, -2):
        Mupper.append(old_div(Moff[k], Msat))
        Bupper.append(B[k])
    for k in range(kmin + 2, len(B)-1):
        Mlower.append(Moff[k] / Msat)
        Blower.append(B[k])
    Iupper = spline.Spline(Bupper, Mupper)  # get splines for upper up and down
    Ilower = spline.Spline(Blower, Mlower)  # get splines for lower
    for b in np.arange(B[0]):  # get range of field values
        Mpos = ((Iupper(b) - Ilower(b)))  # evaluate on both sides of B
        Mneg = ((Iupper(-b) - Ilower(-b)))
        Bdm.append(b)
        deltaM.append(0.5 * (Mpos + Mneg))  # take average delta M
    print('whew')
    for k in range(Npts):
        MadjN.append(old_div(Moff[k], Msat))
        Mnorm.append(old_div(M[k], Msat))
# find Mr : average of two spline fits evaluted at B=0 (times Msat)
    Mr = Msat * 0.5 * (Iupper(0.) - Ilower(0.))
    hpars['hysteresis_mr_moment'] = '%8.3e' % (Mr)
# find Bc (x intercept), interpolate between two bounding points
    Bz = B[Mzero - 1:Mzero + 1]
    Mz = Moff[Mzero - 1:Mzero + 1]
    Baz = B[Mazero - 1:Mazero + 1]
    Maz = Moff[Mazero - 1:Mazero + 1]
    try:
        poly = polyfit(Bz, Mz, 1)  # best fit line through two bounding points
        Bc = old_div(-poly[1], poly[0])  # x intercept
        # best fit line through two bounding points
        poly = polyfit(Baz, Maz, 1)
        Bac = old_div(-poly[1], poly[0])  # x intercept
        hpars['hysteresis_bc'] = '%8.3e' % (0.5 * (abs(Bc) + abs(Bac)))
    except:
        hpars['hysteresis_bc'] = '0'
    return hpars, deltaM, Bdm, B, Mnorm, MadjN