docs: def average_iso_abund_marco(self,mass_range,cycle,stable,i_decay):
          """
          Interface to average over mass_range.
  
          Parameters
          ----------
          mass_range : list
              A 1x2 array required to plot data in a certain mass range.
              Needed for _read_iso_abund_marco.
          cycle : integer
              which cycle from the h5 file?.  Needed for _read_iso_abund_marco
          stable : boolean
              Do you want to plot only stable or not.
          i_decay : integer
              If i_decay is 1, then plot not decayed.  If i_decay is 2,
              then plot decayed.  Make sense only if stable is true.
  
          See Also
          --------
          _read_iso_abund_marco()
  
          """
  
  
          import nuutils as u
  
          if not stable and i_decay == 2:
              print('ERROR: choose i_decay = 1')
              return
  
          #data=mp.se(directory,name_h5_file)
          self._read_iso_abund_marco(mass_range,cycle)
          #print spe
          if i_decay == 2:
              u.stable_specie()
              self.decay(self.mass_frac)
  
          # here I am calculating average mass fraction for all isotopes in given mass range, and then
          # if needed calculating average over decayed.
          # warning: mass_range is bigger than used_masses range, by definition. Should I use it?
          print('average over used_masses range, not over original mass_range')
          print(used_masses[0],used_masses[len(used_masses)-1],'instead of',mass_range[0],mass_range[1])
  
          global average_mass_frac
          average_mass_frac = []
  
          if len(used_masses) >= 2:
              dm_tot = abs(used_masses[len(used_masses)-1]-used_masses[0])
              for j in range(len(u.spe)-1):
                  temp = 0.
                  for i in range(len(used_masses)-1):
                      dm_i = abs(used_masses[i+1]-used_masses[i])
                      temp = float(self.mass_frac[i][j]*dm_i/dm_tot) + temp
                  average_mass_frac.append(temp)
              #print average_mass_frac
          elif  len(used_masses) == 1:
              print('case with 1 mass zone only, not implemented yet')
  
  
  
          somma = 0.
          somma = sum(average_mass_frac)
          print('departure from 1 of sum of average_mass_frac=',abs(1. - somma))
  
          # not let's do it over decayed also, if i_decay = 2
          if i_decay == 2:
              global average_mass_frac_decay
              average_mass_frac_decay = []
              dm_tot = abs(used_masses[len(used_masses)-1]-used_masses[0])
              #
              #print len(decayed_multi_d[0]),decayed_multi_d[0]
              for j in range(len(u.back_ind)):
                  temp = 0.
                  for i in range(len(used_masses)-1):
                      dm_i = abs(used_masses[i+1]-used_masses[i])
                      temp = float(decayed_multi_d[i][j]*dm_i/dm_tot) + temp
                  average_mass_frac_decay.append(temp)
  
              somma = 0.
              somma = sum(average_mass_frac_decay)
              print('departure from 1 of sum of average_mass_frac_decay=',abs(1. - somma))