daily pastebin goal
15%
SHARE
TWEET

Untitled

a guest Aug 13th, 2017 49 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ## bound-free ionization and recombination rate
  2. ---
  3.  
  4. according to H&M2014 P273 eq-9.43, photoionization rate R_ik can be calculated by
  5.  
  6. `R_ik = 4 * PI * HP * alpha(freq, lower) * J(freq) * dfreq / freq `
  7.  
  8. **then sum over all frequency**.
  9.   - this is a elementwise multiplication
  10.   - freq is the frequency mesh of continuum
  11.   - dfreq is the difference of freq, can be obtained by TRAPEZ method, for example :
  12.     - dfreq[i] = freq[i+1] - freq[i-1], for i = 1,2,3,....,n-2
  13.     - dfreq[0] = 0.5 * (freq[1] - freq[0])
  14.     - dfreq[n-1] = 0.5 * ( freq[n-1] - freq[n-2] )
  15.   - alpha(freq, lower) is a function to calculate bound-free photoionization cross section, which is a function of frequency and lower level n.
  16.   - J(freq) is the mean intensity of this continuum transition as a function of frequency
  17.  
  18. and similarly, recombination rate `R_ki` can be calculated by  H&M2014 P273 eq-9.48
  19.  
  20. `R_ki = 4 * PI * HP * alpha(freq, lower) * ( 2*HP*freq**3/CL**2 + J(freq) ) * dfreq / freq `
  21.  
  22. **then sum over all frequency**.
  23.   - this `R_ki` includes spontaneous and stimulated recombination
  24.   - this is also a elementwise multiplication
  25.  
  26. photoionization_cross_section(freq, lower) can be found in H&M2014 p187 eq-7.84 :
  27.  
  28. `alpha = 2.815E29 * Z**4 / ( lower**5 * freq**3 ) * g_bf(freq, lower)`
  29.   - in case of hydrogen, Z=1
  30.   - lower is the lower level n
  31.   - g_bf(freq, lower) calculates gaunt factors for bound-free hydrogen
  32.  
  33. ```
  34. def gaunt_factor(frequency, K):
  35.     """
  36.     gaunt factors for bound-free hydrogen (Gingerich, March 1964)
  37.     Input:
  38.             K: level for gaunt factor (integer)
  39.     """
  40.     W = 2.9979E15 /frequency
  41.    
  42.     if   K==1:
  43.         return 0.9916 + ( (9.068E-3) - 0.2524*W )*W
  44.     elif K==2:
  45.         return 1.105 + ( (-7.922E-2) + (4.536E-3)*W )*W
  46.     elif K==3:
  47.         return 1.101 + ( (-3.290E-2) + (1.152E-3)*W )*W
  48.     elif K==4:
  49.         return 1.101 + ( (-1.923E-2) + (5.110E-4)*W )*W
  50.     elif K==5:
  51.         return 1.102 + ( (-0.01304) + (2.638E-4)*W )*W
  52.     elif K==6:
  53.         return 1.0986 + ( (-0.00902) + (1.367E-4)*W )*W
  54.     else:
  55.         return 1.0
  56. ```
RAW Paste Data
Top