Commit e3daaba9aa0a470a70206341f247d440a7e88c3d

Authored by Jalabert
1 parent 3ae41466
Exists in master

ignored

Showing 1 changed file with 0 additions and 286 deletions   Show diff stats
photoionization_lib.py deleted
@@ -1,286 +0,0 @@ @@ -1,286 +0,0 @@
1 -#!/usr/bin/env python3  
2 -# -*- coding: utf-8 -*-  
3 -"""  
4 -Created on Tue Jun 30 09:54:20 2020  
5 -Library of the photoionization model  
6 -  
7 -@author: sacha.foschino@hotmail.fr  
8 -"""  
9 -import matplotlib.pyplot as plt  
10 -import numpy as np  
11 -from astropy.table import Table  
12 -  
13 -def preparation_Fe(d,Av,wave_range,quiet):  
14 - '''Function adapting the radiation field (RF) to its use by the PE model  
15 - d: float, distance between the region modeled and the UV source, in cm, e.g. 4.413e+17cm for the 0.143pc separating HD200775 and the PDR front,  
16 - Av: float, Extinction in the V band, in mag,  
17 - wave_range: list of 2 elements, limits of the wavelength interval considered in nm, e.g. [91.16,2000],  
18 - quiet: bool, if False, plot the RF at each step of the conversion  
19 -  
20 - returns Fe_E_,E,Fe_l,wave ; i.e. respectively the RF in erg/sec/cm2/eV/sr, the photon energy vector,  
21 - the RF in erg/sec/cm2/nm/sr and the corresponding wavelength vector  
22 - '''  
23 - # =============================================================================  
24 - # preparation of HD200775 radiation field: dillution + extinction  
25 - '''radiation field on HD200775'''  
26 -# wave_rad,spec_rad=np.loadtxt(rf_name,unpack=True,skiprows=4) #if another RF than that of HD200775 is used (e.g. Draine, Habing or Mathis)  
27 - wave_rad, spec_rad=np.loadtxt('./radiation_fields/HD200775_RF.txt', unpack=True)  
28 -  
29 - if not(quiet):  
30 - plt.figure(1)  
31 - plt.loglog(wave_rad,spec_rad,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
32 - plt.title('virgin spectrum')  
33 - plt.xlabel('wavelength (nm)')  
34 - plt.ylabel('Specific intensity (erg/sec/cm$^2$/nm/sr)')  
35 - plt.subplots_adjust(top=0.986,bottom=0.111,left=0.126,right=0.991,hspace=0.2,wspace=0.2)  
36 - plt.legend()  
37 - plt.tight_layout()  
38 -  
39 -  
40 - '''wavelengths ranges for the G0 energies'''  
41 - #selection of the wavelength range  
42 - spec7023=spec_rad[np.where((wave_rad>wave_range[0]) & (wave_rad<wave_range[1]))[0]]  
43 - wave=wave_rad[np.where((wave_rad>wave_range[0]) & (wave_rad<wave_range[1]))[0]]  
44 -  
45 - #intensity of HD200775  
46 - if not(quiet):  
47 - plt.figure(2)  
48 - plt.loglog(wave,spec7023,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
49 - plt.title('virgin spectrum cut')  
50 - plt.xlabel('wavelength (nm)')  
51 - plt.ylabel('Specific intensity (erg/sec/cm$^2$/nm/sr)')  
52 - plt.subplots_adjust(top=0.986,bottom=0.111,left=0.126,right=0.991,hspace=0.2,wspace=0.2)  
53 - plt.legend()  
54 - plt.tight_layout()  
55 -  
56 -  
57 - '''interpolation of the extinction coefficients with HD200775's RF wavelength vector'''  
58 - wave_draine__, albe, cos, cext__ = np.loadtxt('./extinction_curves/draine_curve_55.txt', unpack=True, usecols=[0,1,2,3], skiprows=80)  
59 - Cext=np.interp(wave,wave_draine__[::-1]*1e3,cext__[::-1]) #/!\ reverse the order of the cext__ vector  
60 -  
61 - '''geometrical dillution'''  
62 - R=10*70000000000 #radius of the star in cm, Pilleri et al. 2012  
63 - # d=0.143*3.086e18  
64 - intensity_diluted=spec7023*(R/d)**2  
65 -  
66 - if not(quiet):  
67 - plt.figure(3)  
68 - plt.loglog(wave,intensity_diluted,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
69 - plt.title('diluted spectrum')  
70 - plt.xlabel('wavelength (nm)')  
71 - plt.ylabel('Specific intensity (erg/sec/cm$^2$/nm/sr)')  
72 - plt.subplots_adjust(top=0.986,bottom=0.111,left=0.126,right=0.991,hspace=0.2,wspace=0.2)  
73 - plt.legend()  
74 - plt.tight_layout()  
75 -  
76 - '''extinction'''  
77 - NH=2e21 #conlumn density in cm-2, a column density of 2e21cm-2 implies an Av of 1mag  
78 -# Av=0  
79 - taunu=NH*Av*Cext #opacity  
80 - Fe_l=intensity_diluted*np.exp(-taunu) #intensity diluted and extincted  
81 - if not(quiet):  
82 - plt.figure(4)  
83 - plt.loglog(wave,Fe_l,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
84 - plt.title('ext spectrum')  
85 - plt.xlabel('wavelength (nm)')  
86 - plt.ylabel('Specific intensity (erg/sec/cm$^2$/nm/sr)')  
87 - plt.legend()  
88 - plt.tight_layout()  
89 -  
90 -  
91 -  
92 - '''conversion in erg/sec/cm2/sr/eV'''  
93 - c=299792458000000000 #light speed in nm  
94 - h=4.135667516e-15 #planck constant in eV*sec  
95 -  
96 - E=(h*c/wave)[::-1] #ev  
97 - Fe_E_=Fe_l[::-1]*h*c/(E**2)  
98 -  
99 - if not(quiet):  
100 - plt.figure(5)  
101 - plt.loglog(E,Fe_E_,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
102 -# plt.loglog(E,Fe_E_,label='d={:.0f}mpc, {:.0f}" '.format(d/3.086e15,d/(3.086e18*3.4268e-3)))  
103 - plt.title('energy spectrum')  
104 - plt.xlabel('photon energy (eV)')  
105 - plt.ylabel('Specific intensity (erg/sec/cm$^2$/eV/sr)')  
106 - plt.subplots_adjust(top=0.986,bottom=0.111,left=0.126,right=0.991,hspace=0.2,wspace=0.2)  
107 - plt.legend()  
108 - plt.tight_layout()  
109 -  
110 - return Fe_E_,E,Fe_l,wave  
111 -  
112 -  
113 -  
114 -  
115 -def yield_n_to_p(E,IP):  
116 - from numpy import where, concatenate,zeros,full  
117 - '''adapt the photoelectric yield Neutral -> Cation  
118 - version from Joachims et al 1996  
119 - E: 1xn array, photonenergy vector  
120 - IP: float, first ionization potential of the molecule  
121 -  
122 - returns the yield of the first photoionization  
123 - '''  
124 -  
125 - first_part=where(E<IP)[0]  
126 - scdn_part=where((E>=IP) & (E<(IP+9.2)))[0]  
127 - third_part=where((E>IP+9.2))[0]  
128 -  
129 - Y_1=zeros([len(first_part)])  
130 - Y_2=(E[scdn_part]-IP)/9.2  
131 - Y_3=full(len(third_part), 1)  
132 -  
133 - Yntop=concatenate((Y_1,Y_2,Y_3))  
134 - return Yntop  
135 -  
136 -  
137 -def yield_p_to_2p(E,IP2,Nc,alpha=0.3):  
138 - from numpy import where, concatenate,zeros,full  
139 - '''adapt the photoelectric yield Cation -> Dication to the molecule  
140 - version from Wenzel et al 2020  
141 - E: 1xn array, photonenergy vector  
142 - IP2: float, second ionization potential of the molecule  
143 - Nc: int or float, size of the PAH in number of carbon atoms  
144 - alpha: steepness coefficient, see Wenzel et al. 2020, this parameter should not change  
145 -  
146 - returns the yield of the second photoionization  
147 - '''  
148 -  
149 - if (Nc>=32) & (Nc<50):  
150 - beta=0.59+8.1e-3*Nc  
151 - if Nc>=50:  
152 - beta=1  
153 -  
154 - first_part=where(E<IP2)[0]  
155 - scdn_part=where((E>=IP2) & (E<(11.3)))[0]  
156 - third_part=where((E>=11.3) & (E<(12.9)))[0]  
157 - fourth_part=where((E>=12.9) & (E<(13.6)))[0]  
158 -# print(Nc)  
159 - Y_1=zeros([len(first_part)])  
160 - Y_2=alpha/(11.3-IP2)*(E[scdn_part]-IP2)  
161 - Y_3=full(len(third_part), alpha)  
162 - Y_4=(beta-alpha)/2.1*(E[fourth_part]-12.9)+alpha  
163 - Yptopp=concatenate((Y_1,Y_2,Y_3,Y_4))  
164 -  
165 - return Yptopp  
166 -  
167 -def IP_estimate(Nc,Z):  
168 - '''gives the estimate of the IP from Z to Z+1 according to Eq.1 of Wenzel et al. 2020  
169 - Nc: int of float, size of the PAH in carbon atoms  
170 - Z: positive int, the charge state of the molecule  
171 -  
172 - returns the IP of the molecule'''  
173 -  
174 - eps0=8.85e-21 #en F/nm  
175 - e=1.6e-19 #en C  
176 -  
177 - IP=3.9+(1/(4*np.pi*eps0))*((Z+0.5)*(e**2)/((Nc/468)**(1/3))+(Z+2)*((e**2)/((Nc/468)**(1/3)))*0.03/((Nc/468)**(1/3)))*1/e #en eV  
178 -  
179 - return IP  
180 -  
181 -  
182 -  
183 -  
184 -def cross_secs(Nc,quiet):  
185 - '''==========================|building of the cross section|======================='''  
186 - ''' derives a mean photoabsorption cross section of the molecule considered, in 3 size ranges  
187 - Nc: int, size of the molecule in carbon atoms  
188 - quiet: bool, if False plot the resulting cross sections'''  
189 - #absorption cross section  
190 - if (Nc>=32) & (Nc<40):  
191 - #32 a 40  
192 - eVN_,crossN_1_1=np.loadtxt('./neutrals/ovalene_neutral.txt',unpack=True) #C32  
193 - eVC_,crossC_1_1=np.loadtxt('./cations/ovalene_cation.txt',unpack=True) #C32  
194 - eVDC_,crossDC_1_1=np.loadtxt('./dications/ovalene_dication.txt',unpack=True) #C32  
195 -  
196 - eVN_,crossN_2_1=np.loadtxt('./neutrals/tetrabenzocoronene_neutral.txt',unpack=True) #C36  
197 - eVC_,crossC_2_1=np.loadtxt('./cations/tetrabenzocoronene_cation.txt',unpack=True) #C36  
198 - eVDC_,crossDC_2_1=np.loadtxt('./dications/tetrabenzocoronene_dication.txt',unpack=True) #C36  
199 -  
200 - eVN_,crossN_3_1=np.loadtxt('./neutrals/circumbiphenyl_neutral.txt',unpack=True) #C38  
201 - eVC_,crossC_3_1=np.loadtxt('./cations/circumbiphenyl_cation.txt',unpack=True) #C38  
202 - eVDC_,crossDC_3_1=np.loadtxt('./dications/circumbiphenyl_dication.txt',unpack=True) #C38  
203 -  
204 - cross_N_=(((crossN_1_1/32)+(crossN_2_1/36)+(crossN_3_1/38))/3)*Nc  
205 - cross_C_=(((crossC_1_1/32)+(crossC_2_1/36)+(crossC_3_1/38))/3)*Nc  
206 - cross_DC_=(((crossDC_1_1/32)+(crossDC_2_1/36)+(crossDC_3_1/38))/3)*Nc  
207 -  
208 - E_range=np.where(eVN_<13.6)[0]  
209 -  
210 - if not(quiet):  
211 - plt.figure(121)  
212 - plt.plot(eVN_[E_range],cross_N_[E_range],'b',label='Z=0')  
213 - plt.plot(eVC_[E_range],cross_C_[E_range],'r',label='Z=1')  
214 - plt.plot(eVDC_[E_range],cross_DC_[E_range],'orange',label='Z=2')  
215 - plt.ylabel('$\sigma$ (Mb)')  
216 - plt.xlabel('photon energy (eV)')  
217 - plt.legend()  
218 -# plt.plot(eVN_[E_range],crossN_1_1[E_range])  
219 -# plt.plot(eVN_[E_range],crossN_2_1[E_range])  
220 -# plt.plot(eVN_[E_range],crossN_3_1[E_range])  
221 - elif (Nc>=40) & (Nc<48):  
222 - #40 a 48  
223 -  
224 - eVN_,crossN_1_2=np.loadtxt('./neutrals/circumanthracene_neutral.txt',unpack=True) #C40  
225 - eVC_,crossC_1_2=np.loadtxt('./cations/circumanthracene_cation.txt',unpack=True) #C40  
226 - eVDC_,crossDC_1_2=np.loadtxt('./dications/circumanthracene_dication.txt',unpack=True) #C40  
227 -  
228 - eVN_,crossN_2_2=np.loadtxt('./neutrals/circumpyrene_neutral.txt',unpack=True) #C42  
229 - eVC_,crossC_2_2=np.loadtxt('./cations/circumpyrene_cation.txt',unpack=True) #C42  
230 - eVDC_,crossDC_2_2=np.loadtxt('./dications/circumpyrene_dication.txt',unpack=True) #C42  
231 -  
232 - eVN_,crossN_3_2=np.loadtxt('./neutrals/hexabenzocoronene_neutral.txt',unpack=True) #C42  
233 - eVC_,crossC_3_2=np.loadtxt('./cations/hexabenzocoronene_cation.txt',unpack=True) #C42  
234 - eVDC_,crossDC_3_2=np.loadtxt('./dications/hexabenzocoronene_dication.txt',unpack=True) #C42  
235 -  
236 - E_range=np.where(eVN_<13.6)[0]  
237 -  
238 - cross_N_=(((crossN_1_2/40)+(crossN_2_2/42)+(crossN_3_2/42))/3)*Nc  
239 - cross_C_=(((crossC_1_2/40)+(crossC_2_2/42)+(crossC_3_2/42))/3)*Nc  
240 - cross_DC_=(((crossDC_1_2/40)+(crossDC_2_2/42)+(crossDC_3_2/42))/3)*Nc  
241 - if not(quiet):  
242 -  
243 - plt.figure(121)  
244 - plt.plot(eVN_[E_range],cross_N_[E_range],'b',label='Z=0')  
245 - plt.plot(eVC_[E_range],cross_C_[E_range],'r',label='Z=1')  
246 - plt.plot(eVDC_[E_range],cross_DC_[E_range],'orange',label='Z=2')  
247 - plt.ylabel('$\sigma$ (Mb)')  
248 - plt.xlabel('photon energy (eV)')  
249 - plt.legend()  
250 -# plt.plot(eVN_[E_range],crossN_1_1[E_range])  
251 -# plt.plot(eVN_[E_range],crossN_2_1[E_range])  
252 -# plt.plot(eVN_[E_range],crossN_3_1[E_range])  
253 - elif (Nc>=48):  
254 - #48 a 66  
255 - eVN_,crossN_1_3=np.loadtxt('./neutrals/dicoronylene_neutral.txt',unpack=True) #C48  
256 - eVC_,crossC_1_3=np.loadtxt('./cations/dicoronylene_cation.txt',unpack=True) #C48  
257 - eVDC_,crossDC_1_3=np.loadtxt('./dications/dicoronylene_dication.txt',unpack=True) #C48  
258 -  
259 - eVN_,crossN_2_3=np.loadtxt('./neutrals/circumcoronene_neutral.txt',unpack=True) #C54  
260 - eVC_,crossC_2_3=np.loadtxt('./cations/circumcoronene_cation.txt',unpack=True) #C54  
261 - eVDC_,crossDC_2_3=np.loadtxt('./dications/circumcoronene_dication.txt',unpack=True) #C54  
262 -  
263 - eVN_,crossN_3_3=np.loadtxt('./neutrals/circumovalene_neutral.txt',unpack=True) #C66  
264 - eVC_,crossC_3_3=np.loadtxt('./cations/circumovalene_cation.txt',unpack=True) #C66  
265 - eVDC_,crossDC_3_3=np.loadtxt('./dications/circumovalene_dication.txt',unpack=True) #C66  
266 -  
267 - E_range=np.where(eVN_<13.6)[0]  
268 -  
269 - cross_N_=(((crossN_1_3/48)+(crossN_2_3/54)+(crossN_3_3/66))/3)*Nc  
270 - cross_C_=(((crossC_1_3/48)+(crossC_2_3/54)+(crossC_3_3/66))/3)*Nc  
271 - cross_DC_=(((crossDC_1_3/48)+(crossDC_2_3/54)+(crossDC_3_3/66))/3)*Nc  
272 - if not(quiet):  
273 - plt.figure(121)  
274 - plt.plot(eVN_[E_range],cross_N_[E_range],'b',label='Z=0')  
275 - plt.plot(eVC_[E_range],cross_C_[E_range],'r',label='Z=1')  
276 - plt.plot(eVDC_[E_range],cross_DC_[E_range],'orange',label='Z=2')  
277 - plt.ylabel('$\sigma$ (Mb)')  
278 - plt.xlabel('photon energy (eV)')  
279 - plt.legend()  
280 - return eVN_[E_range],eVC_[E_range],eVDC_[E_range],cross_N_[E_range],cross_C_[E_range],cross_DC_[E_range]  
281 -  
282 -  
283 -  
284 -  
285 -  
286 -  
287 \ No newline at end of file 0 \ No newline at end of file