docs: 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