## bound-free ionization and recombination rate
---

according to H&M2014 P273 eq-9.43, photoionization rate R_ik can be calculated by

`R_ik = 4 * PI * HP * alpha(freq, lower) * J(freq) * dfreq / freq `

**then sum over all frequency**.
  - this is a elementwise multiplication
  - freq is the frequency mesh of continuum
  - dfreq is the difference of freq, can be obtained by TRAPEZ method, for example :
    - dfreq[i] = freq[i+1] - freq[i-1], for i = 1,2,3,....,n-2
    - dfreq[0] = 0.5 * (freq[1] - freq[0])
    - dfreq[n-1] = 0.5 * ( freq[n-1] - freq[n-2] )
  - alpha(freq, lower) is a function to calculate bound-free photoionization cross section, which is a function of frequency and lower level n.
  - J(freq) is the mean intensity of this continuum transition as a function of frequency

and similarly, recombination rate `R_ki` can be calculated by  H&M2014 P273 eq-9.48

`R_ki = 4 * PI * HP * alpha(freq, lower) * ( 2*HP*freq**3/CL**2 + J(freq) ) * dfreq / freq `

**then sum over all frequency**.
  - this `R_ki` includes spontaneous and stimulated recombination
  - this is also a elementwise multiplication

photoionization_cross_section(freq, lower) can be found in H&M2014 p187 eq-7.84 :

`alpha = 2.815E29 * Z**4 / ( lower**5 * freq**3 ) * g_bf(freq, lower)`
  - in case of hydrogen, Z=1
  - lower is the lower level n
  - g_bf(freq, lower) calculates gaunt factors for bound-free hydrogen

```
def gaunt_factor(frequency, K):
    """
    gaunt factors for bound-free hydrogen (Gingerich, March 1964)
    Input:
            K: level for gaunt factor (integer)
    """
    W = 2.9979E15 /frequency

    if   K==1:
        return 0.9916 + ( (9.068E-3) - 0.2524*W )*W
    elif K==2:
        return 1.105 + ( (-7.922E-2) + (4.536E-3)*W )*W
    elif K==3:
        return 1.101 + ( (-3.290E-2) + (1.152E-3)*W )*W
    elif K==4:
        return 1.101 + ( (-1.923E-2) + (5.110E-4)*W )*W
    elif K==5:
        return 1.102 + ( (-0.01304) + (2.638E-4)*W )*W
    elif K==6:
        return 1.0986 + ( (-0.00902) + (1.367E-4)*W )*W
    else:
        return 1.0
```