Commit 2ba9e7abf0c590cc70c2cae04fcb784f11b90483

Authored by Jalabert
1 parent 033199c9
Exists in master

No longer useful

Showing 1 changed file with 0 additions and 343 deletions   Show diff stats
three_levels_model.py deleted
@@ -1,343 +0,0 @@ @@ -1,343 +0,0 @@
1 -# -*- coding: utf-8 -*-  
2 -"""  
3 -Created on Tue Jun 14 10:48:49 2022  
4 -  
5 -@author: frede  
6 -"""  
7 -import numpy as np  
8 -import radiation_fields as rf  
9 -  
10 -'''  
11 -  
12 -This programme seeks to calculate how a neutral gas subjected to a radiation field  
13 -is heated, involving PAHs in the ionization of the gas  
14 -  
15 -This program works with the script "radiation_fields.py"  
16 -which calculates the scale factor of the radiation field G0  
17 -  
18 -Example of execution of the program for the star of 10 solar radius HD200775,  
19 -for a distance to the star at 20pc, for a gas at 750K, an electron density  
20 -at 2.4 cm-3 and a size of PAHs à 54 C atoms :  
21 -  
22 -HeatingGas(filename='HD200775_RF.txt', star_radius=10, t_gas=750,  
23 - n_e=1.6, n_c=54, parsec=20)  
24 -  
25 -If i consider the Interstellar Radiation Field with  
26 -an approach described by Habing (1968), for a gas at 50K and a  
27 -fraction of cosmic carbon locked in PAHs at 5%:  
28 -  
29 -HeatingGas('habing1968.txt', 1, t_gas=100, n_e=1.6, n_c=54,  
30 - parsec=1, fc_pah=0.05, ISRF=True)  
31 -  
32 -While considering ISRF, the value of the parameters  
33 -parsec and star_radius are no longer important  
34 -  
35 -'''  
36 -  
37 -''' Constants '''  
38 -h = 6.62607015e-34 #Planck constant in J s  
39 -c = 299792458 #Light speed in m s-1  
40 -eps_0 = 8.85418782e-12 #epsilon_0, vacuum permitivity in F (Farad) m-1  
41 -z_0 = 0 #the charge state of the neutral molecule  
42 -z_1 = 1 #the charge state of the ionized 1 molecule  
43 -one_in_4_pi_eps_0 = 1/( 4*np.pi*(eps_0/1e9) )  
44 -  
45 -''' Conversions '''  
46 -ev = 1.602176634e-19 # 1 ev = 1.60218e-19 J and value of electron charge  
47 -erg = 1e-7 # 1 erg = 1e-7 J  
48 -ev_to_erg = ev/erg #1 eV in erg  
49 -mb = 1e-18 #1Mb = 1e-18cm2, Mb for Megabarn (unit used to express the cross sectional area of nuclei)  
50 -  
51 -''' Saving parameters '''  
52 -dust_heating_rate = np.zeros([3])  
53 -ionization_rate = np.zeros([2])  
54 -recombination_rate = np.zeros([2])  
55 -gas_heating_rate = np.zeros([2])  
56 -intrinsic_efficiency = np.zeros([2])  
57 -  
58 -class HeatingGas:  
59 -  
60 - """  
61 -  
62 - ----------  
63 -  
64 - Returns total_gas_heating, g_0 , gamma , t_gas, n_e , n_c  
65 -  
66 - total_gas_heating : float,  
67 - gas heating rate  
68 - g_0 : float,  
69 - scaling factor of the radiation field  
70 - gamma : float,  
71 - ( g_0 * sqrt(t_gas) )/n_e ionization parameter  
72 - t_gas : float,  
73 - gas temperature  
74 - n_e : float,  
75 - electron density in cm-3  
76 - n_c : float,  
77 - number of carbon atoms in pah molecules (size of the pah)  
78 -  
79 - -------  
80 -  
81 - """  
82 -  
83 - def __init__(self, filename, star_radius, t_gas, n_e, n_c, parsec, fc_pah=0.1, ISRF=False):  
84 -  
85 - """  
86 -  
87 - ----------  
88 - filename : str,  
89 - name of the file containing wavelength in nm and  
90 - intensity in erg cm-2 s-1 sr-1 nm-1 of a star or of the interstellar medium  
91 - star_radius : float,  
92 - star radius, in unit of solar radius  
93 - n_e : float,  
94 - electron density in cm-3 (n_e = n_h * 1.6e-4, n_h : hydrogen density)  
95 - parsec : float,  
96 - distance in parsec from the star  
97 - fc_pah : float,  
98 - fraction of cosmic carbon locked in PAHs (default: 0.1)  
99 - ISRF : bool,  
100 - Interstellar Radiation Field; if true, then the filename is a file for ISRF  
101 - if false, the radiation field of a star is studied (default: False)  
102 - -------  
103 -  
104 - """  
105 -  
106 - ''' input parameters '''  
107 - self.filename = filename  
108 - self.t_gas = t_gas  
109 - self.n_e = n_e  
110 - self.n_c = n_c  
111 - self.parsec = parsec  
112 - self.star_radius = star_radius  
113 - self.fc_pah = fc_pah  
114 - self.ISRF = ISRF  
115 -  
116 - ''' parameters to be observed '''  
117 - self.g_0 = None  
118 - self.distance = None  
119 - self.wavelength = None  
120 - self.wavelength_intensity = None  
121 - self.energy_intensity = None  
122 - self.energy = None  
123 - self.energy_range = None  
124 - self.energy_neutral = None  
125 - self.energy_charged = None  
126 - self.energy_double_charged = None  
127 - self.pah_cross_n = None  
128 - self.pah_cross_c = None  
129 - self.pah_cross_dc = None  
130 - self.ip_neutral = None  
131 - self.ip_charged = None  
132 - self.yield_of_first_photoionization = None  
133 - self.yield_of_second_photoionization = None  
134 - self.heating_efficiency = None  
135 - self.total_gas_heating = None  
136 - self.frac_neutral = None  
137 - self.frac_charged = None  
138 - self.frac_double_charged = None  
139 -  
140 - def parameters(self):  
141 -  
142 - ''' others parameters to be returned '''  
143 - self.g_0, self.distance, self.wavelength, self.wavelength_intensity, self.energy_intensity, self.energy, self.ISRF, RF_list = rf.radiation_field(self.filename, self.star_radius, self.parsec, self.ISRF, RF_list=False)  
144 - self.gamma = ( self.g_0 * np.sqrt(self.t_gas) ) / self.n_e #ionization parameter  
145 -  
146 - ''' Ionization Potential (IP) estimation '''  
147 - a = (self.n_c/468)**(1/3) #molecule diameter, in nm  
148 - self.ip_neutral = 3.9 + one_in_4_pi_eps_0 * ( ( z_0 + (1/2) ) * (ev**2/a) + ( z_0 + 2 ) * (ev**2/a) *(0.03/a) ) * (1/ev)  
149 - #IP to ionize the molecule : neutral to charged, in ev  
150 - self.ip_charged = 3.9 + one_in_4_pi_eps_0 * ( ( z_1 + (1/2) ) * (ev**2/a) + ( z_1 + 2 ) * (ev**2/a) *(0.03/a) ) * (1/ev)  
151 - #IP to ionize the charged molecule : charged 1 to charged 2, in ev  
152 -  
153 - '''==========================|building of the cross section|======================='''  
154 - ''' derives a mean photoabsorption cross section of the molecule considered, in 3 size ranges'''  
155 -  
156 - ''' small size '''  
157 - self.energy_neutral,crossn_1_case1 = np.loadtxt('./neutrals/ovalene_neutral.txt',unpack=True) #C32  
158 - self.energy_charged,crossc_1_case1 = np.loadtxt('./cations/ovalene_cation.txt',unpack=True) #C32  
159 - self.energy_double_charged,crossdc_1_case1 = np.loadtxt('./dications/ovalene_dication.txt',unpack=True) #C32  
160 -  
161 - self.energy_neutral,crossn_2_case1 = np.loadtxt('./neutrals/tetrabenzocoronene_neutral.txt',unpack=True) #C36  
162 - self.energy_charged,crossc_2_case1 = np.loadtxt('./cations/tetrabenzocoronene_cation.txt',unpack=True) #C36  
163 - self.energy_double_charged,crossdc_2_case1 = np.loadtxt('./dications/tetrabenzocoronene_dication.txt',unpack=True) #C36  
164 -  
165 - self.energy_neutral,crossn_3_case1 = np.loadtxt('./neutrals/circumbiphenyl_neutral.txt',unpack=True) #C38  
166 - self.energy_charged,crossc_3_case1 = np.loadtxt('./cations/circumbiphenyl_cation.txt',unpack=True) #C38  
167 - self.energy_double_charged,crossdc_3_case1 = np.loadtxt('./dications/circumbiphenyl_dication.txt',unpack=True) #C38  
168 -  
169 - ''' medium size '''  
170 - self.energy_neutral,crossn_1_case2=np.loadtxt('./neutrals/circumanthracene_neutral.txt',unpack=True) #C40  
171 - self.energy_charged,crossc_1_case2=np.loadtxt('./cations/circumanthracene_cation.txt',unpack=True) #C40  
172 - self.energy_double_charged,crossdc_1_case2=np.loadtxt('./dications/circumanthracene_dication.txt',unpack=True) #C40  
173 -  
174 - self.energy_neutral,crossn_2_case2 = np.loadtxt('./neutrals/circumpyrene_neutral.txt',unpack=True) #C42  
175 - self.energy_charged,crossc_2_case2 = np.loadtxt('./cations/circumpyrene_cation.txt',unpack=True) #C42  
176 - self.energy_double_charged,crossdc_2_case2 = np.loadtxt('./dications/circumpyrene_dication.txt',unpack=True) #C42  
177 -  
178 - self.energy_neutral,crossn_3_case2 = np.loadtxt('./neutrals/hexabenzocoronene_neutral.txt',unpack=True) #C42  
179 - self.energy_charged,crossc_3_case2 = np.loadtxt('./cations/hexabenzocoronene_cation.txt',unpack=True) #C42  
180 - self.energy_double_charged,crossdc_3_case2 = np.loadtxt('./dications/hexabenzocoronene_dication.txt',unpack=True) #C42  
181 -  
182 - ''' large size '''  
183 - self.energy_neutral,crossn_1_case3 = np.loadtxt('./neutrals/dicoronylene_neutral.txt',unpack=True) #C48  
184 - self.energy_charged,crossc_1_case3 = np.loadtxt('./cations/dicoronylene_cation.txt',unpack=True) #C48  
185 - self.energy_double_charged,crossdc_1_case3 = np.loadtxt('./dications/dicoronylene_dication.txt',unpack=True) #C48  
186 -  
187 - self.energy_neutral,crossn_2_case3 = np.loadtxt('./neutrals/circumcoronene_neutral.txt',unpack=True) #C54  
188 - self.energy_charged,crossc_2_case3 = np.loadtxt('./cations/circumcoronene_cation.txt',unpack=True) #C54  
189 - self.energy_double_charged,crossdc_2_case3 = np.loadtxt('./dications/circumcoronene_dication.txt',unpack=True) #C54  
190 -  
191 - self.energy_neutral,crossn_3_case3 = np.loadtxt('./neutrals/circumovalene_neutral.txt',unpack=True) #C66  
192 - self.energy_charged,crossc_3_case3 = np.loadtxt('./cations/circumovalene_cation.txt',unpack=True) #C66  
193 - self.energy_double_charged,crossdc_3_case3 = np.loadtxt('./dications/circumovalene_dication.txt',unpack=True) #C66  
194 - #for each cross section for each state of the molecule, we have an associated energy  
195 -  
196 - self.energy_range = np.where(self.energy_neutral<13.6)[0]  
197 - #energy_neutral|charged|double_charged are the same  
198 -  
199 - self.pah_cross_n = ( ( (crossn_1_case1/32)+(crossn_2_case1/36)+(crossn_3_case1/38) +\  
200 - (crossn_1_case2/40)+(crossn_2_case2/42)+(crossn_3_case2/42) +\  
201 - (crossn_1_case3/48)+(crossn_2_case3/54)+(crossn_3_case3/66) )/9 ) * self.n_c  
202 - self.pah_cross_c = ( ( (crossc_1_case1/32)+(crossc_2_case1/36)+(crossc_3_case1/38) +\  
203 - (crossc_1_case2/40)+(crossc_2_case2/42)+(crossc_3_case2/42) +\  
204 - (crossc_1_case3/48)+(crossc_2_case3/54)+(crossc_3_case3/66) )/9 ) * self.n_c  
205 - self.pah_cross_dc = ( ((crossdc_1_case1/32)+(crossdc_2_case1/36)+(crossdc_3_case1/38) +\  
206 - (crossdc_1_case2/40)+(crossdc_2_case2/42)+(crossdc_3_case2/42) +\  
207 - (crossdc_1_case3/48)+(crossdc_2_case3/54)+(crossdc_3_case3/66) )/9 ) * self.n_c  
208 - #cross_n is the average cross section of a pah of all types of size  
209 -  
210 - ''' Ranges imposed '''  
211 - self.energy_neutral = self.energy_neutral[self.energy_range]  
212 - self.energy_charged = self.energy_charged[self.energy_range]  
213 - self.energy_double_charged = self.energy_double_charged[self.energy_range]  
214 - self.pah_cross_n = self.pah_cross_n[self.energy_range]  
215 - self.pah_cross_c = self.pah_cross_c[self.energy_range]  
216 - self.pah_cross_dc = self.pah_cross_dc[self.energy_range]  
217 -  
218 - ''' yield from neutral to the first photoionization '''  
219 - first_part = np.where(self.energy_neutral<self.ip_neutral)[0]  
220 -  
221 - second_part = np.where( (self.energy_neutral>=self.ip_neutral) & (self.energy_neutral<(self.ip_neutral+9.2) ) )[0]  
222 - third_part = np.where((self.energy_neutral>self.ip_neutral+9.2))[0]  
223 -  
224 - y_1 = np.zeros([len(first_part)])  
225 - y_2 = ( self.energy_neutral[second_part]-self.ip_neutral )/9.2  
226 - y_3 = np.full(len(third_part), 1)  
227 -  
228 - self.yield_of_first_photoionization = np.concatenate( (y_1,y_2,y_3) )  
229 -  
230 - ''' yield from the first photoionization to the second photoionization '''  
231 - alpha = 0.3 #teepness coefficient, see Wenzel et al. 2020  
232 - if (self.n_c >= 32) & (self.n_c < 50):  
233 - beta = 0.59 + 8.1e-3 * self.n_c  
234 - if self.n_c >= 50:  
235 - beta = 1  
236 -  
237 - first_part = np.where( self.energy_charged<self.ip_charged )[0]  
238 - second_part = np.where( ( self.energy_charged>=self.ip_charged ) & ( self.energy_charged<11.3 ) )[0]  
239 - third_part = np.where( ( self.energy_charged>=11.3 ) & ( self.energy_charged<12.9 ) )[0]  
240 - fourth_part = np.where( ( self.energy_charged>=12.9 ) & ( self.energy_charged<13.6 ) )[0]  
241 -  
242 - y_1 = np.zeros([len(first_part)])  
243 - y_2 = ( alpha/(11.3-self.ip_charged) ) * (self.energy_charged[second_part]-self.ip_charged)  
244 - y_3 = np.full(len(third_part), alpha)  
245 - y_4 = ( (beta-alpha)/2.1 ) * (self.energy_charged[fourth_part]-12.9) + alpha  
246 -  
247 - self.yield_of_second_photoionization = np.concatenate( (y_1,y_2,y_3,y_4) )  
248 -  
249 - ''' adaptating the cross section from 0 to 13.6eV'''  
250 - pah_crossn = np.interp(self.energy,self.energy_neutral,self.pah_cross_n)*mb #in cm2/Carbon (from Mb/C to cm2/C)  
251 - pah_crossc = np.interp(self.energy,self.energy_charged,self.pah_cross_c)*mb #in cm2/Carbon (from Mb/C to cm2/C)  
252 - pah_crossdc = np.interp(self.energy,self.energy_double_charged,self.pah_cross_dc)*mb #in cm2/Carbon (from Mb/C to cm2/C)  
253 -  
254 - yield_n = np.interp(self.energy,self.energy_neutral,self.yield_of_first_photoionization)  
255 - yield_c = np.interp(self.energy,self.energy_charged,self.yield_of_second_photoionization)  
256 - #interpolation  
257 -  
258 - '''===================== dust and gas heating calculation ==================='''  
259 - energy_range_power_absorbed = np.where(self.energy<=13.6)[0]  
260 - energy_range_neutral = np.where(np.logical_and(self.energy >= self.ip_neutral, self.energy <= 13.6))[0]  
261 - energy_range_charged = np.where(np.logical_and(self.energy >= self.ip_charged, self.energy <= 13.6))[0]  
262 -  
263 - ''' photoabsorption of the neutrals, cations and dications molecules '''  
264 - photo_absorption_n = self.energy_intensity*pah_crossn*2*np.pi #erg s-1 eV-1 /!\ the 2*np.pi is the solid angle considered => the RF comes from the star only  
265 - photo_absorption_c = self.energy_intensity*pah_crossc*2*np.pi #erg s-1 eV-1  
266 - photo_absorption_dc = self.energy_intensity*pah_crossdc*2*np.pi #erg s-1 eV-1  
267 -  
268 - ''' power density absorbed for ionization '''  
269 - ionization_absorption_n = yield_n*photo_absorption_n #erg s-1 eV-1  
270 - ionization_absorption_c = yield_c*photo_absorption_c #erg s-1 eV-1  
271 -  
272 - ''' number of ionizations '''  
273 - number_ionization_absorption_n = ionization_absorption_n/(self.energy*ev_to_erg)  
274 - number_ionization_absorption_c = ionization_absorption_c/(self.energy*ev_to_erg)  
275 - #number of ionizations per s per eV for a charge state  
276 -  
277 - ''' photoemission rate '''  
278 - kpe_neutral = np.trapz(number_ionization_absorption_n[energy_range_neutral], (self.energy-self.ip_neutral)[energy_range_neutral])  
279 - kpe_charged = np.trapz(number_ionization_absorption_c[energy_range_charged], (self.energy-self.ip_charged)[energy_range_charged])  
280 - #in s-1  
281 -  
282 - ''' recombination rate '''  
283 - phi = ( 1.85*1e5 )/( self.t_gas*np.sqrt(self.n_c) ) #dimensionless  
284 - krec_neutral = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi )  
285 - krec_charged = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi * (1 + z_1) )  
286 - #in s-1  
287 -  
288 - ''' population fraction computation '''  
289 - #cations  
290 - self.frac_charged = 1/( 1+( krec_neutral/kpe_neutral) + (kpe_charged/krec_charged) )  
291 - #neutrals  
292 - self.frac_neutral = (1 - self.frac_charged * (kpe_charged/krec_charged) )/( 1 + (kpe_neutral/krec_neutral) )  
293 - #dications  
294 - self.frac_double_charged = (1 - self.frac_neutral)/( 1 + (krec_charged/kpe_charged) )  
295 -  
296 - ''' selection of the partition coefficient '''  
297 - partition_coeff = 0.46 #PAH parameter, 0.46 + or - 0.06  
298 -  
299 - ''' spectrum of the gas heating per charge state '''  
300 - neutral_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_neutral) *\  
301 - number_ionization_absorption_n * ev_to_erg #erg s-1 eV-1  
302 - charged_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_charged) *\  
303 - number_ionization_absorption_c * ev_to_erg #erg s-1 eV-1  
304 - #partition_coeff * (E-IP) is the kinetic energy of the photoelectron following absorption of a UV photon of energy E  
305 -  
306 - ''' gas heating rate per charge state '''  
307 - neutral_gas_heating_rate = np.trapz(neutral_heating_rate_spectrum[energy_range_neutral],\  
308 - (self.energy-self.ip_neutral)[energy_range_neutral] ) #erg s-1 molecule-1  
309 - charged_gas_heating_rate = np.trapz(charged_heating_rate_spectrum[energy_range_charged],\  
310 - (self.energy-self.ip_charged)[energy_range_charged] ) #erg s-1 molecule-1  
311 - #powers injected in the gas by photoelectrons ejected from  
312 - #neutral and cationic PAHs  
313 -  
314 - '''======================== heating efficiencies ========================'''  
315 -  
316 - ''' total power injected into the gas via the photoelectrons from pahs '''  
317 - total_injected_power = self.frac_neutral * neutral_gas_heating_rate +\  
318 - self.frac_charged * charged_gas_heating_rate  
319 - #erg s-1 /(molecule of size n_c)  
320 -  
321 - ''' heating rate of the molecule itself '''  
322 - heating_pah_n = np.trapz(photo_absorption_n[energy_range_power_absorbed],\  
323 - self.energy[energy_range_power_absorbed]) #erg s-1  
324 - heating_pah_c = np.trapz(photo_absorption_c[energy_range_power_absorbed],\  
325 - self.energy[energy_range_power_absorbed]) #erg s-1  
326 - heating_pah_dc = np.trapz(photo_absorption_dc[energy_range_power_absorbed],\  
327 - self.energy[energy_range_power_absorbed]) #erg s-1  
328 - #power absorbed by each PAH charge state  
329 -  
330 - ''' total power of the radiation absorbed by PAHs '''  
331 - total_absorbed_radiation_power = self.frac_neutral * heating_pah_n +\  
332 - self.frac_charged * heating_pah_c +\  
333 - self.frac_double_charged * heating_pah_dc  
334 - #erg s-1 /(molecule of size n_c)  
335 -  
336 - ''' heating efficiency '''  
337 - self.heating_efficiency = total_injected_power/total_absorbed_radiation_power  
338 -  
339 - ''' gas heating '''  
340 - self.total_gas_heating = total_injected_power * (self.fc_pah/self.n_c) * 2.7e-4  
341 - #2.7e-4 : elemental abundance of C relative to H (Tielens 2021)  
342 -  
343 - return self.total_gas_heating, self.g_0 , self.gamma , self.t_gas, self.n_e, self.n_c  
344 \ No newline at end of file 0 \ No newline at end of file