docs: def iso_abund(self, cycle, stable=False, amass_range=None,
                    mass_range=None, ylim=[0,0], ref=-1, show=True,
                    log_logic=True, decayed=False, color_plot=True,
                    grid=False, point_set=1, include_title=False,
                    data_provided=False,thedata=None, verbose=True,
                    mov=False,drawfig=None,drawax=None,show_names=True,
                    label=None,colour=None,elemaburtn=False,mypoint=None,plot_type=['-','--','-.',':','-']):
          '''
          plot the abundance of all the chemical species
  
          Parameters
          ----------
          cycle : string, integer or list
              The cycle of interest.  If it is a list of cycles, this
              method will do a plot for each cycle and save them to a
              file.
          stable : boolean, optional
              A boolean of whether to filter out the unstables.  The
              defaults is False.
          amass_range : list, optional
              A 1x2 array containing the lower and upper atomic mass
              range.  If None plot entire available atomic mass range.
              The default is None.
          mass_range : list, optional
              A 1x2 array containing the lower and upper mass range.  If
              this is an instance of abu_vector this will only plot
              isotopes that have an atominc mass within this range.  This
              will throw an error if this range does not make sense ie
              [45,2].  If None, it will plot over the entire range.  The
              defaults is None.
          ylim : list, optional
              A 1x2 array containing the lower and upper Y limits.  If
              it is [0,0], then ylim will be determined automatically.
              The default is [0,0].
          ref : integer or list, optional
              reference cycle.  If it is not -1, this method will plot
              the abundences of cycle devided by the cycle of the same
              instance given in the ref variable.  If ref is a list it
              will be interpreted to have two elements:
              ref=['dir/of/ref/run',cycle] which uses a refernece cycle
              from another run.  If any abundence in the reference cycle
              is zero, it will replace it with 1e-99.  The default is -1.
          show : boolean, optional
              Boolean of if the plot should be displayed.  The default is
              True.
          log_logic : boolean, optional
              Plot abundances in log scale or linear.  The default is
              True.
          decayed : boolean, optional
              If True plot decayed distributions, else plot life
              distribution.  The default is False.
          color_plot : boolean, optional
              Color dots and lines [True/False].  The default is True.
          grid : boolean, optional
              print grid.  The default is False.
          point_set : integer, optional
              Set to 0, 1 or 2 to select one of three point sets, useful
              for multiple abundances or ratios in one plot.  The defalult
              is 1.
          include_title : boolean, optional
              Include a title with the plot.  The default is False.
          drawfig, drawax, mov : optional, not necessary for user to set these variables
              The figure and axes containers to be drawn on, and whether or not a movie is
              being made (only True when se.movie is called, which sets mov to True
              automatically
          elemaburtn : boolean, private
              If true, iso_abund() returns after writing self.***_iso_to_plot for
              use with other plotting routines.f
          mypoint : string, optional
              fix the marker style of all the points in this plot to one type, given
              as a string. If None, multiple point styles are used as per point_set.
              The default is None
          '''
  
          plotType=self._classTest()
          if str(cycle.__class__)=="<type 'list'>":
              self.iso_abundMulti(cycle, stable,amass_range,mass_range,ylim,ref,
                   decayed,include_title,color_plot=color_plot,grid=False,point_set=point_set)
              return
  
          if mass_range!=None and mass_range[0]>mass_range[1]:
              print('Please input a proper mass range')
              print('Returning None')
              return None
          if amass_range!=None and amass_range[0]>amass_range[1]:
              print('Please input a proper Atomic mass range')
              print('Returning None')
              return None
          if plotType=='se':
              if decayed:
                  print('Decay option not yet implemented for mppnp - but it is easy do! Consider investing the time!')
                  return None
  
  
              # get things as arrays
              if not data_provided:
                  cycle=self.se.findCycle(cycle)
                  a_iso_to_plot   = array(self.se.A)
                  abunds          = self.get(cycle,'iso_massf')
                  isotope_to_plot = array(self.se.isotopes)
                  z_iso_to_plot   = array(self.se.Z)
                  isomers_to_plot = array(self.se.isomeric_states)
                  if ref >-1:
                      ref=self.se.findCycle(ref)
                      abundsRef=self.se.get(ref,'iso_massf')
                  masses = self.se.get(cycle,'mass')
              else:
                  cycle=cycle # why so serious?
                  a_iso_to_plot   = array(self.se.A)
                  abunds          = thedata[0]
                  isotope_to_plot = array(self.se.isotopes)
                  z_iso_to_plot   = array(self.se.Z)
                  isomers_to_plot = array(self.se.isomeric_states)
  
                  if ref >-1:
                      raise IOError("No. It's not ready yet.")
                      #ref=self.se.findCycle(ref)
                      #abundsRef=self.se.get(ref,'iso_massf')
                  masses = thedata[1]
  
              if mass_range == None:
                  if verbose:
                      print('Using default mass range')
                  mass_range = [min(masses),max(masses)]
              masses.sort()
              mass_range.sort()
  
              if amass_range == None:
                  amass_range=[int(min(a_iso_to_plot)),int(max(a_iso_to_plot))]
  
              # remove neutrons - this could move in the non- se/PPN specific part below
              if 0 in z_iso_to_plot:
                  ind_neut        = where(z_iso_to_plot==0)[0][0]
                  a_iso_to_plot   = delete(a_iso_to_plot,ind_neut)
                  z_iso_to_plot   = delete(z_iso_to_plot,ind_neut)
                  isomers_to_plot = delete(isomers_to_plot,ind_neut)
                  isotope_to_plot = delete(isotope_to_plot,ind_neut)
                  abunds = delete(abunds,ind_neut,1)
                  if ref >-1:
                      abundsRef = delete(abundsRef,ind_neut,1)
  
              # extract amass_range
              acon=(a_iso_to_plot>=amass_range[0]) & (a_iso_to_plot<=amass_range[1])
              isomers_to_plot = isomers_to_plot[acon]
              isotope_to_plot = isotope_to_plot[acon]
              z_iso_to_plot   = z_iso_to_plot[acon]
              abunds          = abunds.T[acon].T
              if ref >-1:
                  abundsRef = abundsRef.T[acon].T
              a_iso_to_plot   = a_iso_to_plot[acon]
              el_iso_to_plot  = array([x.split('-')[0] for x in isotope_to_plot.tolist()])
              # apply mass range
              if mass_range == None:
                  if verbose:
                      print('Using default mass range')
                  mass_range = [min(masses),max(masses)]
              mass_range.sort()
              aabs = []
              if ref >-1:
                  cyc  = [cycle,ref]
                  abus = [abunds,abundsRef]
              else:
                  cyc  = [cycle]
                  abus = [abunds]
              for cc,aa in zip(cyc,abus):
                  if not data_provided:
                      masses = self.se.get(cc,'mass')
                  else:
                      masses=masses # why so serious?
                  masses.sort()
                  dmass  = masses[1:] - masses[:-1]    # I should check the grid definition
                  dmass  = append(dmass,0.)
                  mcon   = (masses>=mass_range[0]) & (masses<=mass_range[1])
                  dmass  = dmass[mcon]
                  aa = aa[mcon]
  
                  # average over mass range:
                  aa = (aa.T*dmass).T.sum(0)
                  aa = old_div(aa, (mass_range[1] - mass_range[0]))
                  # abunds has now length of isotope_to_plot
                  aabs.append(aa)
              if ref >-1:
                  abunds = old_div(aabs[0],(aabs[1]+1.e-99))
              else:
                  abunds = aabs[0]
  
              self.a_iso_to_plot=a_iso_to_plot
              self.isotope_to_plot=isotope_to_plot
              self.z_iso_to_plot=z_iso_to_plot
              self.el_iso_to_plot=el_iso_to_plot
              self.abunds=abunds
              self.isomers_to_plot=isomers_to_plot
  
              if elemaburtn: return
  #                        self.isotopes = self.se.isotopes
  
          elif plotType=='PPN':
              print("This method adds the following variables to the instance:")
              print("a_iso_to_plot      mass number of plotted range of species")
              print("isotope_to_plot    corresponding list of isotopes")
              print("z_iso_to_plot      corresponding charge numbers")
              print("el_iso_to_plot     corresponding element names")
              print("abunds             corresponding abundances")
              print("isom               isomers and their abundance")
  
              self.get(cycle,decayed=decayed)
              if ref is not -1:
                  if type(ref) is list: # reference cycle from other run
                      import ppn
                      pp=ppn.abu_vector(ref[0])
                      abunds_pp=pp.get(ref[1],decayed=decayed)
                      self.abunds=old_div(self.abunds,pp.abunds)
                  else:
                      abunds=self.abunds
                      self.get(ref,decayed=decayed)
                      self.abunds=old_div(abunds,(self.abunds+1.e-99))
              if amass_range == None:
                  amass_range=[min(self.a_iso_to_plot),max(self.a_iso_to_plot)]
  
              aa=ma.masked_outside(self.a_iso_to_plot,amass_range[0],amass_range[1])
              isotope_to_plot=ma.array(self.isotope_to_plot,mask=aa.mask).compressed()
              z_iso_to_plot=ma.array(self.z_iso_to_plot,mask=aa.mask).compressed()
              el_iso_to_plot=ma.array(self.el_iso_to_plot,mask=aa.mask).compressed()
              abunds=ma.array(self.abunds,mask=aa.mask).compressed()
              a_iso_to_plot=aa.compressed()
              isomers_to_plot=[]
              for i in range(len(self.isom)):
                  if int(self.isom[i][0].split('-')[1])>100:
                      isomers_to_plot.append(self.isom[i])
  
              self.a_iso_to_plot=a_iso_to_plot
              self.isotope_to_plot=isotope_to_plot
              self.z_iso_to_plot=z_iso_to_plot
              self.el_iso_to_plot=el_iso_to_plot
              self.abunds=abunds
              self.isomers_to_plot=isomers_to_plot
          else:
              print('This method, iso_abund, is not supported by this class')
              print('Returning None')
              return None
  
          if verbose:
              print('Using the following conditions:')
              if plotType=='se':
                  print('\tmass_range:', mass_range[0], mass_range[1])
              print('\tAtomic mass_range:', amass_range[0], amass_range[1])
              print('\tcycle:           ',cycle)
              print('\tplot only stable:',stable)
              print('\tplot decayed:    ',decayed)
  
  
          if stable: # remove unstables:
              # For the element that belongs to the isotope at index 5 in isotope_to_plot
              # (C-12) the following gives the mass numbers of stable elements:
              # self.stable_el[self.stable_names.index(el_iso_to_plot[5])][1:]
              ind_delete=[]
              for i in range(len(isotope_to_plot)):
                  if a_iso_to_plot[i] not in self.stable_el[self.stable_names.index(el_iso_to_plot[i])][1:]:
                      ind_delete.append(i)
              a_iso_to_plot   = delete(a_iso_to_plot,  ind_delete)
              z_iso_to_plot   = delete(z_iso_to_plot,  ind_delete)
              isomers_to_plot = delete(isomers_to_plot,ind_delete)
              isotope_to_plot = delete(isotope_to_plot,ind_delete)
              el_iso_to_plot  = delete(el_iso_to_plot, ind_delete)
              abunds          = delete(abunds,         ind_delete)
  
  #        el_list=[] # list of elements in el_iso_to_plot
  #
  #        for el in self.elements_names:
  #            if el in el_iso_to_plot:
  #                el_list.append(el)
  
          # SJONES implicit loop:
          el_list = [el for el in self.elements_names if el in el_iso_to_plot]
  
          abund_plot = []  # extract for each element an abundance and associated
          mass_num   = []  # mass number array, sorted by mass number
  
          for el in el_list:
              numbers = a_iso_to_plot[(el_iso_to_plot==el)]
              abund_plot.append(abunds[(el_iso_to_plot==el)][argsort(numbers)])
              mass_num.append(sort(numbers))
          # now plot:
          #plot_type = ['-','--','-.',':','-'] ##now implemented as an arg
          print(plot_type)
          while len(plot_type)<=4:
              plot_type.append('')
              print(plot_type)
          pl_index = 0
          if mypoint is None:
              points = [['o','^','p','h','*'],['x','+','D','>','s'],['H','v','<','*','3']]
          else:
              points = [ [mypoint]*5 , [mypoint]*5 , [mypoint]*5]
          if color_plot:
              colors = ['g','r','c','m','k']
          elif colour is not None:
              colors = [colour]*5
          else:
              colors = ['k','k','k','k','k']
  
          ylim1 =  1.e99
          ylim2 = -1.e99
  
          # initialise movie-related things:
          if mov:
              artists=[]
              ax=drawax
              fig=drawfig
          elif drawax is not None:
              ax=drawax
          else:
              ax=pl.axes()
  
          if drawfig is not None:
              fig=drawfig
  
          for j in range(len(abund_plot)):        #Loop through the elements of interest
  #            for l in xrange(len(abund_plot[j])):
  #                if abund_plot[j][l] == 0:
  #                    abund_plot[j][l] = 1e-99
  
              abund_plot[j] = np.maximum(abund_plot[j],1.e-99) # SJONES instead of looping
  
  #            a_dum=zeros(len(abund_plot[j]))   # this I (FH) have to do because for some
              if log_logic == False:            # reason log10(abu_abund[j]) does not work
                  a_dum = abund_plot[j]         # although abu_abund[j] is a numpy array?!?
              else:
  #                for ii in range(len(abund_plot[j])):
  #                    a_dum[ii]=log10(abund_plot[j][ii])
                  a_dum=np.log10(abund_plot[j]) # SJONES this seems to work fine for me
              if type(colors[0]) is str:
                  this_label=str(colors[pl_index]+points[point_set][pl_index]+\
                             plot_type[pl_index])
              else:
                  this_label=None
              if mov:
                  artist1,=ax.plot(mass_num[j],a_dum,this_label,markersize=6,
                                   markeredgecolor='None')
              else:
                  if this_label is not None:
                      if label is not None and j==0:
                          pl.plot(mass_num[j],a_dum,this_label,markersize=6,
                                  label=label,markeredgecolor='None')
                          pl.legend(loc='best').draw_frame(False)
                      else:
                          pl.plot(mass_num[j],a_dum,this_label,markersize=6,
                                  markeredgecolor='None')
                  else:
                      if label is not None and j==0:
                          pl.plot(mass_num[j],a_dum,
                                  color=colors[pl_index],
                                  marker=points[point_set][pl_index],
                                  linestyle=plot_type[pl_index],
                                  markersize=6,label=label,
                                  markeredgecolor='None')
                          pl.legend(loc='best').draw_frame(False)
                      else:
                          pl.plot(mass_num[j],a_dum,
                                  color=colors[pl_index],
                                  marker=points[point_set][pl_index],
                                  linestyle=plot_type[pl_index],
                                  markersize=6,markeredgecolor='None')
  
              abu_max = max(a_dum)
              max_index=where(a_dum==abu_max)[0][0]
              coordinates=[mass_num[j][max_index],abu_max]
              if mov:
                  artist2=ax.text(coordinates[0]+0.1,1.05*coordinates[1],el_list[j],clip_on=True)
              else:
                  if show_names:
  #                    pl.text(coordinates[0]+0.1,1.05*coordinates[1],el_list[j],clip_on=True)
                      pl.text(coordinates[0],np.log10(2.2*10.**coordinates[1]),
                              el_list[j],clip_on=True,
                              horizontalalignment='center')
  
              pl_index+=1
              if pl_index > 4:
                  pl_index = 0
              ylim1=min(ylim1,min(a_dum))
              ylim2=max(ylim2,max(a_dum))
  
              if mov:
                  artists.extend([artist1,artist2])
  
          # now trimming the ylims
          if log_logic:
              dylim=0.05*(ylim2-ylim1)
              ylim1 = ylim1 -dylim
              ylim2 = ylim2 +dylim
              if ref is not -1:
                  ylim2 = min(ylim2,4)
                  ylim1 = max(ylim1,-4)
              else:
                  ylim2 = min(ylim2,0.2)
                  ylim1 = max(ylim1,-13)
          else:
              ylim1 = ylim1 *0.8
              ylim2 = ylim2 *1.1
          if include_title:
              if plotType=='se':
                  if ref == -1:
                      title = str('Range %4.2f' %mass_range[0]) + str('-%4.2f' %mass_range[1]) +\
                          str(' for cycle %d' %int(cycle))
                  else:
                      title = str('Range %4.2f' %mass_range[0]) + \
                          str('-%4.2f' %mass_range[1]) + str(' for cycle %d' %int(cycle))+\
                          str(' relative to cycle %d'  %int(ref))
              else:
                  if ref == -1:
                      title = str('Cycle %d' %int(cycle))
                  else:
                      title = str('Cycle %d' %int(cycle))+\
                              str(' relative to cycle %d'  %int(ref))
              print("including title: ...")
              if mov:
                  artist1,=ax.title(title)
                  artists.append(artist1)
              else:
                  pl.title(title)
          if ylim[0] == 0 and ylim[1] == 0:
              pl.ylim(ylim1,ylim2)
          else:
              pl.ylim(ylim[0],ylim[1])
          pl.xlim([amass_range[0]-.5,amass_range[1]+.5])
          pl.xlabel('mass number (A)',fontsize=14)
          if ref is not -1:
              if log_logic:
                  pl.ylabel(r'log abundance ratio',fontsize=14)
              else:
                  pl.ylabel(r'abundance ratio',fontsize=14)
          else:
              if log_logic:
                  pl.ylabel(r'log mass fraction ',fontsize=14)
              else:
                  pl.ylabel(r'mass fraction',fontsize=14)
  
  
          if amass_range != None:
              minimum_mass = amass_range[0]
              maximum_mass = amass_range[1]
  
          elif mass_range != None:
              minimum_mass = mass_range[0]
              maximum_mass = mass_range[1]
  
          else:
              minimum_mass = 0
              maximum_mass = 200
  
          if log_logic == False:
              if mov:
                  artist1,=ax.plot([amass_range[0]-.5,amass_range[1]+.5],[1,1],'k-')
                  artists.append(artist1)
              else:
                  pl.plot([amass_range[0]-.5,amass_range[1]+.5],[1,1],'k-')
          else:
              if mov:
                  artist1,=ax.plot([amass_range[0]-.5,amass_range[1]+.5],[0,0],'k-')
                  artists.append(artist1)
              else:
                  pl.plot([amass_range[0]-.5,amass_range[1]+.5],[0,0],'k-')
  
          labelsx=[]
          if (maximum_mass-minimum_mass) > 100:
              delta_labelsx = 10
          else:
              delta_labelsx = 5
  
          iii = amass_range[0]%delta_labelsx
          if iii == 0:
              labelsx.append(str(amass_range[0]))
          else:
              labelsx.append(' ')
          iii = iii+1
          kkk = 0
          for label1 in range(amass_range[1]-amass_range[0]):
              if iii == 5:
                  kkk = kkk+1
                  labelsx.append(str((iii*kkk)+amass_range[0]-(amass_range[0]%5)))
                  iii = 0
                  iii = iii+1
              else:
                  labelsx.append(' ')
                  iii = iii+1
  
          if delta_labelsx == 5:
              xticks = arange(amass_range[0],amass_range[1],1)
              pl.xticks(xticks,labelsx)
          else:
              pl.xticks()
  
          # SJONES moved the pl.grid and pl.show to the very end
          if grid:
              pl.grid()
          if show:
              pl.show()
  
  ##!!FOR!!###### print 'LEN LABELS= ', len(labelsx)
  ##DEBUGGING####
  ####!!!######## for bbb in range (len(labelsx)):
  ###############     print labelsx[bbb]
          if mov:
              return artists