docs: def abu_profiles(p,ifig=1,xlm=xlm,ylm=(-8,0),show=False,abunds='All',xaxis=xaxis_type, figsize1=(8,8)):
      '''Four panels of abundance plots
  
      Parameters
      ----------
  
      p : instance
        mesa_profile instance
  
      xlm : tuple
        xlimits: mass_min, mass_max
        
      abus : 'All' plots many 'commonly used' isotopes up to Fe if they are in your mesa output.
      otherwise provide a list of lists of desired abus 
  
      show : Boolean
        False  for batch use
        True   for interactive use
  
      xaxis : character
        Lagrangian    mass is radial mass coordinate
        Eulerian      radius is radial coordinate, in Mm
      '''
  
  
      matplotlib.rc('figure',facecolor='white',figsize=figsize1)
     
      # create subplot structure
      f, ([ax1,ax2],[ax3,ax4]) = pl.subplots(2, 2, sharex=False, sharey=True, figsize=figsize1)
      
      
      # define 4 groups of elements, one for each of the 4 subplots 
      all_isos=[['h1','he3','he4','li6','c12','c13','n13','n14','n15','o16','o17','o18','f19'],['ne20','ne21','ne22','na22','na23','mg24','mg25','mg26','al26','al27','si28','si29','si30'], ['p31', 's32','s33', 's34','s36','cl35','cl37','ar36', 'ar38','ar40', 'k39', 'k40','k41'],
  ['ca40','ca42','ca48','sc45','ti46','ti48','ti50','v50','v51','cr52','cr54','mn55','fe56']]
      
      if abunds == 'All':
          abus=[[],[],[],[]]
          j=0
          for i, row in enumerate(all_isos):
              for iso in row:
                  if iso in p.cols:
                      abus[i].append(iso)
                      j+=1  
                    
                      
          abus1=[]
          abus2 =[[],[],[],[]]
          for l in range(len(abus)):
              for k in range(len(abus[l])):
                  abus1.append(abus[l][k])
  
                  
          is_small_isos = False                
          for i in range(len(abus)):
              if len(abus[i]) < 5:
                  is_small_isos = True
                  print("Missing isotopes from the default list. Distributing the ones you have over the panels.")
          
          if is_small_isos:
              n=4
              quo, rem = divmod(len(abus1), n)
              for i in range(len(abus2)):
                  for k in range(i*quo,(i+1)*quo+rem):
                      abus2[i].append(abus1[k])
           
                      abus = abus2
          #print(abus)        
                  
      else:
          abus = abus    
      ax = [ax1,ax2,ax3,ax4]
      xxx = p.get('radius') if xaxis is "Eulerian" else p.get('mass') 
      mass = p.get('mass')                      # in units of Msun
      radius = p.get('radius')*ast.rsun_cm/1.e8  # in units of Mm
      if xaxis is "Eulerian":
          xxx = radius
  
          if xlm[0] == 0 and xlm[1] == 0:
              indtop = 0
              indbot = len(mass)-1
          else: 
              indbot = np.where(radius>=xlm[0])[0][-1]
              indtop = np.where(radius<xlm[1])[0][0]
  	   
  
          xll = (radius[indbot],radius[indtop])
          xxlabel = "Radius (Mm)"
     
      elif xaxis is "Lagrangian": 
          xxx = mass
          xll = xlm
          xxlabel = "$M / \mathrm{M_{sun}}$"
      else:
          print("Error: don't understand xaxis choice, must be Lagrangian or Eulerian")
  
  
      for i in range(4):
          for thing in abus[i]:
              ind = abus[i].index(thing)
              ax[i].plot(xxx, np.log10(p.get(thing)), ls=u.linestylecb(ind,a,b)[0],\
              marker=u.linestylecb(ind,a,b)[1], color=u.linestylecb(ind,a,b)[2],\
              markevery=50,label=thing)
      # set x and y lims and labels
          ax[i].set_ylim(ylm)
          ax[i].set_xlim(xll)
          ax[i].legend(loc=1)
          ax[i].set_xlabel(xxlabel)
          if i%2 == 0:
              ax[i].set_ylabel('log X')
         # ax[i].set_aspect('equal')
     
      title_str = "Abundance plot: "+'t ='+str(title_format%p.header_attr['star_age'])\
                +' dt ='+str(title_format%p.header_attr['time_step'])\
                +'model number = '+str(int(p.header_attr['model_number']))
      f.suptitle(title_str, fontsize=12)
      f.tight_layout()
      f.subplots_adjust(left=0.1, bottom=0.1, right=0.95, top=0.9, wspace=0, hspace=0.1)
      f.savefig('abuprof'+str(int(p.header_attr['model_number'])).zfill(6)+'.png')