old_contents: from scipy.interpolate import RectBivariateSpline
  import xxdata_11
  
  recombination = xxdata_11.read_acd('acd96_c.dat')
  ionisation = xxdata_11.read_scd('scd96_c.dat')
  
  #recombination = xxdata_11.read_acd('acd89_ar.dat')
  #ionisation = xxdata_11.read_scd('scd89_ar.dat')
  
  nuclear_charge = recombination['charge']
  
  logalpha, logS = [], []
  for i in xrange(nuclear_charge):
      d = recombination
      sp = RectBivariateSpline(d['temperature'], d['density'], d['coeff_table'][i])
      logalpha.append(sp)
  
      d = ionisation
      sp = RectBivariateSpline(d['temperature'], d['density'], d['coeff_table'][i])
      logS.append(sp)
  
  
  import matplotlib.pyplot as plt
  import numpy as np
  
  logdens = 13 # log10(n_e[cm-3])
  
  temperature = np.logspace(0, 3, 300)
  logtemp = np.log10(temperature)
  
  y = np.zeros((nuclear_charge+1, len(temperature)))
  y[0] = np.ones_like(temperature)
  for i in xrange(nuclear_charge):
      alpha = logalpha[i](logtemp, logdens)
      alpha = np.power(10, alpha).squeeze()
  
      S = logS[i](logtemp, logdens)
      S = np.power(10, S).squeeze()
      y[i+1] = y[i] * S / alpha
  
  y /= y.sum(0) # fractional abundance
  
  plt.figure(1); plt.clf()
  ax = plt.gca()
  
  ax.loglog(temperature, y.T)
  
  ax.set_xlabel('$T_\mathrm{e}\ [\mathrm{eV}]$')
  ax.set_ylim(0.05, 1)
  
  max_pos = y.argmax(axis=-1)
  for i in xrange(nuclear_charge+1):
      index = max_pos[i]
      xy = temperature[index], y[i, index]
      s = '$%d$' % (i+1,)
      ax.annotate(s, xy, ha='center', va='top')
  
  
  plt.draw()
  plt.show()
  
new_contents: from scipy.interpolate import RectBivariateSpline
  import xxdata_11
  
  recombination = xxdata_11.read_acd('acd96_c.dat')
  ionisation = xxdata_11.read_scd('scd96_c.dat')
  
  #recombination = xxdata_11.read_acd('acd89_ar.dat')
  #ionisation = xxdata_11.read_scd('scd89_ar.dat')
  
  nuclear_charge = recombination['charge']
  
  logalpha, logS = [], []
  for i in xrange(nuclear_charge):
      d = recombination
      sp = RectBivariateSpline(d['temperature'], d['density'], d['coeff_table'][i])
      logalpha.append(sp)
  
      d = ionisation
      sp = RectBivariateSpline(d['temperature'], d['density'], d['coeff_table'][i])
      logS.append(sp)
  
  
  import matplotlib.pyplot as plt
  import numpy as np
  
  logdens = 13 # log10(n_e[cm-3])
  
  temperature = np.logspace(0, 3, 300)
  logtemp = np.log10(temperature)
  
  y = np.zeros((nuclear_charge+1, len(temperature)))
  y[0] = np.ones_like(temperature)
  for i in xrange(nuclear_charge):
      alpha = logalpha[i](logtemp, logdens)
      alpha = np.power(10, alpha).squeeze()
  
      S = logS[i](logtemp, logdens)
      S = np.power(10, S).squeeze()
      y[i+1] = y[i] * S / alpha
  
  y /= y.sum(0) # fractional abundance
  
  plt.figure(1); plt.clf()
  ax = plt.gca()
  
  ax.loglog(temperature, y.T)
  
  ax.set_xlabel('$T_\mathrm{e}\ [\mathrm{eV}]$')
  ax.set_ylim(0.05, 1.5)
  
  max_pos = y.argmax(axis=-1)
  for i in xrange(nuclear_charge+1):
      index = max_pos[i]
      xy = temperature[index], y[i, index]
      s = '$%d$' % (i+1,)
      ax.annotate(s, xy, ha='center')
  
  
  plt.draw()
  plt.show()
  
subject: Add some space to the top.