def abu_chart(self, cycle, mass_range=None ,ilabel=True,
                  imlabel=True, imlabel_fontsize=8, imagic=False,
                  boxstable=True, lbound=(-12, 0),
                  plotaxis=[0, 0, 0, 0], show=True, color_map='jet',
                  ifig=None,data_provided=False,thedata=None,
                  savefig=False,drawfig=None,drawax=None,mov=False,
                  path=None):
        '''
        Plots an abundance chart

        Parameters
        ----------
        cycle : string, integer or list
            The cycle we are looking in. If it is a list of cycles,
            this method will then do a plot for each of these cycles
            and save them all to a file.
        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 atomic mass within this range.  This
            will throw an error if this range does not make sence ie
            [45,2] if None, it will plot over the entire range.  The
            default is None.
        ilabel : boolean, optional
            Elemental labels off/on.  The default is True.
        imlabel : boolean, optional
            Label for isotopic masses off/on.  The default is True.
        imlabel_fontsize : integer, optional
            Fontsize for isotopic mass labels.  The default is 8.
        imagic : boolean, optional
            Turn lines for magic numbers off/on.  The default is False.
        boxstable : boolean, optional
            Plot the black boxes around the stable elements.  The
            defaults is True.
        lbound : tuple, optional
            Boundaries for colour spectrum ploted.  The default is
            (-12,0).
        plotaxis : list, optional
            Set axis limit.  If [0, 0, 0, 0] the complete range in (N,Z)
            will be plotted.  It equates to [xMin, xMax, Ymin, Ymax].
            The default is [0, 0, 0, 0].
        show : boolean, optional
            Boolean of if the plot should be displayed.  Useful with
            saving multiple plots using abu_chartMulti.  The default is
            True.
        color_map : string, optional
            Color map according to choices in matplotlib
            (e.g. www.scipy.org/Cookbook/Matplotlib/Show_colormaps).
            The default is 'jet'.
        ifig : integer, optional
            Figure number, if ifig is None it wiil be set to the cycle
            number.  The defaults is None.
        savefig : boolean, optional
            Whether or not to save the figure.
            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
        path: path where to save figure

        '''

        if ifig == None and not mov:
            ifig=cycle

        if type(cycle)==type([]):
            self.abu_chartMulti(cycle, mass_range,ilabel,imlabel,imlabel_fontsize,imagic,boxstable,\
                                lbound,plotaxis,color_map, path=path)
            return
        plotType=self._classTest()

        if mass_range!=None and mass_range[0]>mass_range[1]:
            raise IOError("Please input a proper mass range")

        if plotType=='se':
            if not data_provided:
                cycle=self.se.findCycle(cycle)
#                nin=zeros(len(self.se.A))
#                zin=zeros(len(self.se.Z))
                yin=self.get(cycle, 'iso_massf')
                isom=self.se.isomeric_states
                masses = self.se.get(cycle,'mass')
            else:
                cycle=cycle # why so serious?
                yin=thedata[0]
                isom=self.se.isomeric_states
                masses = thedata[1]

#            for i in xrange(len(nin)):
#                zin[i]=self.se.Z[i]
#                nin[i]=self.se.A[i]-zin[i]
            # SJONES implicit loop instead:
            zin=array([el for el in self.se.Z])
            nin=array([el for el in self.se.A])-zin

            #Test if the mass cell order is inverted
            #and hence mass[-1] the center.
            if masses[0]>masses[-1]:
                #invert
                print('Inverted order of mass cells will be taken into account.')
                yin=yin[::-1]
                masses=masses[::-1]

            if mass_range != None:
                # trim out only the zones needed:
                tmpyps=[]
                masses.sort() # SJ: not sure why this sort if necessary
#                for i in xrange(len(masses)):
#                    if (masses[i] >mass_range[0] and masses[i]<mass_range[1]) or\
#                            (masses[i]==mass_range[0] or masses[i]==mass_range[1]):
#                        tmpyps.append(yin[i])
#                yin=tmpyps
                # find lower and upper indices and slice instead:
                idxl=np.abs(masses-mass_range[0]).argmin()
                if masses[idxl] < mass_range[0]: idxl+=1
                idxu=np.abs(masses-mass_range[1]).argmin()
                if masses[idxu] > mass_range[1]: idxu-=1
                yin=yin[idxl:idxu+1]




            #tmp=zeros(len(yin[0]))
            #for i in xrange(len(yin)):
            #    for j in xrange(len(yin[i])):
            #        tmp[j]+=yin[i][j]

            tmp2=sum(yin,axis=0) # SJONES sum along axis instead of nested loop
            tmp=old_div(tmp2,len(yin))

            yin=tmp

        elif plotType=='PPN':

            ain=self.get('A',cycle)
            zin=self.get('Z',cycle)
            nin=ain-zin
            yin=self.get('ABUNDANCE_MF',cycle)
            isom=self.get('ISOM',cycle)

            if mass_range != None:
                tmpA=[]
                tmpZ=[]
                tmpIsom=[]
                tmpyps=[]
                for i in range(len(nin)):
                    if (ain[i] >mass_range[0] and ain[i]<mass_range[1])\
                    or (ain[i]==mass_range[0] or ain[i]==mass_range[1]):
                        tmpA.append(nin[i])
                        tmpZ.append(zin[i])
                        tmpIsom.append(isom[i])
                        tmpyps.append(yin[i])
                zin=tmpZ
                nin=tmpA
                yin=tmpyps
                isom=tmpIsom

        else:
            raise IOError("This method, abu_chart, is not supported by this class")

        # in case we call from ipython -pylab, turn interactive on at end again
        turnoff=False
        if not show:
            try:
                ioff()
                turnoff=True
            except NameError:
                turnoff=False

        nnmax = int(max(nin))+1
        nzmax = int(max(zin))+1
        nzycheck = zeros([nnmax,nzmax,3])

        for i in range(len(nin)):
            if isom[i]==1:
                ni = int(nin[i])
                zi = int(zin[i])

                nzycheck[ni,zi,0] = 1
                nzycheck[ni,zi,1] = yin[i]

        #######################################################################
        # elemental names: elname(i) is the name of element with Z=i

        elname=self.elements_names

        #### create plot
        ## define axis and plot style (colormap, size, fontsize etc.)
        if plotaxis==[0,0,0,0]:
            xdim=10
            ydim=6
        else:
            dx = plotaxis[1]-plotaxis[0]
            dy = plotaxis[3]-plotaxis[2]
            ydim = 6
            xdim = ydim*dx/dy


        params = {'axes.labelsize':  12,
                  'text.fontsize':   12,
                  'legend.fontsize': 12,
                  'xtick.labelsize': 12,
                  'ytick.labelsize': 12,
                  'text.usetex': True}
        #pl.rcParams.update(params) #May cause Error, someting to do with tex
        if mov:
            fig=drawfig
            fig.set_size_inches(xdim,ydim)
            artists=[]
        else:
            fig=pl.figure(ifig,figsize=(xdim,ydim),dpi=100)
        axx = 0.10
        axy = 0.10
        axw = 0.85
        axh = 0.8
        if mov:
            ax=drawax
        else:
            ax=pl.axes([axx,axy,axw,axh])
        # Tick marks
        xminorlocator = MultipleLocator(1)
        xmajorlocator = MultipleLocator(5)
        ax.xaxis.set_major_locator(xmajorlocator)
        ax.xaxis.set_minor_locator(xminorlocator)
        yminorlocator = MultipleLocator(1)
        ymajorlocator = MultipleLocator(5)
        ax.yaxis.set_major_locator(ymajorlocator)
        ax.yaxis.set_minor_locator(yminorlocator)

        # color map choice for abundances

        cmapa = cm.get_cmap(name=color_map)
        # color map choice for arrows
        cmapr = cm.autumn
        # if a value is below the lower limit its set to white
        cmapa.set_under(color='w')
        cmapr.set_under(color='w')
        # set value range for abundance colors (log10(Y))
        norma = colors.Normalize(vmin=lbound[0],vmax=lbound[1])
        # set x- and y-axis scale aspect ratio to 1
        ax.set_aspect('equal')
        #print time,temp and density on top
        temp = ' '#'%8.3e' %ff['temp']
        time = ' '#'%8.3e' %ff['time']
        dens = ' '#'%8.3e' %ff['dens']

        #May cause Error, someting to do with tex
        '''
        #box1 = TextArea("t : " + time + " s~~/~~T$_{9}$ : " + temp + "~~/~~$\\rho_{b}$ : " \
        #             + dens + ' g/cm$^{3}$', textprops=dict(color="k"))
        anchored_box = AnchoredOffsetbox(loc=3,
                        child=box1, pad=0.,
                        frameon=False,
                        bbox_to_anchor=(0., 1.02),
                        bbox_transform=ax.transAxes,
                        borderpad=0.,
                        )
        ax.add_artist(anchored_box)
        '''
        ## Colour bar plotted

        patches = []
        color = []
        for i in range(nzmax):
            for j in range(nnmax):
                if nzycheck[j,i,0]==1:
                    xy = j-0.5,i-0.5

                    rect = Rectangle(xy,1,1,)

                    # abundance
                    yab = nzycheck[j,i,1]
                    if yab == 0:

                        yab=1e-99


                    col =log10(yab)

                    patches.append(rect)
                    color.append(col)

        p = PatchCollection(patches, cmap=cmapa, norm=norma)
        p.set_array(array(color))
        p.set_zorder(1)
        if mov:
            artist1=ax.add_collection(p)
            artists.append(artist1)
        else:
            ax.add_collection(p)
        if not mov:
            cb = pl.colorbar(p)

            # colorbar label
            cb.set_label('log$_{10}$(X)')

        # plot file name
        graphname = 'abundance-chart'+str(cycle)

        # Add black frames for stable isotopes
        if boxstable:
            for i in range(len(self.stable_el)):
                if i == 0:
                    continue


                tmp = self.stable_el[i]
                try:
                    zz= self.elements_names.index(tmp[0]) #charge
                except:
                    continue

                for j in range(len(tmp)):
                    if j == 0:
                        continue

                    nn = int(tmp[j]) #atomic mass
                    nn=nn-zz

                    xy = nn-0.5,zz-0.5
                    rect = Rectangle(xy,1,1,ec='k',fc='None',fill='False',lw=3.)
                    rect.set_zorder(2)
                    ax.add_patch(rect)


        # decide which array to take for label positions
        iarr = 0

        # plot element labels
        if ilabel:
            for z in range(nzmax):
                try:
                    nmin = min(argwhere(nzycheck[:,z,iarr]))[0]-1
                    ax.text(nmin,z,elname[z],horizontalalignment='center',verticalalignment='center',\
                            fontsize='x-small',clip_on=True)
                except ValueError:
                    continue

        # plot mass numbers
        if imlabel:
            for z in range(nzmax):
                for n in range(nnmax):
                    a = z+n
                    if nzycheck[n,z,iarr]==1:
                        ax.text(n,z,a,horizontalalignment='center',verticalalignment='center',\
                                fontsize=imlabel_fontsize,clip_on=True)

        # plot lines at magic numbers
        if imagic:
            ixymagic=[2, 8, 20, 28, 50, 82, 126]
            nmagic = len(ixymagic)
            for magic in ixymagic:
                if magic<=nzmax:
                    try:
                        xnmin = min(argwhere(nzycheck[:,magic,iarr]))[0]
                        xnmax = max(argwhere(nzycheck[:,magic,iarr]))[0]
                        line = ax.plot([xnmin,xnmax],[magic,magic],lw=3.,color='r',ls='-')
                    except ValueError:
                        dummy=0
                if magic<=nnmax:
                    try:
                        yzmin = min(argwhere(nzycheck[magic,:,iarr]))[0]
                        yzmax = max(argwhere(nzycheck[magic,:,iarr]))[0]
                        line = ax.plot([magic,magic],[yzmin,yzmax],lw=3.,color='r',ls='-')
                    except ValueError:
                        dummy=0

        # set axis limits
        if plotaxis==[0,0,0,0]:

            xmax=max(nin)
            ymax=max(zin)
            ax.axis([-0.5,xmax+0.5,-0.5,ymax+0.5])
        else:
            ax.axis(plotaxis)

        # set x- and y-axis label
        ax.set_xlabel('neutron number (A-Z)')
        ax.set_ylabel('proton number Z')
        if not mov:
            pl.title('Isotopic Chart for cycle '+str(int(cycle)))
        if savefig:
            if path is not None:
                graphname = os.path.join(path, graphname)
            fig.savefig(graphname)
            print(graphname,'is done')
        if show:
            pl.show()
        if turnoff:
            ion()

        if mov:
            return p,artists
        else:
            return