docs: def thermal_fit_result(fit_result, v_residual=None,
                         v_label='Unit-cell volume $(\mathrm{\AA}^3)$',
                         temp_fitline=np.asarray(
                             [300., 1000., 1500., 2000., 2500., 3000.]),
                         figsize=(5, 5), height_ratios=(3, 1), ms_data=50,
                         p_err=None, v_err=None, cbar_loc=(0.99, 0.1, .01, 0.82),
                         pdf_filen=None, title='Fit result'):
      """
      plot P-V-T EOS curve fitting result
  
      :param fit_result: lmfit result object, see example jnb file for detail
      :param v_label: label for volume axis
      :param temp_fitline: temperatures to calculate isothermal compression
          curves, default = [300., 1000., 1500., 2000., 2500., 3000.]
      :param figsize: figure size, default = (7,7)
      :param height_ratios: height ratio between the main and residue plots,
          default = (3,1)
      :param ms_data: marker size for data points
      :param p_err: pressure error bar
      :param v_err: volume error bar
      :param cbar_loc: location of color bar
      :param pdf_filen: name of pdf output file
      :param title: title of the figure
      :return: None
      """
      # basic figure setup
      f, ax = plt.subplots(2, 1, sharex=True, figsize=figsize,
                           gridspec_kw={'height_ratios': height_ratios})
      for ax_i in ax:
          ax_i.tick_params(direction='in')
      # read data to plot
      v_data = fit_result.userkws['v']
      temp_data = fit_result.userkws['temp']
      p_data = fit_result.data
      p_datafit = fit_result.best_fit
      v0 = uct.ufloat(fit_result.params['st_v0'].value,
                      fit_result.params['st_v0'].stderr)
      sm = plt.cm.ScalarMappable(cmap=c_map,
                                 norm=plt.Normalize(
                                     vmin=300., vmax=temp_data.max()))
      a = sm.to_rgba(temp_fitline)
      v_fitline = np.linspace(v0.n, min(v_data), 1000)
      fitmodel_copy = copy.deepcopy(fit_result)
      for a_i, temp_i in zip(a, temp_fitline):
          p_fitline = fitmodel_copy.eval(v=v_fitline,
                                         temp=np.ones_like(v_fitline) * temp_i)
          ax[0].plot(p_fitline, v_fitline, c=a_i)
      # error range here does not make a lot sense, so not supported
      # if (p_err is not None) and (v_err is not None):
      ax[0].errorbar(p_data, v_data, xerr=p_err, yerr=v_err, fmt=' ', c='k',
                     capsize=0, elinewidth=0.5, label='Data', zorder=0)
      points = ax[0].scatter(p_data, v_data, marker='o', s=ms_data, c=temp_data,
                             cmap=c_map, vmin=300., vmax=temp_data.max(),
                             zorder=1)
      if v_residual is None:
          ax[1].scatter(p_data, p_data - p_datafit, marker='o', s=ms_data,
                        c=temp_data, cmap=c_map, vmin=300.,
                        vmax=temp_data.max(), zorder=1)
          ax[1].errorbar(p_data, p_data - p_datafit, yerr=p_err, fmt=' ', c='k',
                         capsize=0, elinewidth=0.5, label='Data', zorder=0)
          ax[1].set_ylabel('$P_{obs} - P_{fit}$')
      else:
          ax[1].scatter(p_data, v_residual, marker='o', s=ms_data, c=temp_data,
                        cmap=c_map, vmin=300., vmax=temp_data.max(), zorder=1)
          ax[1].errorbar(p_data, v_residual, yerr=p_err, fmt=' ', c='k',
                         capsize=0, elinewidth=0.5, label='Data', zorder=0)
          ax[1].set_ylabel('$V_{obs} - V_{fit}$')
      # ax[0].legend()
      position = f.add_axes(cbar_loc)
      f.colorbar(points, orientation="vertical", cax=position,
                 ticks=temp_fitline)
      ax[1].axhline(0, c='k', ls='--')
      ax[1].set_xlabel('Pressure (GPa)')
      ax[0].set_ylabel(v_label)
      ax[0].set_title(title)
      plt.tight_layout()
      if pdf_filen is not None:
          f.savefig(pdf_filen)