Various small corrections suggested by L. Klijn

[ladsantos]

Lars Klijn kindly sent a list of suggestions to apply in the code, which I list below from their email:

p-winds code:

parker:

• Line 12: tools imported but not used.
• Lines 213 and 330: doesn’t matter much, but it is a small inconsistency: comments say that ‘x0’ is set at 2 for r>1, but it is actually at 2.1. The line “v_init = (np.array(r > 1, dtype=int) 2 + 0.1)” could be changed to “v_init = (np.array(r > 1, dtype=int) 1.9 + 0.1)” to be consistent with comments.

Hydrogen:

• Line 40: “Number density profile for the atmosphere, in units of 1 / cm 3.” Should be “Density profile for the atmosphere, in units of g / cm 3.”
• Line 80: comment says “Proton mass in unit of kg” but should be “Proton mass in unit of g”.
• Exact_phi is missing in the documentation of the ion_fraction function.

Helium/Oxygen/Carbon:

• Lines 600 and 909 in Helium, line 345 in Oxygen and line 433 in Carbon: “dr = np.concatenate((dr, np.array([r[-1], ])))” should be “dr = np.concatenate((dr, np.array([dr[-1], ])))” as it is in the hydrogen module, where lines 399 and 400 read:

  dr_norm = np.diff(r_norm)
  dr_norm = np.concatenate((dr_norm, np.array([dr_norm[-1], ])))

Carbon:

• Line 404: 7 values are returned, but only 4 of them are normalized. These 3 lines should be added:

  phi_cii = phi_cii / phi_unit
  a_cii = a_cii / (rs * 7.1492E+09) ** 2
  a_h_cii = a_h_cii / (rs * 7.1492E+09) ** 2

• Lines 517/518: I think both term21 and term22 defined below these lines should also be included in dfcii_dr, with a minus sign. So probably move lines 517/518 below the current line 523 and have them read: dfcii_dr = (term11 + term12 + term13 + term14 - term15 - term16 + term17 + term18 + term19 – term21 – term22) / _v

Oxygen:

• Line 242: “carbon (11.26 eV, or 1101 Angstrom)” should be “oxygen (13.62 eV, or 910 Angstrom)”
• Lines 289-293 should be replaced by:

f_oii_r (``numpy.ndarray``):
        Fraction of singly-ionized oxygen in function of radius.

Energetics:

Line 325: “tau_he_plus = he_plus_cross he_column” should be “tau_he_plus = he_plus_cross he_plus_column”

p-winds notebooks:

Exospheric metals:

• Cell 5: line “n_he = (rho_array rhos he_fraction / (1 + 4 he_fraction) / m_h)” should be ““n_he = (rho_array rhos he_fraction / (h_fraction + 4 he_fraction) / m_h)”
• Cell 6: line “n_c = (rho_array rhos c_fraction / (1 + 6 c_fraction) / m_h)” should have (h_fraction+4he_fraction) in the denominator instead of “1”.
• Cell 7: same as for cell 6, concerning line “n_o = (rho_array rhos o_fraction / (1 + 8 * o_fraction) / m_h)”.

Quickstart example:

• Cell 10/14 (and also cell 8 of Roche lobe effects): same as exospheric metals cell 5, line “n_he = (rho_array rhos he_fraction / (1 + 4 * he_fraction) / m_h)” should include h_fraction.
• After cell 10: “For now we will assume an impact factor of zero (so the planet is in the dead center of the star)”. This is an inconsistency, because the actual impact parameter is used in cell 11.

[LarsKlijn]

Thank you for taking a look at my email! I would like to add one more to the list that I discovered since I wrote you: In the helium, oxygen and carbon code you incorporate a safeguard against numerical instabilities in the calculation of the population fractions, but this is absent from the hydrogen code. I discovered this because one of my solutions had the ionization fraction oscillate around 1. So probably add the lines f_r[f_r < 0] = 1E-15 f_r[f_r > 1.0] = 1.0 after first finding f_r and in the relaxation loop. If I do this I get this plot:

[LarsKlijn]

Today I noticed two more things while looking through the oxygen and carbon codes again:

• The central wavelengths of the carbon and oxygen lines hard coded in the lines module, and used to produce the transit plots in the exospheric metal notebook are quoted in the documentation of lines.py and in the x-axis labels of the notebook to be the wavelengths in air, while they are actually the wavelengths in vacuum, as the NIST database has vacuum lines for lambda<185nm
• For carbon (similarly for oxygen): in the notebook you have this line: m_C = 6 * 1.67262192369e-27 # Carbon atomic mass in kg and in carbon.py there are these lines: k1 = h_fraction / (h_fraction + 4 * he_fraction + 6 * c_fraction) / m_h k2 = he_fraction / (h_fraction + 4 * he_fraction + 6 * c_fraction) / m_h k3 = c_fraction / (h_fraction + 4 * he_fraction + 6 * c_fraction) / m_h where the number "6" is used. Now this is the atomic number, but shouldn't the mass number be used here? As it is for helium of course. So 12 for carbon and 16 for oxygen instead of 6 and 8 respectively.

[ladsantos]

Thanks, Lars. Everything has been corrected in commit 8aa035c. One item that did not require correction: the number densities are given in units of 1 / cm^{-3} (number of atoms per unit volume, and not mass per unit volume).

[LarsKlijn]

Hi Leonardo, I still think that line 40 in hydrogen.py should be “Density profile for the atmosphere, in units of g / cm 3.” instead of “Number density profile for the atmosphere, in units of 1 / cm 3.”. This density parameter is used in line 87 to find the number density: n_tot = density / mu / m_h And when the function radiative_processes_exact() is called later in the module, a density (g/cm3) is passed as the density parameter. So it is used correctly, just the documentation is incorrect.

[ladsantos]

Gotcha, this is now corrected in commit 7f797de.