def flux_chart(file_name, plotaxis, plot_type, which_flux=None,
               I_am_the_target=None, prange=None):
    '''
    Plots a chart with fluxes

    Parameters
    ----------
    file_name : string
        Name of the file of fluxes we are looking at.
    plotaxis : list
        [xmin, xmax, ymin, ymax], where on x axis there is neutron
        number and on y axis there is Z.
    plot_types : integer
        Set to 0 for standard flux plot.  Set to 1 if fluxes focused
        on one specie.
    which_flux : integer, optional
        Set to 0 for nucleosynthesis flux plot.  Set to 1 is for energy
        flux plot.  Seting to None is the same a 0.  The default is
        None.
    I_am_the_target : list, optional
        A 2xArray used only if plot_type=1, and is given by [neutron
        number, proton number].  The default is None.
    prange : integer, optional
        The range of fluxes to be considered.  If prange is None, then
        8 fluxes are consdered.  The default is None.

    Notes
    -----
    This script is terribly slow and needs to be improved.  For now I
    put here in data_plot:

    [1]: import data_plot

    [2]: data_plot.flux_chart('file_name', [xmin, xmax, ymin, ymax],
                              int, which_flux, I_am_the_target, prange)

    The pdf is created, but an error bumped up and the gui is empty.
    To avoid this, I had to set 'text.usetex': False.  See below.  Also,
    for the same reason no label in x axys is written using
    'text.usetex': True.

    Note also that the GUI works really slow with this plot.  So, we
    need to optimize from the graphic point of view.  This need to be
    included in ppn.py I think, and set in multi option too, in case
    we want to read more flux files at the same time.

    Finally, you need to have stable.dat to read in to make it work ...

    '''

    import numpy as np
    import matplotlib.pyplot as plt
    #from matplotlib.mpl import colors,cm # deppreciated in mpl ver 1.3
                                          # use line below instead
    from matplotlib import colors,cm
    from matplotlib.patches import Rectangle, Arrow
    from matplotlib.collections import PatchCollection
    from matplotlib.offsetbox import AnchoredOffsetbox, TextArea
    import sys

    print_max_flux_in_plot =  True


    f = open(file_name)
    lines = f.readline()
    lines = f.readlines()
    f.close()

    # 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 == 2:
        print('chart for timescales')
        line_to_read = 11
    elif which_flux > 2:
        print("you have only option 0, 1 or 2, not larger than 2")

    single_line = []
    for i in range(len(lines)):
        single_line.append(lines[i].split())
        coord_y_1.append(float(single_line[i][1]))
        coord_x_1.append(float(single_line[i][2])-coord_y_1[i])
        coord_y_2.append(float(single_line[i][5]))
        coord_x_2.append(float(single_line[i][6])-coord_y_2[i])
        coord_y_3.append(float(single_line[i][7]))
        coord_x_3.append(float(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(np.log10(flux_read[i]+1.0e-99))

    print('file read!')

    # I need to select smaller sample, with only fluxes inside plotaxis.
    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(coord_y_1[i])
            coord_x_1_small.append(coord_x_1[i])
            coord_y_2_small.append(coord_y_2[i])
            coord_x_2_small.append(coord_x_2[i])
            coord_y_3_small.append(coord_y_3[i])
            coord_x_3_small.append(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(coord_y_1[i])
            coord_x_1_small.append(coord_x_1[i])
            coord_y_2_small.append(coord_y_2[i])
            coord_x_2_small.append(coord_x_2[i])
            coord_y_3_small.append(coord_y_3[i])
            coord_x_3_small.append(coord_x_3[i])
            flux_log10_small.append(flux_log10[i])


    if plot_type == 1:
        print('I_am_the_target=',I_am_the_target)
        #I_am_the_target = [56.-26.,26.]
    # here below need for plotting
    # plotaxis = [xmin,xmax,ymin,ymax]
    #plotaxis=[1,20,1,20]
    #plotaxis=[0,0,0,0]

    # 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.

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

    # 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)
    max_flux_small = max(flux_log10_small)
    min_flux = min(flux_log10)
    ind_min_flux = flux_log10.index(min_flux)
    min_flux_small = min(flux_log10_small)

    #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
    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


    nzycheck = np.zeros([nnmax,nzmax,3])
    #coord_x_out = np.zeros(len(coord_x_2))
    #coord_y_out = np.zeros(len(coord_y_2))
    #for i in range(len(flux_log10)):
    #       nzycheck[coord_x_1[i],coord_y_1[i],0] = 1
    #       nzycheck[coord_x_1[i],coord_y_1[i],1] = flux_log10[i]
    #       if coord_x_2[i] >= coord_x_3[i]:
    #                coord_x_out[i] = coord_x_2[i]
    #                coord_y_out[i] = coord_y_2[i]
    #               nzycheck[coord_x_out[i],coord_y_out[i],0] = 1
    #               nzycheck[coord_x_out[i],coord_y_out[i],1] = flux_log10[i]
    #        elif coord_x_2[i] < coord_x_3[i]:
    #                coord_x_out[i] = coord_x_3[i]
    #                coord_y_out[i] = coord_y_3[i]
    #               nzycheck[coord_x_out[i],coord_y_out[i],0] = 1
    #               nzycheck[coord_x_out[i],coord_y_out[i],1] = flux_log10[i]
    #       if flux_log10[i]>max_flux-prange:
    #               nzycheck[coord_x_1[i],coord_y_1[i],2] = 1
    #               nzycheck[coord_x_out[i],coord_y_out[i],2] = 1
    coord_x_out = np.zeros(len(coord_x_2_small))
    coord_y_out = np.zeros(len(coord_y_2_small))
    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 which_flux == None or which_flux < 2 and 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
        elif which_flux == 2 and flux_log10_small[i]<min_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

    #######################################################################
    # elemental names: elname(i) is the name of element with Z=i
    elname= ('none','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe',
    'Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb',
    'Te', 'I','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os',
    'Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U','Np','Pu')


    #### 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

    format = 'pdf'
    # note that I had to set 'text.usetex': False, to avoid Exception in Tkinter callback.
    # and to make the GUI work properly. Why? some missing package?
    params = {'axes.labelsize':  15,
      'text.fontsize':   15,
      'legend.fontsize': 15,
      'xtick.labelsize': 15,
      'ytick.labelsize': 15,
      'text.usetex': False}
    plt.rcParams.update(params)
    fig=plt.figure(figsize=(xdim,ydim),dpi=100)
    axx = 0.10
    axy = 0.10
    axw = 0.85
    axh = 0.8
    ax=plt.axes([axx,axy,axw,axh])

    # color map choice for abundances
    cmapa = cm.jet
    # color map choice for arrows
    if which_flux == None or which_flux < 2:
        cmapr = cm.autumn
    elif  which_flux == 2:
        cmapr = cm.autumn_r
    # 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=-20,vmax=0)
    # 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']

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

    # Add black frames for stable isotopes
    f = open('stable.dat')

    head = f.readline()
    stable = []

    for line in f.readlines():
        tmp = line.split()
        zz = int(tmp[2])
        nn = int(tmp[3])
        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)

    apatches = []
    acolor = []
    m = old_div(0.8,prange)#0.8/prange
    if which_flux == None or which_flux < 2:
        vmax=np.ceil(max(flux_log10_small))
        vmin=max(flux_log10_small)-prange
        b=-vmin*m+0.1
    elif which_flux == 2:
        vmin=min(flux_log10_small)
        vmax=np.ceil(min(flux_log10_small)+prange)
        b=vmax*m+0.1
    if which_flux == None or which_flux < 3:
        normr = colors.Normalize(vmin=vmin,vmax=vmax)
        print('vmin and vmax =',vmin,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 plot_type == 0:
            if which_flux == None or which_flux < 2:
                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)
            elif which_flux == 2:
                if flux_log10_small[i]<=vmax:
                    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)
        elif plot_type == 1 and which_flux != 2:
            if x==I_am_the_target[0] and y==I_am_the_target[1] and 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)
            if x+dx==I_am_the_target[0] and y+dy==I_am_the_target[1] and 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)
        elif plot_type == 1 and which_flux == 2:
            if x==I_am_the_target[0] and y==I_am_the_target[1] and flux_log10_small[i]<=vmax:
                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)
            if x+dx==I_am_the_target[0] and y+dy==I_am_the_target[1] and flux_log10_small[i]<=vmax:
                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)

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

    #for i in range(len(flux_log10)):
    #       x = coord_x_1[i]
    #       y = coord_y_1[i]
    #       dx = coord_x_out[i]-coord_x_1[i]
    #       dy = coord_y_out[i]-coord_y_1[i]
    #       if plot_type == 0:
    #               if flux_log10[i]>=vmin:
    #                       arrowwidth = flux_log10[i]*m+b
    #                       arrow = Arrow(x,y,dx,dy, width=arrowwidth)
    #                       if xmax<x:
    #                               xmax=x
    #                       if ymax<y:
    #                               ymax=y
    #                       acol = flux_log10[i]
    #                       apatches.append(arrow)
    #                       acolor.append(acol)
    #       elif plot_type == 1:
    #               if x==I_am_the_target[0] and y==I_am_the_target[1] and flux_log10[i]>=vmin:
    #                       arrowwidth = flux_log10[i]*m+b
    #                       arrow = Arrow(x,y,dx,dy, width=arrowwidth)
    #                       if xmax<x:
    #                               xmax=x
    #                       if ymax<y:
    #                               ymax=y
    #                       acol = flux_log10[i]
    #                       apatches.append(arrow)
    #                       acolor.append(acol)
    #               if x+dx==I_am_the_target[0] and y+dy==I_am_the_target[1] and flux_log10[i]>=vmin:
    #                       arrowwidth = flux_log10[i]*m+b
    #                       arrow = Arrow(x,y,dx,dy, width=arrowwidth)
    #                       if xmax<x:
    #                               xmax=x
    #                       if ymax<y:
    #                               ymax=y
    #                       acol = flux_log10[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.)
        rect.set_zorder(2)
        ax.add_patch(rect)
        xy = x+dx-0.5,y+dy-0.5
        rect = Rectangle(xy,1,1,ec='k',fc='None',fill='False',lw=1.)
        rect.set_zorder(2)
        ax.add_patch(rect)


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

    # colorbar label
    if which_flux == None or which_flux == 0:
        cb.set_label('log$_{10}$(f)')
    elif which_flux ==1:
        cb.set_label('log$_{10}$(E)')
    elif which_flux ==2:
        cb.set_label('log$_{10}$(timescale)')


    # plot file name
    graphname = 'flow-chart.'+format

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

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

    # plot mass numbers
    if imlabel==1:
        for z in range(nzmax):
            for n in range(nnmax):
                a = z+n
                if nzycheck[n,z,iarr-2]==1:
                    ax.text(n,z,a,horizontalalignment='center',verticalalignment='center',fontsize='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(np.argwhere(nzycheck[:,magic,iarr-2]))[0]
                    xnmax = max(np.argwhere(nzycheck[:,magic,iarr-2]))[0]
                    line = ax.plot([xnmin,xnmax],[magic,magic],lw=3.,color='r',ls='-')
                except ValueError:
                    dummy=0
            if magic<=nnmax:
                try:
                    yzmin = min(np.argwhere(nzycheck[magic,:,iarr-2]))[0]
                    yzmax = max(np.argwhere(nzycheck[magic,:,iarr-2]))[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]:
        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')
    ax.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))
    elif which_flux == 2:
        min_flux_label="min timescale [s] = "+str('{0:.4f}'.format(min_flux))
    if print_max_flux_in_plot:
        if which_flux == None or which_flux < 2:
            ax.text(plotaxis[1]-1.8,plotaxis[2]+0.1,max_flux_label,fontsize=10.)
        elif which_flux == 2:
            ax.text(plotaxis[1]-1.8,plotaxis[2]+0.1,min_flux_label,fontsize=10.)



    fig.savefig(graphname)
    print(graphname,'is done')
    if which_flux == None or which_flux < 2:
        print(max_flux_label,'for reaction =',ind_max_flux+1)
    elif which_flux == 2:
        print(min_flux_label,'for reaction =',ind_min_flux+1)

    plt.show()