def abu_flux_chart(self, cycle, ilabel=True, imlabel=True,
                       imagic=False, boxstable=True, lbound=(-12,0),
                       plotaxis=[0,0,0,0], which_flux=None, prange=None,
                       profile='charged', show=True):
        '''
        Plots an abundance and flux 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.
        ilabel : boolean, optional
            Elemental labels off/on.  The default is True.
        imlabel : boolean, optional
            Label for isotopic masses off/on.  The default is True.
        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].
        which_flux : integer, optional
            Set to 0 for nucleosynthesis flux plot.  Set to 1 for
            energy flux plot.  Setting wich_flux to 0 is equivelent to
            setting it to 0.  The default is None.
        prange : integer, optional
            Range of fluxes to be considered, if prange is None then
            the plot range is set to 8.  The default is None.
        profile : string, optional
            'charged' is ideal setting to show charged particle
            reactions flow.  'neutron' is ideal setting for neutron
            captures flows.  The default is 'charged'.
        show : boolean, optional
            Boolean of if the plot should be displayed.  Useful with
            saving multiple plots using abu_chartMulti.  The default is
            True.

        '''
        #######################################################################
        #### plot options
        # Set axis limit: If default [0,0,0,0] the complete range in (N,Z) will
        # be plotted, i.e. all isotopes, else specify the limits in
        # plotaxis = [xmin,xmax,ymin,ymax]

        #######################################################################

        # read data file
        #inpfile = cycle
        #ff = fdic.ff(inpfile)
        # with the flux implementation I am not using mass range for now.
        # It may be introduced eventually.
        mass_range = None
        if str(cycle.__class__)=="<type 'list'>":
            self.abu_chartMulti(cycle, mass_range,ilabel,imlabel,imlabel_fontsize,imagic,boxstable,\
                                lbound,plotaxis)
            return
        plotType=self._classTest()

        #if mass_range!=None and mass_range[0]>mass_range[1]:
            #print 'Please input a proper mass range'
            #print 'Returning None'
            #return None

        if plotType=='se':
            cycle=self.se.findCycle(cycle)
            nin=zeros(len(self.se.A))
            zin=zeros(len(self.se.Z))
            for i in range(len(nin)):
                nin[i]=self.se.A[i]
                zin[i]=self.se.Z[i]
            for i in range(len(nin)):
                nin[i]=nin[i]-zin[i]
            yin=self.get(cycle, 'iso_massf')
            isom=self.se.isomeric_states

            masses = self.se.get(cycle,'mass')
            if mass_range != None:
                masses = self.se.get(cycle,'mass')
                masses.sort()

            if mass_range != None:
                tmpyps=[]
                masses = self.se.get(cycle,'mass')
                masses = self.se.get(cycle,'mass')
                masses.sort()
                for i in range(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


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

            tmp=old_div(tmp,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:
            print('This method, abu_chart, is not supported by this class')
            print('Returning None')
            return None
        # 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
        nnmax_plot = nnmax
        nzmax_plot = nzmax
        nzycheck = zeros([nnmax,nzmax,3])
        nzycheck_plot = 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]
                nzycheck_plot[ni,zi,0] = 1



        #######################################################################
        # 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':  15,
                  'text.fontsize':   12,
                  'legend.fontsize': 15,
                  'xtick.labelsize': 15,
                  'ytick.labelsize': 15,
                  'text.usetex': True}
        #pl.rcParams.update(params) #May cause Error, someting to do with tex
        #fig=pl.figure(figsize=(xdim,ydim),dpi=100)
        fig=pl.figure()
        if profile == 'charged':
            ax1 = fig.add_subplot(1, 2, 1)
        elif profile == 'neutron':
            ax1 = fig.add_subplot(2, 1, 1)
        #axx = 0.10
        #axy = 0.10
        #axw = 0.85
        #axh = 0.8
        #ax1=pl.axes([axx,axy,axw,axh])
        # Tick marks
        xminorlocator = MultipleLocator(1)
        xmajorlocator = MultipleLocator(5)
        ax1.xaxis.set_major_locator(xmajorlocator)
        ax1.xaxis.set_minor_locator(xminorlocator)
        yminorlocator = MultipleLocator(1)
        ymajorlocator = MultipleLocator(5)
        ax1.yaxis.set_major_locator(ymajorlocator)
        ax1.yaxis.set_minor_locator(yminorlocator)

        # color map choice for abundances
        #cmapa = cm.jet
        cmapa = cm.summer
        # color map choice for arrows
        cmapr = cm.summer
        # 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
        #ax1.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)
        ax1.add_collection(p)
        cb = pl.colorbar(p)

        # colorbar label
        if profile == 'neutron':
            cb.set_label('log$_{10}$(X)')

        # plot file name
        graphname = 'abundance-flux-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=4.)
                    rect.set_zorder(2)
                    ax1.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
                    nmax = max(argwhere(nzycheck[:,z,iarr]))[0]+1
                    ax1.text(nmin,z,elname[z],horizontalalignment='center',verticalalignment='center',\
                            fontsize='small',clip_on=True)
                    ax1.text(nmax,z,elname[z],horizontalalignment='center',verticalalignment='center',\
                            fontsize='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:
                        ax1.text(n,z,a,horizontalalignment='center',verticalalignment='center',\
                                fontsize='x-small',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 = ax1.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 = ax1.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)
            ax1.axis([-0.5,xmax+0.5,-0.5,ymax+0.5])
        else:
            ax1.axis(plotaxis)

        # set x- and y-axis label
        ax1.set_ylabel('Proton number')
        if profile == 'charged':
            ax1.set_xlabel('Neutron number')
        #pl.title('Isotopic Chart for cycle '+str(int(cycle)))

        #
        # here below I read data from the flux_*****.DAT file.
        #
        file_name = 'flux_'+str(cycle).zfill(5)+'.DAT'
        print(file_name)
        f = open(file_name)
        lines = f.readline()
        lines = f.readlines()
        f.close()

        print_max_flux_in_plot =  False
        # color map choice for fluxes
        #cmapa = cm.jet
        cmapa = cm.autumn
        # color map choice for arrows
        cmapr = cm.autumn
        # starting point of arrow
        coord_x_1 = []
        coord_y_1 = []
        # ending point of arrow (option 1)
        coord_x_2 = []
        coord_y_2 = []
        # ending point of arrow (option 2)
        coord_x_3 = []
        coord_y_3 = []
        # fluxes
        flux_read = []
        flux_log10 = []

        if which_flux == None or which_flux == 0:
            print('chart for nucleosynthesis fluxes [dYi/dt]')
            line_to_read = 9
        elif which_flux == 1:
            print('chart for energy fluxes')
            line_to_read = 10
        elif which_flux > 1:
            print("you have only option 0 or 1, not larger than 1")

        single_line = []
        for i in range(len(lines)):
            single_line.append(lines[i].split())
            coord_y_1.append(int(single_line[i][1]))
            coord_x_1.append(int(single_line[i][2])-coord_y_1[i])
            coord_y_2.append(int(single_line[i][5]))
            coord_x_2.append(int(single_line[i][6])-coord_y_2[i])
            coord_y_3.append(int(single_line[i][7]))
            coord_x_3.append(int(single_line[i][8])-coord_y_3[i])
            try:
                flux_read.append(float(single_line[i][line_to_read]))
            except ValueError: # this is done to avoid format issues like 3.13725-181...
                flux_read.append(1.0E-99)
            flux_log10.append(log10(flux_read[i]+1.0e-99))

        print(file_name,' read!')

        # I need to select smaller sample, with only fluxes inside plotaxis.
        if plotaxis!=[0,0,0,0]:
            coord_y_1_small=[]
            coord_x_1_small=[]
            coord_y_2_small=[]
            coord_x_2_small=[]
            coord_y_3_small=[]
            coord_x_3_small=[]
            flux_log10_small = []
            for i in range(len(flux_log10)):
                I_am_in = 0
                if coord_y_1[i] > plotaxis[2] and coord_y_1[i] < plotaxis[3] and coord_x_1[i] > plotaxis[0] and coord_x_1[i] < plotaxis[1]:
                    I_am_in = 1
                    coord_y_1_small.append(int(coord_y_1[i]))
                    coord_x_1_small.append(int(coord_x_1[i]))
                    coord_y_2_small.append(int(coord_y_2[i]))
                    coord_x_2_small.append(int(coord_x_2[i]))
                    coord_y_3_small.append(int(coord_y_3[i]))
                    coord_x_3_small.append(int(coord_x_3[i]))
                    flux_log10_small.append(flux_log10[i])
                if coord_y_3[i] > plotaxis[2] and coord_y_3[i] < plotaxis[3] and coord_x_3[i] > plotaxis[0] and coord_x_3[i] < plotaxis[1] and I_am_in == 0:
                    I_am_in = 1
                    coord_y_1_small.append(int(coord_y_1[i]))
                    coord_x_1_small.append(int(coord_x_1[i]))
                    coord_y_2_small.append(int(coord_y_2[i]))
                    coord_x_2_small.append(int(coord_x_2[i]))
                    coord_y_3_small.append(int(coord_y_3[i]))
                    coord_x_3_small.append(int(coord_x_3[i]))
                    flux_log10_small.append(flux_log10[i])



        # elemental labels off/on [0/1]
        ilabel = 1

        # label for isotopic masses off/on [0/1]
        imlabel = 1

        # turn lines for magic numbers off/on [0/1]
        imagic = 0

        # flow is plotted over "prange" dex. If flow < maxflow-prange it is not plotted
        if prange == None:
            print('plot range given by default')
            prange = 8.

        #############################################
        #print flux_log10_small
        # we should scale prange on plot_axis range, not on max_flux!
        max_flux = max(flux_log10)
        ind_max_flux = flux_log10.index(max_flux)
        if plotaxis!=[0,0,0,0]:
            max_flux_small = max(flux_log10_small)

        if plotaxis==[0,0,0,0]:
            nzmax = int(max(max(coord_y_1),max(coord_y_2),max(coord_y_3)))+1
            nnmax = int(max(max(coord_x_1),max(coord_x_2),max(coord_x_3)))+1
            coord_x_1_small = coord_x_1
            coord_x_2_small = coord_x_2
            coord_x_3_small = coord_x_3
            coord_y_1_small = coord_y_1
            coord_y_2_small = coord_y_2
            coord_y_3_small = coord_y_3
            flux_log10_small= flux_log10
            max_flux_small  = max_flux
        else:
            nzmax = int(max(max(coord_y_1_small),max(coord_y_2_small),max(coord_y_3_small)))+1
            nnmax = int(max(max(coord_x_1_small),max(coord_x_2_small),max(coord_x_3_small)))+1

        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,)
                    patches.append(rect)


        nzycheck = zeros([nnmax_plot,nzmax,3])
        coord_x_out = zeros(len(coord_x_2_small), dtype='int')
        coord_y_out = zeros(len(coord_y_2_small),dtype='int')
        for i in range(len(flux_log10_small)):
            nzycheck[coord_x_1_small[i],coord_y_1_small[i],0] = 1
            nzycheck[coord_x_1_small[i],coord_y_1_small[i],1] = flux_log10_small[i]
            if coord_x_2_small[i] >= coord_x_3_small[i]:
                coord_x_out[i] = coord_x_2_small[i]
                coord_y_out[i] = coord_y_2_small[i]
                nzycheck[coord_x_out[i],coord_y_out[i],0] = 1
                nzycheck[coord_x_out[i],coord_y_out[i],1] = flux_log10_small[i]
            elif coord_x_2_small[i] < coord_x_3_small[i]:
                coord_x_out[i] = coord_x_3_small[i]
                coord_y_out[i] = coord_y_3_small[i]
                nzycheck[coord_x_out[i],coord_y_out[i],0] = 1
                nzycheck[coord_x_out[i],coord_y_out[i],1] = flux_log10_small[i]
            if flux_log10_small[i]>max_flux_small-prange:
                nzycheck[coord_x_1_small[i],coord_y_1_small[i],2] = 1
                nzycheck[coord_x_out[i],coord_y_out[i],2] = 1

        #### create plot
        if profile == 'charged':
            ax2 = fig.add_subplot(1, 2, 2)
        elif profile == 'neutron':
            ax2 = fig.add_subplot(2, 1, 2)
        # Tick marks
        xminorlocator = MultipleLocator(1)
        xmajorlocator = MultipleLocator(5)
        ax2.xaxis.set_major_locator(xmajorlocator)
        ax2.xaxis.set_minor_locator(xminorlocator)
        yminorlocator = MultipleLocator(1)
        ymajorlocator = MultipleLocator(5)
        ax2.yaxis.set_major_locator(ymajorlocator)
        ax2.yaxis.set_minor_locator(yminorlocator)
        ## 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

        format = 'pdf'
        # set x- and y-axis scale aspect ratio to 1
        #ax2.set_aspect('equal')

        # Add black frames for stable isotopes
        # 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=4.)
                    rect.set_zorder(2)
                    ax2.add_patch(rect)


        apatches = []
        acolor = []
        m = old_div(0.8,prange)
        vmax=ceil(max(flux_log10_small))
        vmin=max(flux_log10_small)-prange
        b=-vmin*m+0.1
        normr = colors.Normalize(vmin=vmin,vmax=vmax)
        ymax=0.
        xmax=0.

        for i in range(len(flux_log10_small)):
            x = coord_x_1_small[i]
            y = coord_y_1_small[i]
            dx = coord_x_out[i]-coord_x_1_small[i]
            dy = coord_y_out[i]-coord_y_1_small[i]
            if flux_log10_small[i]>=vmin:
                arrowwidth = flux_log10_small[i]*m+b
                arrow = Arrow(x,y,dx,dy, width=arrowwidth)
                if xmax<x:
                    xmax=x
                if ymax<y:
                    ymax=y
                acol = flux_log10_small[i]
                apatches.append(arrow)
                acolor.append(acol)
            xy = x-0.5,y-0.5
            rect = Rectangle(xy,1,1,ec='k',fc='None',fill='False',lw=1.)
            patches.append(rect)
            xy = x+dx-0.5,y+dy-0.5
            rect = Rectangle(xy,1,1,ec='k',fc='None',fill='False',lw=1.)
            patches.append(rect)


        p = PatchCollection(patches,norm=0,facecolor='w')
        p.set_zorder(1)
        ax2.add_collection(p)




        a = PatchCollection(apatches, cmap=cmapr, norm=normr)
        a.set_array(array(acolor))
        a.set_zorder(3)
        ax2.add_collection(a)
        cb = pl.colorbar(a)

        # colorbar label
        cb.set_label('log$_{10}$($x$)')
        if profile == 'neutron':
            cb.set_label('log$_{10}$(f)')

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

        # plot element labels
        for z in range(nzmax):
            try:
                nmin = min(argwhere(nzycheck_plot[:,z,iarr-2]))[0]-1
                nmax = max(argwhere(nzycheck_plot[:,z,iarr-2]))[0]+1
                ax2.text(nmin,z,elname[z],horizontalalignment='center',verticalalignment='center',fontsize='small',clip_on=True)
                ax2.text(nmax,z,elname[z],horizontalalignment='center',verticalalignment='center',fontsize='small',clip_on=True)
            except ValueError:
                continue

        # plot mass numbers
        if imlabel:
            for z in range(nzmax):
                for n in range(nnmax_plot):
                    a = z+n
                    if nzycheck_plot[n,z,iarr-2]==1:
                        ax2.text(n,z,a,horizontalalignment='center',verticalalignment='center',fontsize='x-small',clip_on=True)

        # plot lines at magic numbers
        if imagic==1:
            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-2]))[0]
                        xnmax = max(argwhere(nzycheck[:,magic,iarr-2]))[0]
                        line = ax2.plot([xnmin,xnmax],[magic,magic],lw=3.,color='r',ls='-')
                    except ValueError:
                        dummy=0
                if magic<=nnmax:
                    try:
                        yzmin = min(argwhere(nzycheck[magic,:,iarr-2]))[0]
                        yzmax = max(argwhere(nzycheck[magic,:,iarr-2]))[0]
                        line = ax2.plot([magic,magic],[yzmin,yzmax],lw=3.,color='r',ls='-')
                    except ValueError:
                        dummy=0

        # set axis limits
        if plotaxis==[0,0,0,0]:
            ax2.axis([-0.5,xmax+0.5,-0.5,ymax+0.5])
        else:
            ax2.axis(plotaxis)

        # set x- and y-axis label
        ax2.set_xlabel('Neutron number')
        if profile == 'neutron':
            ax2.set_ylabel('Proton number')
        if which_flux == None or which_flux == 0:
            max_flux_label="max flux = "+str('{0:.4f}'.format(max_flux))
        elif which_flux == 1:
            max_flux_label="max energy flux = "+str('{0:.4f}'.format(max_flux))
        if print_max_flux_in_plot:
            ax2.text(plotaxis[1]-1.8,plotaxis[2]+0.1,max_flux_label,fontsize=10.)

        #fig.savefig(graphname)
        print(graphname,'is done')
        if show:
            pl.show()
        if turnoff:
            ion()
        return