docs: def _kip_cont2(self, sparse, cycle_start=0, cycle_end=0,
                    plot=['dcoeff'], thresholds=[1.0E+12],
                    xax='log_time_left', alphas=[1.0], yllim=0., yulim=0.,
                    y_res=2000, xllim=0., xulim=0., age='seconds',
                    sparse_intrinsic=20, engen=False,
                    netnuc_name='eps_nuc', engenalpha=0.6,
                    outfile='plot.pdf', annotation='',
                    KEPLER=False):
          """
          !! EXPERIMENTAL FEATURE (flagged as private) !!
  
          This function creates a Kippenhahn diagram as a contour plot of
          the .se.h5 or .out.h5 files using any continuous variable
          (columns in the hdf5 cycle data). Multiple columns may be
          plotted, their name indicated in the list "plot", and their
          thresholds in the list "thresholds".
  
          Currently, this is only designed to take one threshold for each
          variable but future versions will hopefully be able to plot with
          multiple thresholds, however you may circumvent this issue by
          repeating the variable in "plots" and entering a second
          threshold for it in "thresholds".
  
          Parameters
          ----------
          sparse : integer
              x-axis (timestep) sparsity.  The true sparsity is
              sparse*sparse_intrinsic.  Try 100 or 500 for .se.h5 and 20
              for .out.h5 files and for preliminaxllimry plots.
          cycle_start : integer, optional
              Cycle from which you wish to plot.  The default is 0.
          cycle_end : integer, optional
              Maximum cycle that you wish to plot.  If cycle_end is 0,
              then it will plot up to the last cycle available.  The
              default is 0.
          plot : list, optional
              1-D array containing the variables to be plotted (as
              strings, e.g. plots=['dcoeff','C-13'].  I recommend always
              plotting 'dcoeff' as plots[0]).  The default is ['dcoeff'].
          thresholds : list, optional
              1-D array containing the thresholds corresponding to the
              variables in "plots".  The default is [1.0E+12].
          xax : string, optional
              x-axis quantity; either 'log_time_left' or 'cycles'.  The
              default is 'log_time_left'.
          alphas : list, optional
              Array containing the opacity (0 to 1) of the contour for
              each variable.  The default is [1.0].
          yllim : float, optional
              Lower plot limit for y-axis (mass co-ordinate).  The default
              is 0..
          yulim : float, optional
              Upper plot limit for y-axis (mass co-ordinate).  The default
              is 0..
          y_res : integer, optional
              y-axis resolution. Defaults to 2000 but increasing to as
              much as 10000 does not significantly affect the plotting
              time.  The default is 2000.
          xllim : float, optional
              Lower plot limit for x-axis.  The default is 0..
          xulim : float, optional
              Upper plot limit for x-axis.  The default is 0..
          age : string, optional
              Either 'years' or 'seconds', depending on the data.  The
              default is 'seconds'.
          sparse_intrinsic : integer, optional
              Sparsity of timesteps in the data provided (usually 20 for
              .out.h5 files and 1 for .se.h5 files).  The default is 20.
          engen : boolean, optional
              Whether the user would like to plot Kippenhahn of convective
              zones and energy generation.  If True, please still include
              plots=['dcoeff'] and thresholds=[1.0E+12'] in your call.
              This will require the data to have an 'eps_nuc' column, so
              the plot is only working for .se.h5 files from MESA in the
              current se library.  This is the most recent addition, so
              probably the most buggy.  The plot script will automatically
              calculate and assign multiple thresholds according to the
              model.  The default is False.
          netnuc_name : string, optional
              The name of the column containing (eps_nuc-eps_neu).  If you
              do not have available eps_neu, then you can give
              netnuc_name="eps_nuc" to just plot energy generation.  The
              default is "eps_nuc".
          engenalpha : float, optional
              Opacity of the energy generation contours.  The default
              is 0.6.
          outfile : string, optional
              The name to save the plot as.  The default is 'plot.pdf'.
          annotation : string, optional
              Some optional annotation to add to the plot.   The default
              is ''.
          KEPLER : boolean, optional
              The default is False.
  
          """
  
          # Organize cycles and ages:
          original_cyclelist = self.se.cycles
          if cycle_end==0.:
              cycle_end = original_cyclelist[-1]
          cycle_end = old_div(int(cycle_end),sparse_intrinsic) - 1
          if cycle_start==0:
              pass
          else:
              cycle_start = old_div(int(cycle_start),sparse_intrinsic) - 1
          cyclelist = original_cyclelist[cycle_start:cycle_end:sparse]
          # fix for KEPLER restart counting at O burning:
          original_ages= self.se.ages
          age_at_restart = 0.
          age_at_restart_idx = 0
          if KEPLER == True:
              for i in range(1,len(original_ages)):
                  if (original_ages[i]-original_ages[i-1]) < 0.:
                      age_at_restart_idx = i-1
                      age_at_restart = original_ages[i-1]
                      print('age restart found at cycle = '+str(age_at_restart_idx)+', age = '+str(age_at_restart))
                      KEPLER = True
                      break
  
              for i in range(age_at_restart_idx+1,len(original_ages)):
                  original_ages[i] = original_ages[i] + age_at_restart
  
  
          # Figure:
          fig = pl.figure()
          ax = pl.axes()
          params = {'axes.labelsize':  15,
            'text.fontsize':   15,
            'legend.fontsize': 15,
            'xtick.labelsize': 15,
            'ytick.labelsize': 15,
            'text.usetex': False}
          fsize=18
          pl.rcParams.update(params)
          # X-axis:
          if xax == 'log_time_left':
              if KEPLER == True:
                  xxtmp = original_ages[cycle_start:cycle_end:sparse]
              else:
                  xxtmp = self.se.ages[cycle_start:cycle_end:sparse]
              if age == 'years':
                  if KEPLER == True:
                      pass
                  else:
                      xxtmp = self.se.ages[cycle_start:cycle_end:sparse]
              elif age == 'seconds':
                  if KEPLER == True:
                      for i in range(len(cyclelist)):
                          xxtmp[i] = old_div(original_ages[cycle_start:cycle_end:sparse][i],31558149.984)
                  else:
                      for i in range(len(cyclelist)):
                          xxtmp[i] = old_div(self.se.ages[cycle_start:cycle_end:sparse][i],31558149.984)
          if xax == 'cycles':
              xx = cyclelist
              xxtmp = cyclelist
          # Set up x-axis according to whether ages are in years or seconds and
          # re-write as log(time left). The last entry will always be -inf in this
          # way so we calculate it by extrapolating the anti-penultimate and
          # penultimate entries.
          # Modified from the GENEC gdic (R. Hirschi)
          if xax == 'log_time_left':
              xx=np.zeros(len(xxtmp))
              agemin = max(old_div(abs(xxtmp[-1]-xxtmp[-2]),5.),1.e-10)
              for i in np.arange(len(xxtmp)):
                  if xxtmp[-1]-xxtmp[i]>agemin:
                      xx[i]=np.log10(xxtmp[-1]-xxtmp[i]+agemin)
                  else :
                      xx[i]=np.log10(agemin)
              ax.set_xlabel('$log_\mathrm{10}(t_\mathrm{end}-t)\;[\mathrm{yr}]$',fontsize=fsize-1)
          if xax == 'cycles':
              ax.set_xlabel('$\mathrm{CYCLE}$',fontsize=fsize-1)
          # Y-axis limits and resolution:
          totalmass = []
          try:
              m_ini = float(self.se.get('mini'))
          except:
              mini=m.se.get(0,'total_mass')
              mini=old_div(mini,constants.mass_sun)
              print('getting Mini from 1st cycle')
  
          if yulim==0.:
              yulim = m_ini
          dy = old_div(m_ini,float(y_res))
          vlinetol = 1.0E-8
          # Set up (y-axis) vector and a 3-D (hist) array to store all of the
          # contours.
          y = np.arange(0., m_ini, dy)
          if engen == True:
              Z = np.zeros([len(y),len(xxtmp),len(plot)+2],float)
          else:
              Z = np.zeros([len(y),len(xxtmp),len(plot)],float)
  
          # Define function extracting the contour boundaries which will be
          # called for every cycle in cyclelist, for every variable to be plotted
          # along with its corresponding threshold(s).
          def getlims(variable_array,thresh,massco_array):
              """This function returns the variable boundaries (in mass) for a cycle,
              given the cycle's variable and mass columns, ensuring that the boundaries
              are ordered centre to surface (as some .se.h5 files are the opposite)."""
              plotlims = []
              # Just a fix for some get problem I was having:
              if len(massco_array) == 2:
                  massco_array = massco_array[0]
              if len(variable_array) == 2:
                  variable_array = variable_array[0]
              if massco_array[0] > massco_array[-1]:
                  for j in range(-1,-len(variable_array)-1,-1):
                      if j == -1:
                          if variable_array[j] >= thresh:
                              plotlims.append(massco_array[j])
                          else:
                              pass
                      elif (variable_array[j]-thresh)*(variable_array[j+1]-thresh) < 0:
                          plotlims.append(massco_array[j])
                      if j == -len(variable_array):
                          if variable_array[j] >= thresh:
                              plotlims.append(massco_array[j])
                  return plotlims
              else:
                  for j in range(len(variable_array)):
                      if j == 0:
                          if variable_array[j] >= thresh:
                              plotlims.append(massco_array[j])
                          else:
                              pass
                      elif (variable_array[j]-thresh)*(variable_array[j-1]-thresh) < 0:
                          plotlims.append(massco_array[j])
                      if j == len(variable_array)-1:
                          if variable_array[j] >= thresh:
                              plotlims.append(massco_array[j])
                  return plotlims
          # Flag preventing plotting any other variables on an energy generation
          # Kippenhahn plot:
          if engen == True:
              plot = ['dcoeff']
          # This loop gets the mass co-ordinate array and the variable arrays,
          # calls to get the boundaries in order, and populates the contour array.
          ypscoeff = [-1,-1,-1] # this should have same length as plot - quick fix for yps.
          total_massco = []
          for i in range(len(cyclelist)):
  #            print 'CYCLE: ', cyclelist[i]
              massco = self.se.get(cyclelist[i],'mass')
              total_massco.append(max(massco))
              plotlimits=[]
              for j in range(len(plot)):
                  if plot[j][1] == '-' or plot[j][2] == '-':
                      # Assume file has yps, not iso_massf
                      ypsthere = True
                      try: variables = self.se.get(cyclelist[i],'yps')
                      # If this is not the case, do the usual se calls for iso_massf
                      except KeyError:
                          variables = self.se.get(cyclelist[i],plot[j])
                          ypsthere = False
                      # If yps is there, ask which indices correspond to the
                      # elements that are to be plotted, one by one.
                      if ypsthere == True:
                          if ypscoeff[j] == -1:
                              ypscoeff[j] = int(input("What integer is your element "+str(plot[j])+" in the 'yps' array? "))
                          else:
                              pass
                          variables = self.se.get(cyclelist[i],'yps')[:,ypscoeff[j]]
                  else:
                      variables = self.se.get(cyclelist[i],plot[j])
                  plotlims = getlims(variables,thresholds[j],massco)
                  plotlimits.append(plotlims)
              percent = int(i*100/len(cyclelist))
              sys.stdout.flush()
              sys.stdout.write("\rcreating color map " + "...%d%%" % percent)
              for g in range(len(plot)):
                  for k in range(0,len(plotlimits[g]),2):
                      llimit = plotlimits[g][k]
                      ulimit = plotlimits[g][k+1]
                      #if xx[i] >= 0:
                      for f in range(y_res):
                          if llimit<=y[f] and ulimit>y[f]:
                              Z[f,i,g]=1.
                      #else:
                      #   ax.axvline(xx[i],ymin=llimit/m_ini,ymax=ulimit/m_ini,color='#8B8386',alpha=alphas[0],linewidth=0.5)
          # This function determines the adjacent two mass cells to a point in the
          # y-vector (which contains mass co-ordinates centre to surface, split into
          # y_res chunks), returning their index in the mass co-ordinate vector for
          # that timestep/cycle.
          def find_nearest(array,value):
              """
              Returns [lower,upper] indexes locating adjacent mass cells
              (in the massco vector) around y-value (one of y_res points
              equally spaced between centre and surface).
  
              """
              idx=(np.abs(array-value)).argmin()
              lims=np.zeros([2],int)
              if idx == len(array)-1: # SJONES post-mod
                  lims[0] = idx - 1   # SJONES post-mod
                  lims[1] = idx       # SJONES post-mod
                  return lims
              if array[idx] < value:
                  if array[idx]-array[idx+1] < 0.:
                      lims[0] = idx
                      lims[1] = idx-1
                      return lims
                  else:
                      lims[0] = idx
                      lims[1] = idx+1
                      return lims
              elif array[idx] > value:
                  if array[idx]-array[idx+1] < 0.:
                      lims[0] = idx+1
                      lims[1] = idx
                      return lims
                  else:
                      lims[0] = idx-1
                      lims[1] = idx
                      return lims
          # This flag enebles the loop below it to populate the contour array for
          # energy generation. It does not take threshold arguments, as the array
          # contains the log of the energy generation rather than "above" or "below".
          # Because of this, contour boundaries are automatically calculated
          # according to the max energy generation in the model.
          dummy_engen=[]
          engen_signs = []
          if engen == True:
          # Requires eps_nuc array in the data. Produces energy generation contour
          # by linearly interpolating eps_nuc between mass co-ordinates according
          # to the y-resolution:
              for i in range(len(cyclelist)):
  #                print 'CYCLE: ', cyclelist[i]
                  max_energy_gen = 0.
                  min_energy_gen = 0.
                  massco = self.se.get(cyclelist[i],'mass')
                  if len(massco) <= 10:
                      massco=massco[0]
                  dummy_engen = self.se.get(cyclelist[i],netnuc_name)
                  if len(dummy_engen) <= 10:
                      dummy_engen = dummy_engen[0]
                  for f in range(len(dummy_engen)):
                      # make all values absolute, but note in engen_signs which were negative:
                      if dummy_engen[f] == 0.:
                          engen_signs.append(1.)
                      else:
                          engen_signs.append(old_div(dummy_engen[f],abs(dummy_engen[f])))
                      if abs(engen_signs[f]) != 1.:
                          print('engen sign not +/- 1!!')
                          print('engen_signs['+str(f)+'] = ',engen_signs[f])
                          print('dummy_engen[f] = ', dummy_engen[f])
                          sys.exit()
                      dummy_engen[f] = abs(dummy_engen[f])
                  log_epsnuc = np.log10(dummy_engen)
                  # now insert the correct signs again:
                  for f in range(len(log_epsnuc)):
                      log_epsnuc[f] = log_epsnuc[f]*engen_signs[f]
                  #for f in range(len(log_epsnuc)):
                      #if str(log_epsnuc[f]) == 'nan': log_epsnuc[f] = 0.
  #                print log_epsnuc
                  percent = int(i*100/len(cyclelist))
                  sys.stdout.flush()
                  sys.stdout.write("\rcreating color map " + "...%d%%" % percent)
                  for j in range(len(y)):
                      if j == len(y)-1:
                          energy_here = 0.
                      elif j == 0:
                          energy_here = log_epsnuc[-1]
                      elif y[j] > max(massco):
                          energy_here = 0.
                      else:
                          lims = find_nearest(massco,y[j])
                          frac = old_div((y[j]-massco[lims[0]]),(massco[lims[1]]-massco[lims[0]]))
                          energy_here = frac*(log_epsnuc[lims[1]]-log_epsnuc[lims[0]]) + log_epsnuc[lims[0]]
                      if energy_here > max_energy_gen:
                          max_energy_gen = energy_here
                      if energy_here < min_energy_gen:
                          min_energy_gen = energy_here
                      if abs(max_energy_gen) > 100.:
                          print(y[j])
                          print(engen_signs[f], log_epsnuc[f], frac, lims[0], lims[1], massco[lims[0]], massco[lims[1]])
                          print((massco[lims[1]]-massco[lims[0]]), (y[j]-massco[lims[0]]))
                          print(max_energy_gen)
                          print('exit due to energy generation > 100')
                          sys.exit()
  #                    print energy_here
  #                    print max_energy_gen
  #                    if energy_here >= 0.:
                      #Z[j,i,1] = 10**energy_here #SJONES comment
                      if energy_here < 0.:
                          energy_here = 0.
                      Z[j,i,1] = energy_here
  #                    if energy_here < 0.:
  #                        Z[j,i,2] = 10**energy_here
  
  
          # Here we define the colourmap for the energy generation and an array
          # containing a list of colours in which to plot each variable (in the
          # order that the variables appear in "plots") iso_colours is obsolete
          # but was for when we tried plotting isotopes with just their boundary
          # lines as opposed to shading (for clarity). Colourmaps of these choices
          # are written to cmap (array).
          engen_cmap=mpl.cm.get_cmap('Blues')
          engen_cmap.set_under(color='w',alpha=engenalpha)
  
          enloss_cmap=mpl.cm.get_cmap('Reds')
          colours = ['#8B8386','m','g','b']
          iso_colours = ['b','r','y']
          cmap = []
          for i in range(len(plot)):
              cmap.append(mpl.colors.ListedColormap(['w',colours[i]]))
  
          print('plotting contours')
  
          if xllim==0. and xulim==0.:
              ax.axis([float(xx[0]),float(xx[-1]),yllim,yulim])
          else:
              ax.axis([xllim,xulim,yllim,yulim])
  
          # Plot all of the contours. Levels indicates to only plot the shaded
          # regions and not plot the white regions, so that they are essentially
          # transparent. If engen=True, then the energy generation levels
          # (boundary values) are calculated (in dex) from 2 to the maximum in
          # steps of 2.
          for i in range(len(plot)):
              ax.contourf(xx,y,Z[:,:,i],levels=[0.5,1.5],colors=colours[i], alpha=alphas[i])
          if engen == True:
              #ceiling = int(max_energy_gen+1)
              #floor = int(min_energy_gen+1)
              #cburn = ax.contourf(xx,y,Z[:,:,1],cmap=engen_cmap,locator=mpl.ticker.LogLocator(),alpha=engenalpha) # SJONES comment
              cburn = ax.contourf(xx,y,Z[:,:,1],cmap=engen_cmap,alpha=engenalpha,levels=list(range(5,32,5)))
              cbarburn = pl.colorbar(cburn)
  #            if min_energy_gen != 0:
  #                closs = ax.contourf(xx,y,Z[:,:,2],cmap=enloss_cmap,locator=mpl.ticker.LogLocator(),alpha=engenalpha)
  #                cbarloss = pl.colorbar(closs)
  
  #            cbarburn.set_label('$\epsilon_\mathrm{nuc}-\epsilon_{\\nu} \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$, > 0',fontsize=fsize)
              cbarburn.set_label('$log_{10}(\epsilon_\mathrm{nuc}) \; (\mathrm{erg\,g}^{-1}\mathrm{\,s}^{-1})$, > 0',fontsize=fsize)
  
  
          pl.plot(xx,total_massco,color='k')
          pl.text(0.9,0.9,annotation,horizontalalignment='right',transform = ax.transAxes,fontsize=fsize)
          pl.ylabel('$\mathrm{Mass}\;[M_\odot]$',fontsize=fsize-1)
          pl.savefig(outfile)
          print(outfile+' is done.')
          pl.show()