#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Tue Jul 7 15:00:58 2020 @author: sylviaploeckinger """ import numpy as np # constants m_e = 9.1094e-28 # electron rest mass [g] h = 6.626e-27 # planck constant [erg s] kB_eV = 8.617e-5 # boltzmann constant [eV K^-1] kB_erg = 1.380649e-16 # boltzmann constant [erg K^-1] E_hyd = 13.605 # eV # T [K], n [cm-3] def helper(n, T): return 1./n * np.power((2. * np.pi * m_e * kB_erg * T)/(h*h), 3./2.) * \ np.exp(-E_hyd / (kB_eV * T)) def boltzmann_fraction(n_low, n_high, T): g_high = np.float(n_high*n_high) g_low = np.float(n_low *n_low) E_high = - E_hyd / np.float(n_high * n_high) E_low = - E_hyd / np.float(n_low * n_low ) return g_high/g_low * np.exp(- (E_high-E_low) / (kB_eV * T)) logT = np.linspace(3, 4.5, 1000) T = np.power(10., logT) n = 1.e14 # typical density in a stellar atmosphere [cm-3] x1 = - helper(n,T)/2. + np.sqrt( np.power(helper(n,T)/2., 2) + helper(n,T) ) import matplotlib.pyplot as plt fig = plt.figure() #plt.gca().invert_xaxis() fig.subplots_adjust(left = 0.14) ax = plt.subplot(111) ax.set_title ('Hydrogen') ax.set_yscale('log') ax.set_ylim(1.e-12, 10) ax.set_xlabel('Temperature [K]') ax.set_ylabel('N$_{\mathrm{a}}$ / N$_{\mathrm{b}}$') ax.plot(T, x1, label = 'N$_{\mathrm{HII}}$ / (N$_{\mathrm{HI}}$ + N$_{\mathrm{HII}}$)',\ color = 'black', linestyle = 'solid' ) ax.plot(T, 1. - x1, label = 'N$_{\mathrm{HI}}$ / (N$_{\mathrm{HI}}$ + N$_{\mathrm{HII}}$)',\ color = 'black', linestyle = 'dashed') ax.plot (T, boltzmann_fraction(1, 2, T), \ label = 'N$_{\mathrm{n}=2}$ / N$_{\mathrm{n}=1}$', \ color = 'grey') ax.plot(T, (1. - x1) * boltzmann_fraction(1, 2, T), color = 'firebrick', \ label='[N$_{\mathrm{HI}}$ / (N$_{\mathrm{HI}}$ + N$_{\mathrm{HII}}$)] [N$_{\mathrm{n}=2}$ / N$_{\mathrm{n}=1}$]') ax.legend(loc = 'lower right') plt.close() fig.savefig('saha_hydrogen.png', dpi = 150)