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')