docs: def _kip_vline(self, modstart, modstop, sparse, outfile,
                    xlims=[0.,0.], ylims=[0.,0.], ixaxis='log_time_left',
                    mix_zones=5, burn_zones=50):
          """
          *** DEPRECIATED and hence UNSUPPORTED ***
          This function creates a Kippenhahn plot with energy flux using
          vertical lines, better thermal pulse resolution.
  
          For a more comprehensive plot, your history.data or star.log
          file should contain columns called "mix_type_n", "mix_qtop_n",
          "burn_type_n" and "burn_qtop_n".  The number of columns
          (i.e. the bbiggest value of n) is what goes in the arguments as
          mix_zones and burn_zones.
  
          DO NOT WORRY! if you do not have these columns, just leave the
          default values alone and the script should recognise that you
          do not have these columns and make the most detailed plot that
          is available to you.
  
          Parameters
          ----------
          modstart : integer
              Model from which you want to plot (be careful if your
              history.data or star.log output is sparse...).
          modstop : integer
              Model to which you wish to plot.
          sparse : integer
              x-axis sparsity.
          outfile : string
              'filename + extension' where you want to save the figure.
          xlims, ylims : list, optional
              plot limits, however these are somewhat obsolete now that
              we have modstart and modstop.  Leaving them as 0. is
              probably no slower, and you can always zoom in afterwards
              in mpl.  The default is [0., 0.,].
          ixaxis : string, optional
              Either 'log_time_left', 'age', or 'model_number'.  The
              default "log_time_left".
          mix_zones, burn_zones : integer
              As described above, if you have more detailed output about
              your convection and energy generation boundaries in columns
              mix_type_n, mix_qtop_n, burn_type_n and burn_qtop_n, you
              need to specify the total number of columns for mixing zones
              and burning zones that you have.  Can't work this out from
              your history.data or star.log file?  Check the
              history_columns.list that you used, it'll be the number
              after "mixing regions" and "burning regions".  Can't see
              these columns?  leave it and 2 conv zones and 2 burn zones
              will be drawn using other data that you certainly should
              have in your history.data or star.log file.  The default for
              mix_zones is 5, the defalut for burn_zones is 50.
  
          """
  
  
          xxyy=[self.get('star_age')[modstart:modstop],self.get('star_age')[modstart:modstop]]
          mup = max(float(self.get('star_mass')[0])*1.02,1.0)
          nmodels=len(self.get('model_number')[modstart:modstop])
  
          Msol=1.98892E+33
  
          engenstyle = 'full'
  
          dx = sparse
          x = np.arange(0, nmodels, dx)
  
          btypemax = 20
          btypemin = -20
          btypealpha=0.
  
          ########################################################################
          #----------------------------------plot--------------------------------#
          fig = pl.figure()
  #       fig.set_size_inches(16,9)
          fsize=15
          ax=pl.axes()
  
          if ixaxis == 'log_time_left':
          # log of time left until core collapse
              gage= self.get('star_age')
              lage=np.zeros(len(gage))
              agemin = max(old_div(abs(gage[-1]-gage[-2]),5.),1.e-10)
              for i in np.arange(len(gage)):
                  if gage[-1]-gage[i]>agemin:
                      lage[i]=np.log10(gage[-1]-gage[i]+agemin)
                  else :
                      lage[i]=np.log10(agemin)
              xxx = lage[modstart:modstop]
              print('plot versus time left')
              ax.set_xlabel('$\mathrm{log}_{10}(t^*) \, \mathrm{(yr)}$',fontsize=fsize)
          elif ixaxis =='model_number':
              xxx= self.get('model_number')[modstart:modstop]
              print('plot versus model number')
              ax.set_xlabel('Model number',fontsize=fsize)
          elif ixaxis =='age':
              xxx= old_div(self.get('star_age')[modstart:modstop],1.e6)
              print('plot versus age')
              ax.set_xlabel('Age [Myr]',fontsize=fsize)
          else:
              print('ixaxis must be one of: log_time_left, age or model_number')
              sys.exit()
  
          if xlims == [0.,0.]:
              xlims[0] = xxx[0]
              xlims[1] = xxx[-1]
          if ylims == [0.,0.]:
              ylims[0] = 0.
              ylims[1] = mup
  
  
          print('plotting patches')
          ax.plot(xxx[::dx],self.get('star_mass')[modstart:modstop][::dx],'k-')
  
          print('plotting abund boundaries')
          ax.plot(xxx,self.get('h1_boundary_mass')[modstart:modstop],label='H boundary')
          ax.plot(xxx,self.get('he4_boundary_mass')[modstart:modstop],label='He boundary')
  #       ax.plot(xxx,self.get('c12_boundary_mass')[modstart:modstop],label='C boundary')
  
          ax.axis([xlims[0],xlims[1],ylims[0],ylims[1]])
  
          ax.set_ylabel('Mass [M$_\odot$]')
  
          ########################################################################
  
          try:
              self.get('burn_qtop_1')
          except:
              engenstyle = 'twozone'
          old_percent = 0
          if engenstyle == 'full':
              for i in range(len(x)):
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map1 " + "...%d%%" % percent)
                      old_percent = percent
  
                  for j in range(1,burn_zones+1):
                      ulimit=self.get('burn_qtop_'+str(j))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                      if j==1:
                          llimit=0.0
                      else:
                          llimit=self.get('burn_qtop_'+str(j-1))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                      btype=float(self.get('burn_type_'+str(j))[modstart:modstop][i*dx])
                      if llimit!=ulimit:
                          if btype>0.:
                              #btypealpha = btype/btypemax
                              #ax.axvline(xxx[i*dx],ymin=(llimit-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimit-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=btypealpha)
                              pass
                          if btype<0.:
                              #btypealpha = (btype/btypemin)/5
                              #ax.axvline(xxx[i*dx],ymin=(llimit-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimit-ylims[0])/(ylims[1]-ylims[0]),color='r',alpha=btypealpha)
                              pass
  
              print(' \n')
  
          old_percent = 0
  
          if engenstyle == 'twozone':
              for i in range(len(x)):
                  # writing status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map1 " + "...%d%%" % percent)
                      old_percent = percent
  
                  llimitl1=old_div(self.get('epsnuc_M_1')[modstart:modstop][i*dx],Msol)
                  ulimitl1=old_div(self.get('epsnuc_M_4')[modstart:modstop][i*dx],Msol)
                  llimitl2=old_div(self.get('epsnuc_M_5')[modstart:modstop][i*dx],Msol)
                  ulimitl2=old_div(self.get('epsnuc_M_8')[modstart:modstop][i*dx],Msol)
                  llimith1=old_div(self.get('epsnuc_M_2')[modstart:modstop][i*dx],Msol)
                  ulimith1=old_div(self.get('epsnuc_M_3')[modstart:modstop][i*dx],Msol)
                  llimith2=old_div(self.get('epsnuc_M_6')[modstart:modstop][i*dx],Msol)
                  ulimith2=old_div(self.get('epsnuc_M_7')[modstart:modstop][i*dx],Msol)
                  # lower thresh first, then upper thresh:
                  #if llimitl1!=ulimitl1:
                      #ax.axvline(xxx[i*dx],ymin=(llimitl1-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimitl1-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=1.)
                  #if llimitl2!=ulimitl2:
                      #ax.axvline(xxx[i*dx],ymin=(llimitl2-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimitl2-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=1.)
                  #if llimith1!=ulimith1:
                      #ax.axvline(xxx[i*dx],ymin=(llimith1-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimith1-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=4.)
                  #if llimith2!=ulimith2:
                      #ax.axvline(xxx[i*dx],ymin=(llimith2-ylims[0])/(ylims[1]-ylims[0]),ymax=(ulimith2-ylims[0])/(ylims[1]-ylims[0]),color='b',alpha=4.)
  
              print(' \n')
  
          mixstyle = 'full'
          try:
              self.get('mix_qtop_1')
          except:
              mixstyle = 'twozone'
  
          old_percent = 0
  
          if mixstyle == 'full':
              for i in range(len(x)):
              # writing reading status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map2 " + "...%d%%" % percent)
                      old_percent = percent
  
                  for j in range(1,mix_zones+1):
                      ulimit=self.get('mix_qtop_'+str(j))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                      if j==1:
                          llimit=0.0
                      else:
                          llimit=self.get('mix_qtop_'+str(j-1))[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                      mtype=self.get('mix_type_'+str(j))[modstart:modstop][i*dx]
                      if llimit!=ulimit:
                          if mtype == 1:
                              ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=3., linewidth=.5)
  
              print(' \n')
  
          old_percent = 0
  
          if mixstyle == 'twozone':
              for i in range(len(x)):
              # writing reading status
                  percent = int(i*100/(len(x) - 1))
                  if percent >= old_percent + 5:
                      sys.stdout.flush()
                      sys.stdout.write("\r creating color map2 " + "...%d%%" % percent)
                      old_percent = percent
  
                  ulimit=self.get('conv_mx1_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  llimit=self.get('conv_mx1_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  if llimit!=ulimit:
                      ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=5.,linewidth=.5)
                  ulimit=self.get('conv_mx2_top')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  llimit=self.get('conv_mx2_bot')[modstart:modstop][i*dx]*self.get('star_mass')[modstart:modstop][i*dx]
                  if llimit!=ulimit:
                      ax.axvline(xxx[i*dx],ymin=old_div((llimit-ylims[0]),(ylims[1]-ylims[0])),ymax=old_div((ulimit-ylims[0]),(ylims[1]-ylims[0])),color='k',alpha=3.,linewidth=.5)
  
              print(' \n')
  
          print('engenstyle was ', engenstyle)
          print('mixstyle was ', mixstyle)
          print('\n finished preparing color map')
  
          #fig.savefig(outfile)
          pl.show()