docs: 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()