docs: def plot_hdd(HDD, B, M, s):
      """
      Function to make hysteresis, deltaM and DdeltaM plots
      Parameters:
      _______________
      Input
          HDD :  dictionary with figure numbers for the keys:
              'hyst' : hysteresis plot  normalized to maximum value
              'deltaM' : Delta M plot
              'DdeltaM' : differential of Delta M plot
          B : list of field values in tesla
          M : list of magnetizations in arbitrary units
          s : specimen name string
      Ouput
        hpars : dictionary of hysteresis parameters with keys:
          'hysteresis_xhf', 'hysteresis_ms_moment', 'hysteresis_mr_moment', 'hysteresis_bc'
  
      """
      hpars, deltaM, Bdm = plot_hys(
          HDD['hyst'], B, M, s)  # Moff is the "fixed" loop data
      DdeltaM = []
      Mhalf = ""
      for k in range(2, len(Bdm)):
          # differnential
          DdeltaM.append(
              old_div(abs(deltaM[k] - deltaM[k - 2]), (Bdm[k] - Bdm[k - 2])))
      for k in range(len(deltaM)):
          if old_div(deltaM[k], deltaM[0]) < 0.5:
              Mhalf = k
              break
      try:
          Bhf = Bdm[Mhalf - 1:Mhalf + 1]
          Mhf = deltaM[Mhalf - 1:Mhalf + 1]
          # best fit line through two bounding points
          poly = np.polyfit(Bhf, Mhf, 1)
          Bcr = old_div((.5 * deltaM[0] - poly[1]), poly[0])
          hpars['hysteresis_bcr'] = '%8.3e' % (Bcr)
          hpars['magic_method_codes'] = "LP-BCR-HDM"
          if HDD['deltaM'] != 0:
              plot_delta_m(HDD['deltaM'], Bdm, deltaM, Bcr, s)
              plt.axhline(0, color='k')
              plt.axvline(0, color='k')
              plot_d_delta_m(HDD['DdeltaM'], Bdm, DdeltaM, s)
      except:
          hpars['hysteresis_bcr'] = '0'
          hpars['magic_method_codes'] = ""
      return hpars