docs: def solar(filename_solar, solar_factor):
      '''
      read solar abundances from filename_solar.
  
      Parameters
      ----------
      filename_solar : string
          The file name.
      solar_factor : float
          The correction factor to apply, in case filename_solar is not
          solar, but some file used to get initial abundances at
          metallicity lower than solar. However, notice that this is
          really rude, since alpha-enahncements and things like that are
          not properly considered.  Only H and He4 are not multiplied. So,
          for publications PLEASE use proper filename_solar at...solar,
          and use solar_factor = 1. Marco
  
      '''
  
      f0=open(filename_solar)
      sol=f0.readlines()
      f0.close
      sol[0].split("         ")
  
      # Now read in the whole file and create a hashed array:
      global names_sol
      names_sol=[]
      global z_sol
      z_sol=[]
      yps=np.zeros(len(sol))
      mass_number=np.zeros(len(sol))
      for i in range(len(sol)):
          z_sol.append(int(sol[i][1:3]))
          names_sol.extend([sol[i].split("         ")[0][4:]])
          yps[i]=float(sol[i].split("         ")[1]) * solar_factor
          try:
              mass_number[i]=int(names_sol[i][2:5])
          except ValueError:
              print("WARNING:")
              print("This initial abundance file uses an element name that does")
              print("not contain the mass number in the 3rd to 5th position.")
              print("It is assumed that this is the proton and we will change")
              print("the name to 'h   1' to be consistent with the notation used in")
              print("iniab.dat files")
              names_sol[i]='h   1'
              mass_number[i]=int(names_sol[i][2:5])
          if mass_number[i] == 1 or mass_number[i] == 4:
              yps[i] = old_div(yps[i],solar_factor)
      #  convert 'h   1' in prot, not needed any more??
      #names_sol[0] = 'prot '
  
  
      # now zip them together:
      global solar_abundance
      solar_abundance={}
      for a,b in zip(names_sol,yps):
          solar_abundance[a] = b
  
  
  
      z_bismuth = 83
      global solar_elem_abund
      solar_elem_abund = np.zeros(z_bismuth)
  
  
      for i in range(z_bismuth):
          dummy = 0.
          for j in range(len(solar_abundance)):
              if z_sol[j] == i+1:
                  dummy = dummy + float(solar_abundance[names_sol[j]])
          solar_elem_abund[i] = dummy