plot_program.py
10 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 21 14:58:11 2022
@author: frede
"""
import main_program as mp
import matplotlib.pyplot as plt
import numpy as np
'''================================ initialization of parameters ================================'''
filename = './converted/kp00_10000.txt'
star_radius = 5
t_gas = 1000
n_e = 1.6
n_c = 54
parsec = 6e-9
fc_pah = 0.1
ISRF = True
test = mp.HeatingGas(filename, star_radius, t_gas, n_e, n_c, parsec, fc_pah, ISRF)
test.parameters()
# ''' parameter '''
# calculus = 0
# distance_list = np.array([parsec])
# for j in range(1,1001,1):
# d = distance_list[0]
# distance_list = np.append(distance_list, j * d) #in parsec
# for i in distance_list:
# test = mp.HeatingGas(filename, star_radius, t_gas, n_e, n_c, i, fc_pah, ISRF)
# test.parameters()
# calculus += test.total_gas_heating*1e-7/(test.g_0 * test.heating_efficiency * 0.1)
calculus = test.total_gas_heating*1e-7/(test.g_0 * test.heating_efficiency * 0.1)
# print('value = ',calculus/len(distance_list))
print('value = ',calculus)
# ''' conversion from mb/C atom to cm²/C atom '''
# test.pah_cross_n = test.pah_cross_n*1e-18
# test.pah_cross_c = test.pah_cross_c*1e-18
# test.pah_cross_dc = test.pah_cross_dc*1e-18
# '''================================ plots ================================'''
# ''' intensity per wavelength '''
# plt.figure()
# plt.xlim([0,1500])
# plt.plot(test.wavelength,test.wavelength_intensity*1e-3,'k')
# plt.ylabel('intensity spectrum (W m$^{-2}$ sr$^{-1}$ nm$^{-1}$)')
# plt.xlabel('wavelength (nm)')
# ''' average of 3 sizes categories of pah of photoabsorption sigma
# and photoionization sigma_ion cross sections per C atom of pahs
# as functions of photon energy in eV'''
# ''' small n_c '''
# test_small_n_c_1 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 32, parsec, fc_pah, ISRF)
# test_small_n_c_2 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 36, parsec, fc_pah, ISRF)
# test_small_n_c_3 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 38, parsec, fc_pah, ISRF)
# test_small_n_c_1.parameters()
# test_small_n_c_2.parameters()
# test_small_n_c_3.parameters()
# ''' medium n_c '''
# test_medium_n_c_1 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 40, parsec, fc_pah, ISRF)
# test_medium_n_c_2 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 42, parsec, fc_pah, ISRF)
# test_medium_n_c_1.parameters()
# test_medium_n_c_2.parameters()
# ''' large n_c '''
# test_large_n_c_1 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 48, parsec, fc_pah, ISRF)
# test_large_n_c_2 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 54, parsec, fc_pah, ISRF)
# test_large_n_c_3 = mp.HeatingGas(filename, star_radius, t_gas, n_e, 66, parsec, fc_pah, ISRF)
# test_large_n_c_1.parameters()
# test_large_n_c_2.parameters()
# test_large_n_c_3.parameters()
# ''' average of first and second photoionization yield '''
# yield_of_first_photoionization = (test_small_n_c_1.yield_of_first_photoionization +\
# test_small_n_c_2.yield_of_first_photoionization +\
# test_small_n_c_3.yield_of_first_photoionization +\
# test_medium_n_c_1.yield_of_first_photoionization +\
# 2*test_medium_n_c_2.yield_of_first_photoionization +\
# test_large_n_c_1.yield_of_first_photoionization +\
# test_large_n_c_2.yield_of_first_photoionization +\
# test_large_n_c_3.yield_of_first_photoionization )/9
# yield_of_second_photoionization = (test_small_n_c_1.yield_of_second_photoionization +\
# test_small_n_c_2.yield_of_second_photoionization +\
# test_small_n_c_3.yield_of_second_photoionization +\
# test_medium_n_c_1.yield_of_second_photoionization +\
# 2*test_medium_n_c_2.yield_of_second_photoionization +\
# test_large_n_c_1.yield_of_second_photoionization +\
# test_large_n_c_2.yield_of_second_photoionization +\
# test_large_n_c_3.yield_of_second_photoionization )/9
# ''' average of ionization cross-section for the cases of neutral molecules, and
# cations '''
# ionization_cross_n = test.pah_cross_n * yield_of_first_photoionization
# ionization_cross_c = test.pah_cross_c * yield_of_second_photoionization
# #ionization_cross_dc is too small to be traced
# ''' average of energy per cross-section for the cases of neutral molecules,
# cations and dications '''
# energy = (test.energy_neutral + test.energy_charged + test.energy_double_charged)/3
# ''' plot sigma and sigma_ion in cm²/(C atom) as a function of photon energy in eV
# for the cases of neutral, charged and double charged molecules '''
# plt.figure(figsize=(3.97,10.17))
# plt.subplot(3,1,1)
# plt.plot(test.energy_neutral,test.pah_cross_n,'k', label='$\sigma(E,Z)$')
# plt.plot(energy,ionization_cross_n,'k',ls='--', label='$\sigma_{ion}(E,Z)$')
# plt.title('Z=0', fontsize = 10)
# plt.legend(fontsize = 9)
# plt.subplot(3,1,2)
# plt.plot(test.energy_neutral,test.pah_cross_c,'k')
# plt.plot(energy,ionization_cross_c,'k',ls='--')
# plt.title('Z=1', fontsize = 10)
# plt.ylabel('Cross sections (cm$^2$/C)')
# plt.subplot(3,1,3)
# plt.plot(test.energy_neutral,test.pah_cross_dc,'k')
# plt.title('Z=2', fontsize = 10)
# plt.xlabel('Photon energy (eV)')
# plt.subplots_adjust(hspace = 0.25)
# if ISRF == False:
# ''' evolution of the different populations of neutral pah molecules,
# cations and dications as a function of the distance to the star considered '''
# ''' parameter '''
# population_frac = np.zeros([1, len(distance_list)])
# ''' indentation '''
# j=0
# gamma_list = []
# for i in distance_list:
# population_plot_test = mp.HeatingGas(filename, star_radius, t_gas, n_e, n_c, i, fc_pah, ISRF)
# population_plot_test.parameters()
# population_frac[0,j] =(population_plot_test.frac_charged +\
# population_plot_test.frac_double_charged)/\
# (population_plot_test.frac_neutral +\
# population_plot_test.frac_charged +\
# population_plot_test.frac_double_charged)
# # population_frac[0,j] = population_plot_test.frac_neutral
# # population_frac[1,j] = population_plot_test.frac_charged
# # population_frac[2,j] = population_plot_test.frac_double_charged
# gamma_list.append(population_plot_test.gamma)
# j=j+1
# ''' population fractions as function of the distance '''
# plt.figure(figsize=(5.90,3.93))
# plt.plot(gamma_list,population_frac[0,:],'k')
# plt.xlim([1e1, 1e7])
# # plt.plot(gamma_list, population_frac[0,:], color = 'green', label='Z=0')
# # plt.plot(gamma_list, population_frac[1,:], color = 'blue', label='Z=1')
# # plt.plot(gamma_list, population_frac[2,:], color = 'red', label='Z=2')
# plt.xlabel('$\gamma$ ($G_0\sqrt{T}/n_e$)')
# plt.ylabel('population fraction R$_i$')
# plt.legend()
# ''' photoelectric efficiencies with the pah heating model '''
# gamma_list_T_one = []
# heating_efficiency_list_T_one = []
# gamma_list_T_two = []
# heating_efficiency_list_T_two = []
# # for i in distance_list:
# # test_for_efficiencies_T_one = mp.HeatingGas(filename, star_radius, 100, n_e, 54, i, fc_pah, ISRF)
# # test_for_efficiencies_T_one.parameters()
# # gamma_list_T_one.append(test_for_efficiencies_T_one.gamma)
# # heating_efficiency_list_T_one.append(test_for_efficiencies_T_one.heating_efficiency)
# # test_for_efficiencies_T_two = mp.HeatingGas(filename, star_radius, 1000, n_e, 54, i, fc_pah, ISRF)
# # test_for_efficiencies_T_two.parameters()
# # gamma_list_T_two.append(test_for_efficiencies_T_two.gamma)
# # heating_efficiency_list_T_two.append(test_for_efficiencies_T_two.heating_efficiency)
# # plt.figure()
# # plt.xlim([1e1, 1e7])
# # plt.loglog(gamma_list_T_one, heating_efficiency_list_T_one, 'darkblue', label='$\epsilon_{PAH}$ (T = 100K)')
# # plt.loglog(gamma_list_T_two, heating_efficiency_list_T_two, 'darkred', ls='--' , label='$\epsilon_{PAH}$ (T = 1000K)')
# # plt.xlabel('$\gamma(G_0\sqrt{T}/n_e)$')
# # plt.ylabel('$\epsilon_{PAH}$')
# # plt.legend()
# # ''' Total photoelectric heating rates of the gas '''
# # gamma_list_fc_one = []
# # total_gas_heating_per_rf_list_fc_one = []
# # gamma_list_fc_two = []
# # total_gas_heating_per_rf_list_fc_two = []
# # for i in distance_list:
# # test_for_total_heating_rates_fc_one = mp.HeatingGas(filename, star_radius, t_gas, n_e, n_c, i, 0.05, ISRF)
# # test_for_total_heating_rates_fc_one.parameters()
# # gamma_list_fc_one.append(test_for_total_heating_rates_fc_one.gamma)
# # total_gas_heating_per_rf_list_fc_one.append( (test_for_total_heating_rates_fc_one.total_gas_heating*1e-7) / test_for_total_heating_rates_fc_one.g_0)
# # test_for_total_heating_rates_fc_two = mp.HeatingGas(filename, star_radius, t_gas, n_e, n_c, i, 0.1, ISRF)
# # test_for_total_heating_rates_fc_two.parameters()
# # gamma_list_fc_two.append(test_for_total_heating_rates_fc_two.gamma)
# # total_gas_heating_per_rf_list_fc_two.append( (test_for_total_heating_rates_fc_two.total_gas_heating*1e-7) / test_for_total_heating_rates_fc_two.g_0)
# # plt.figure()
# # plt.xlim([1e1, 1e7])
# # plt.loglog(gamma_list_fc_one, total_gas_heating_per_rf_list_fc_one, 'red', label='This work (5% of C in PAHs)')
# # plt.loglog(gamma_list_fc_two, total_gas_heating_per_rf_list_fc_two, 'blue', label='This work (10% of C in PAHs)')
# # plt.xlabel('$\gamma(G_0\sqrt{T}/n_e)$')
# # plt.ylabel('$\Gamma_{tot}$ (W H$^{-1}$ $G_0^{-1}$)')
# # plt.legend()