From d08e41f0a727058a5e2364b46011bf97729d49cd Mon Sep 17 00:00:00 2001 From: Jalabert Date: Fri, 22 Jul 2022 10:11:05 +0200 Subject: [PATCH] feat files --- .gitignore.txt | 3 ++- .idea/.gitignore | 8 ++++++++ README.md | 7 +++++++ README.rtf | 36 ------------------------------------ main_program.py | 343 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- three_levels_model.py | 343 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 360 insertions(+), 380 deletions(-) create mode 100644 .idea/.gitignore delete mode 100644 README.rtf delete mode 100644 main_program.py create mode 100644 three_levels_model.py diff --git a/.gitignore.txt b/.gitignore.txt index 441740a..4170e84 100644 --- a/.gitignore.txt +++ b/.gitignore.txt @@ -10,4 +10,5 @@ photoionization_lib.py plot_program.py reproj_lib.py test_plot.py -Programs.zip \ No newline at end of file +Programs.zip +fov_data_NGC7023_H2.py \ No newline at end of file diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..13566b8 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,8 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/README.md b/README.md index e69de29..bbfe525 100644 --- a/README.md +++ b/README.md @@ -0,0 +1,7 @@ +## Requirement + +pip install numpy +pip install astropy +pip install matplotlib + +## Programs \ No newline at end of file diff --git a/README.rtf b/README.rtf deleted file mode 100644 index f4cfc4b..0000000 --- a/README.rtf +++ /dev/null @@ -1,36 +0,0 @@ -{\rtf1\ansi\ansicpg1252\cocoartf1671\cocoasubrtf600 -{\fonttbl\f0\fswiss\fcharset0 Helvetica;} -{\colortbl;\red255\green255\blue255;} -{\*\expandedcolortbl;;} -\paperw11900\paperh16840\margl1440\margr1440\vieww14300\viewh7720\viewkind0 -\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\pardirnatural\partightenfactor0 - -\f0\fs24 \cf0 Welcome in the directory concerning the study of the influence of PAH on the gas heating by photoelectric effect.\ -\ -Description of the files:\ - - coupes_fits.py: library defining a function allowing to cut a map between 2 points. Used in the code \'ab\'a0photoelectric_heating_model_ngc7023nw.py\'a0\'bb to obtain the cuts along 3 lines of sight with HD200775 as a reference.\ - \ - -fov_data_NGC7023_H2.py : code extracting the useful information from all the maps stored in \'abthe \'ab\'a0data\'a0\'bb directory. Read its documentation to know how it produces a dictionary containing the information on the common field of view between IRS and PACS observations.\ -\ - - NGC7023_map_analysis_f0evsg.py : code in which some computation and plots are done, in particular to compute the PE heating efficiency, ionization fraction\'85\ -\ - - \'ab\'a0photoelectric_heating_model_ngc7023nw.py\'a0\'bb : Mol\'e9cular model of the photoelectric heating on PAH, using data of laboratory experiments and quantum chemistry computations. \ -All the directories here are related to this code, see below.\ -\ - - \'ab\'a0 photoionization_lib.py\'a0\'bb : library of functions used is the code \'ab\'a0photoelectric_heating_model_ngc7023nw.py\'a0\'bb.\ -\ - - reproj_lib.py : library of the functions used in the code \'abfov_data_NGC7023_H2.py\'a0\'bb \'a0to produce maps downgraded to the coordinate frame and spatial resolution of the [Cii] 158mu PACS map (11\'92\'92).\ -\ -Look at the \'ab\'a0reproj_lib.py\'a0\'bb description for the methodology used to produce the maps.\ - \ -\ -Description of the directories:\ - - \'ab\'a0neutrals\'a0\'bb, \'ab\'a0cations\'a0\'bb and \'ab\'a0dication\'a0\'bb: directories where the photo absorption cross sections are stored (Spectral database of Theoretical PAH, Malocci et al 2007)\ - - images_molecules : png image of the molecules considered (from the database)\ - -yields: now useless in the code of the model (the yields are generated by a function) but contains the photoionization yields of several studies (see titles and descriptions of these files)\ - -extinction_curves : Obvious purpose, from Weingartner and Draine 2001\ - -radiation_fields : contains the radiation fields of Habing (1968), Mathis (1983) and Draine (noted ISRF). They come from the following website: https://home.strw.leidenuniv.nl/~ewine/photo/radiation_fields.html\ -The HD200775_RF.txt is the radiation field of the star HD200775\ -\ - \ -} \ No newline at end of file diff --git a/main_program.py b/main_program.py deleted file mode 100644 index 20840ff..0000000 --- a/main_program.py +++ /dev/null @@ -1,343 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Jun 14 10:48:49 2022 - -@author: frede -""" -import numpy as np -import radiation_fields as rf - -''' - -This programme seeks to calculate how a neutral gas subjected to a radiation field -is heated, involving pahs in the ionization of the gas - -This program works with the script "radiation_fields.py" -which calculates the scale factor of the radiation field G0 - -Example of execution of the program for the star of 10 solar radius HD200775, -for a distance to the star at 20pc, for a gas at 750K, an electron density -at 2.4 cm-3 and a size of PAHs à 54 C atoms : - -HeatingGas(filename='HD200775_RF.txt', star_radius=10, t_gas=750, - n_e=1.6, n_c=54, parsec=20) - -If i consider the Interstellar Radiation Field with -an approach described by Habing (1968), for a gas at 50K and a -fraction of cosmic carbon locked in PAHs at 5%: - -HeatingGas('habing1968.txt', 1, t_gas=100, n_e=1.6, n_c=54, - parsec=1, fc_pah=0.05, ISRF=True) - -While considering ISRF, the value of the parameters -parsec and star_radius are no longer important - -''' - -''' Constants ''' -h = 6.62607015e-34 #Planck constant in J s -c = 299792458 #Light speed in m s-1 -eps_0 = 8.85418782e-12 #epsilon_0, vacuum permitivity in F (Farad) m-1 -z_0 = 0 #the charge state of the neutral molecule -z_1 = 1 #the charge state of the ionized 1 molecule -one_in_4_pi_eps_0 = 1/( 4*np.pi*(eps_0/1e9) ) - -''' Conversions ''' -ev = 1.602176634e-19 # 1 ev = 1.60218e-19 J and value of electron charge -erg = 1e-7 # 1 erg = 1e-7 J -ev_to_erg = ev/erg #1 eV in erg -mb = 1e-18 #1Mb = 1e-18cm2, Mb for Megabarn (unit used to express the cross sectional area of nuclei) - -''' Saving parameters ''' -dust_heating_rate = np.zeros([3]) -ionization_rate = np.zeros([2]) -recombination_rate = np.zeros([2]) -gas_heating_rate = np.zeros([2]) -intrinsic_efficiency = np.zeros([2]) - -class HeatingGas: - - """ - - ---------- - - Returns total_gas_heating, g_0 , gamma , t_gas, n_e , n_c - - total_gas_heating : float, - gas heating rate - g_0 : float, - scaling factor of the radiation field - gamma : float, - ( g_0 * sqrt(t_gas) )/n_e ionization parameter - t_gas : float, - gas temperature - n_e : float, - electron density in cm-3 - n_c : float, - number of carbon atoms in pah molecules (size of the pah) - - ------- - - """ - - def __init__(self, filename, star_radius, t_gas, n_e, n_c, parsec, fc_pah=0.1, ISRF=False): - - """ - - ---------- - filename : str, - name of the file containing wavelength in nm and - intensity in erg cm-2 s-1 sr-1 nm-1 of a star or of the interstellar medium - star_radius : float, - star radius, in unit of solar radius - n_e : float, - electron density in cm-3 (n_e = n_h * 1.6e-4, n_h : hydrogene density) - parsec : float, - distance in parsec from the star - fc_pah : float, - fraction of cosmic carbon locked in PAHs (default: 0.1) - ISRF : bool, - Interstellar Radiation Field; if true, then the filename is a file for ISRF - if false, the radiation field of a star is studied (default: False) - ------- - - """ - - ''' input parameters ''' - self.filename = filename - self.t_gas = t_gas - self.n_e = n_e - self.n_c = n_c - self.parsec = parsec - self.star_radius = star_radius - self.fc_pah = fc_pah - self.ISRF = ISRF - - ''' parameters to be observed ''' - self.g_0 = None - self.distance = None - self.wavelength = None - self.wavelength_intensity = None - self.energy_intensity = None - self.energy = None - self.energy_range = None - self.energy_neutral = None - self.energy_charged = None - self.energy_double_charged = None - self.pah_cross_n = None - self.pah_cross_c = None - self.pah_cross_dc = None - self.ip_neutral = None - self.ip_charged = None - self.yield_of_first_photoionization = None - self.yield_of_second_photoionization = None - self.heating_efficiency = None - self.total_gas_heating = None - self.frac_neutral = None - self.frac_charged = None - self.frac_double_charged = None - - def parameters(self): - - ''' others parameters to be returned ''' - 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) - self.gamma = ( self.g_0 * np.sqrt(self.t_gas) ) / self.n_e #ionization parameter - - ''' Ionization Potential (IP) estimation ''' - a = (self.n_c/468)**(1/3) #molecule diameter, in nm - 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) - #IP to ionize the molecule : neutral to charged, in ev - 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) - #IP to ionize the charged molecule : charged 1 to charged 2, in ev - - '''==========================|building of the cross section|=======================''' - ''' derives a mean photoabsorption cross section of the molecule considered, in 3 size ranges''' - - ''' small size ''' - self.energy_neutral,crossn_1_case1 = np.loadtxt('./neutrals/ovalene_neutral.txt',unpack=True) #C32 - self.energy_charged,crossc_1_case1 = np.loadtxt('./cations/ovalene_cation.txt',unpack=True) #C32 - self.energy_double_charged,crossdc_1_case1 = np.loadtxt('./dications/ovalene_dication.txt',unpack=True) #C32 - - self.energy_neutral,crossn_2_case1 = np.loadtxt('./neutrals/tetrabenzocoronene_neutral.txt',unpack=True) #C36 - self.energy_charged,crossc_2_case1 = np.loadtxt('./cations/tetrabenzocoronene_cation.txt',unpack=True) #C36 - self.energy_double_charged,crossdc_2_case1 = np.loadtxt('./dications/tetrabenzocoronene_dication.txt',unpack=True) #C36 - - self.energy_neutral,crossn_3_case1 = np.loadtxt('./neutrals/circumbiphenyl_neutral.txt',unpack=True) #C38 - self.energy_charged,crossc_3_case1 = np.loadtxt('./cations/circumbiphenyl_cation.txt',unpack=True) #C38 - self.energy_double_charged,crossdc_3_case1 = np.loadtxt('./dications/circumbiphenyl_dication.txt',unpack=True) #C38 - - ''' medium size ''' - self.energy_neutral,crossn_1_case2=np.loadtxt('./neutrals/circumanthracene_neutral.txt',unpack=True) #C40 - self.energy_charged,crossc_1_case2=np.loadtxt('./cations/circumanthracene_cation.txt',unpack=True) #C40 - self.energy_double_charged,crossdc_1_case2=np.loadtxt('./dications/circumanthracene_dication.txt',unpack=True) #C40 - - self.energy_neutral,crossn_2_case2 = np.loadtxt('./neutrals/circumpyrene_neutral.txt',unpack=True) #C42 - self.energy_charged,crossc_2_case2 = np.loadtxt('./cations/circumpyrene_cation.txt',unpack=True) #C42 - self.energy_double_charged,crossdc_2_case2 = np.loadtxt('./dications/circumpyrene_dication.txt',unpack=True) #C42 - - self.energy_neutral,crossn_3_case2 = np.loadtxt('./neutrals/hexabenzocoronene_neutral.txt',unpack=True) #C42 - self.energy_charged,crossc_3_case2 = np.loadtxt('./cations/hexabenzocoronene_cation.txt',unpack=True) #C42 - self.energy_double_charged,crossdc_3_case2 = np.loadtxt('./dications/hexabenzocoronene_dication.txt',unpack=True) #C42 - - ''' large size ''' - self.energy_neutral,crossn_1_case3 = np.loadtxt('./neutrals/dicoronylene_neutral.txt',unpack=True) #C48 - self.energy_charged,crossc_1_case3 = np.loadtxt('./cations/dicoronylene_cation.txt',unpack=True) #C48 - self.energy_double_charged,crossdc_1_case3 = np.loadtxt('./dications/dicoronylene_dication.txt',unpack=True) #C48 - - self.energy_neutral,crossn_2_case3 = np.loadtxt('./neutrals/circumcoronene_neutral.txt',unpack=True) #C54 - self.energy_charged,crossc_2_case3 = np.loadtxt('./cations/circumcoronene_cation.txt',unpack=True) #C54 - self.energy_double_charged,crossdc_2_case3 = np.loadtxt('./dications/circumcoronene_dication.txt',unpack=True) #C54 - - self.energy_neutral,crossn_3_case3 = np.loadtxt('./neutrals/circumovalene_neutral.txt',unpack=True) #C66 - self.energy_charged,crossc_3_case3 = np.loadtxt('./cations/circumovalene_cation.txt',unpack=True) #C66 - self.energy_double_charged,crossdc_3_case3 = np.loadtxt('./dications/circumovalene_dication.txt',unpack=True) #C66 - #for each cross section for each state of the molecule, we have an associated energy - - self.energy_range = np.where(self.energy_neutral<13.6)[0] - #energy_neutral|charged|double_charged are the same - - self.pah_cross_n = ( ( (crossn_1_case1/32)+(crossn_2_case1/36)+(crossn_3_case1/38) +\ - (crossn_1_case2/40)+(crossn_2_case2/42)+(crossn_3_case2/42) +\ - (crossn_1_case3/48)+(crossn_2_case3/54)+(crossn_3_case3/66) )/9 ) * self.n_c - self.pah_cross_c = ( ( (crossc_1_case1/32)+(crossc_2_case1/36)+(crossc_3_case1/38) +\ - (crossc_1_case2/40)+(crossc_2_case2/42)+(crossc_3_case2/42) +\ - (crossc_1_case3/48)+(crossc_2_case3/54)+(crossc_3_case3/66) )/9 ) * self.n_c - self.pah_cross_dc = ( ((crossdc_1_case1/32)+(crossdc_2_case1/36)+(crossdc_3_case1/38) +\ - (crossdc_1_case2/40)+(crossdc_2_case2/42)+(crossdc_3_case2/42) +\ - (crossdc_1_case3/48)+(crossdc_2_case3/54)+(crossdc_3_case3/66) )/9 ) * self.n_c - #cross_n is the average cross section of a pah of all types of size - - ''' Ranges imposed ''' - self.energy_neutral = self.energy_neutral[self.energy_range] - self.energy_charged = self.energy_charged[self.energy_range] - self.energy_double_charged = self.energy_double_charged[self.energy_range] - self.pah_cross_n = self.pah_cross_n[self.energy_range] - self.pah_cross_c = self.pah_cross_c[self.energy_range] - self.pah_cross_dc = self.pah_cross_dc[self.energy_range] - - ''' yield from neutral to the first photoionization ''' - first_part = np.where(self.energy_neutral=self.ip_neutral) & (self.energy_neutral<(self.ip_neutral+9.2) ) )[0] - third_part = np.where((self.energy_neutral>self.ip_neutral+9.2))[0] - - y_1 = np.zeros([len(first_part)]) - y_2 = ( self.energy_neutral[second_part]-self.ip_neutral )/9.2 - y_3 = np.full(len(third_part), 1) - - self.yield_of_first_photoionization = np.concatenate( (y_1,y_2,y_3) ) - - ''' yield from the first photoionization to the second photoionization ''' - alpha = 0.3 #teepness coefficient, see Wenzel et al. 2020 - if (self.n_c >= 32) & (self.n_c < 50): - beta = 0.59 + 8.1e-3 * self.n_c - if self.n_c >= 50: - beta = 1 - - first_part = np.where( self.energy_charged=self.ip_charged ) & ( self.energy_charged<11.3 ) )[0] - third_part = np.where( ( self.energy_charged>=11.3 ) & ( self.energy_charged<12.9 ) )[0] - fourth_part = np.where( ( self.energy_charged>=12.9 ) & ( self.energy_charged<13.6 ) )[0] - - y_1 = np.zeros([len(first_part)]) - y_2 = ( alpha/(11.3-self.ip_charged) ) * (self.energy_charged[second_part]-self.ip_charged) - y_3 = np.full(len(third_part), alpha) - y_4 = ( (beta-alpha)/2.1 ) * (self.energy_charged[fourth_part]-12.9) + alpha - - self.yield_of_second_photoionization = np.concatenate( (y_1,y_2,y_3,y_4) ) - - ''' adaptating the cross section from 0 to 13.6eV''' - pah_crossn = np.interp(self.energy,self.energy_neutral,self.pah_cross_n)*mb #in cm2/Carbon (from Mb/C to cm2/C) - pah_crossc = np.interp(self.energy,self.energy_charged,self.pah_cross_c)*mb #in cm2/Carbon (from Mb/C to cm2/C) - pah_crossdc = np.interp(self.energy,self.energy_double_charged,self.pah_cross_dc)*mb #in cm2/Carbon (from Mb/C to cm2/C) - - yield_n = np.interp(self.energy,self.energy_neutral,self.yield_of_first_photoionization) - yield_c = np.interp(self.energy,self.energy_charged,self.yield_of_second_photoionization) - #interpolation - - '''===================== dust and gas heating calculation ===================''' - energy_range_power_absorbed = np.where(self.energy<=13.6)[0] - energy_range_neutral = np.where(np.logical_and(self.energy >= self.ip_neutral, self.energy <= 13.6))[0] - energy_range_charged = np.where(np.logical_and(self.energy >= self.ip_charged, self.energy <= 13.6))[0] - - ''' photoabsorption of the neutrals, cations and dications molecules ''' - 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 - photo_absorption_c = self.energy_intensity*pah_crossc*2*np.pi #erg s-1 eV-1 - photo_absorption_dc = self.energy_intensity*pah_crossdc*2*np.pi #erg s-1 eV-1 - - ''' power density absorbed for ionization ''' - ionization_absorption_n = yield_n*photo_absorption_n #erg s-1 eV-1 - ionization_absorption_c = yield_c*photo_absorption_c #erg s-1 eV-1 - - ''' number of ionizations ''' - number_ionization_absorption_n = ionization_absorption_n/(self.energy*ev_to_erg) - number_ionization_absorption_c = ionization_absorption_c/(self.energy*ev_to_erg) - #number of ionizations per s per eV for a charge state - - ''' photoemission rate ''' - kpe_neutral = np.trapz(number_ionization_absorption_n[energy_range_neutral], (self.energy-self.ip_neutral)[energy_range_neutral]) - kpe_charged = np.trapz(number_ionization_absorption_c[energy_range_charged], (self.energy-self.ip_charged)[energy_range_charged]) - #in s-1 - - ''' recombination rate ''' - phi = ( 1.85*1e5 )/( self.t_gas*np.sqrt(self.n_c) ) #dimensionless - krec_neutral = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi ) - krec_charged = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi * (1 + z_1) ) - #in s-1 - - ''' population fraction computation ''' - #cations - self.frac_charged = 1/( 1+( krec_neutral/kpe_neutral) + (kpe_charged/krec_charged) ) - #neutrals - self.frac_neutral = (1 - self.frac_charged * (kpe_charged/krec_charged) )/( 1 + (kpe_neutral/krec_neutral) ) - #dications - self.frac_double_charged = (1 - self.frac_neutral)/( 1 + (krec_charged/kpe_charged) ) - - ''' selection of the partition coefficient ''' - partition_coeff = 0.46 #PAH parameter, 0.46 + or - 0.06 - - ''' spectrum of the gas heating per charge state ''' - neutral_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_neutral) *\ - number_ionization_absorption_n * ev_to_erg #erg s-1 eV-1 - charged_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_charged) *\ - number_ionization_absorption_c * ev_to_erg #erg s-1 eV-1 - #partition_coeff * (E-IP) is the kinetic energy of the photoelectron following absorption of a UV photon of energy E - - ''' gas heating rate per charge state ''' - neutral_gas_heating_rate = np.trapz(neutral_heating_rate_spectrum[energy_range_neutral],\ - (self.energy-self.ip_neutral)[energy_range_neutral] ) #erg s-1 molecule-1 - charged_gas_heating_rate = np.trapz(charged_heating_rate_spectrum[energy_range_charged],\ - (self.energy-self.ip_charged)[energy_range_charged] ) #erg s-1 molecule-1 - #powers injected in the gas by photoelectrons ejected from - #neutral and cationic PAHs - - '''======================== heating efficiencies ========================''' - - ''' total power injected into the gas via the photoelectrons from pahs ''' - total_injected_power = self.frac_neutral * neutral_gas_heating_rate +\ - self.frac_charged * charged_gas_heating_rate - #erg s-1 /(molecule of size n_c) - - ''' heating rate of the molecule itself ''' - heating_pah_n = np.trapz(photo_absorption_n[energy_range_power_absorbed],\ - self.energy[energy_range_power_absorbed]) #erg s-1 - heating_pah_c = np.trapz(photo_absorption_c[energy_range_power_absorbed],\ - self.energy[energy_range_power_absorbed]) #erg s-1 - heating_pah_dc = np.trapz(photo_absorption_dc[energy_range_power_absorbed],\ - self.energy[energy_range_power_absorbed]) #erg s-1 - #power absorbed by each PAH charge state - - ''' total power of the radiation absorbed by PAHs ''' - total_absorbed_radiation_power = self.frac_neutral * heating_pah_n +\ - self.frac_charged * heating_pah_c +\ - self.frac_double_charged * heating_pah_dc - #erg s-1 /(molecule of size n_c) - - ''' heating efficiency ''' - self.heating_efficiency = total_injected_power/total_absorbed_radiation_power - - ''' gas heating ''' - self.total_gas_heating = total_injected_power * (self.fc_pah/self.n_c) * 2.7e-4 - #2.7e-4 : elemental abundance of C relative to H (Tielens 2021) - - return self.total_gas_heating, self.g_0 , self.gamma , self.t_gas, self.n_e, self.n_c \ No newline at end of file diff --git a/three_levels_model.py b/three_levels_model.py new file mode 100644 index 0000000..20840ff --- /dev/null +++ b/three_levels_model.py @@ -0,0 +1,343 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jun 14 10:48:49 2022 + +@author: frede +""" +import numpy as np +import radiation_fields as rf + +''' + +This programme seeks to calculate how a neutral gas subjected to a radiation field +is heated, involving pahs in the ionization of the gas + +This program works with the script "radiation_fields.py" +which calculates the scale factor of the radiation field G0 + +Example of execution of the program for the star of 10 solar radius HD200775, +for a distance to the star at 20pc, for a gas at 750K, an electron density +at 2.4 cm-3 and a size of PAHs à 54 C atoms : + +HeatingGas(filename='HD200775_RF.txt', star_radius=10, t_gas=750, + n_e=1.6, n_c=54, parsec=20) + +If i consider the Interstellar Radiation Field with +an approach described by Habing (1968), for a gas at 50K and a +fraction of cosmic carbon locked in PAHs at 5%: + +HeatingGas('habing1968.txt', 1, t_gas=100, n_e=1.6, n_c=54, + parsec=1, fc_pah=0.05, ISRF=True) + +While considering ISRF, the value of the parameters +parsec and star_radius are no longer important + +''' + +''' Constants ''' +h = 6.62607015e-34 #Planck constant in J s +c = 299792458 #Light speed in m s-1 +eps_0 = 8.85418782e-12 #epsilon_0, vacuum permitivity in F (Farad) m-1 +z_0 = 0 #the charge state of the neutral molecule +z_1 = 1 #the charge state of the ionized 1 molecule +one_in_4_pi_eps_0 = 1/( 4*np.pi*(eps_0/1e9) ) + +''' Conversions ''' +ev = 1.602176634e-19 # 1 ev = 1.60218e-19 J and value of electron charge +erg = 1e-7 # 1 erg = 1e-7 J +ev_to_erg = ev/erg #1 eV in erg +mb = 1e-18 #1Mb = 1e-18cm2, Mb for Megabarn (unit used to express the cross sectional area of nuclei) + +''' Saving parameters ''' +dust_heating_rate = np.zeros([3]) +ionization_rate = np.zeros([2]) +recombination_rate = np.zeros([2]) +gas_heating_rate = np.zeros([2]) +intrinsic_efficiency = np.zeros([2]) + +class HeatingGas: + + """ + + ---------- + + Returns total_gas_heating, g_0 , gamma , t_gas, n_e , n_c + + total_gas_heating : float, + gas heating rate + g_0 : float, + scaling factor of the radiation field + gamma : float, + ( g_0 * sqrt(t_gas) )/n_e ionization parameter + t_gas : float, + gas temperature + n_e : float, + electron density in cm-3 + n_c : float, + number of carbon atoms in pah molecules (size of the pah) + + ------- + + """ + + def __init__(self, filename, star_radius, t_gas, n_e, n_c, parsec, fc_pah=0.1, ISRF=False): + + """ + + ---------- + filename : str, + name of the file containing wavelength in nm and + intensity in erg cm-2 s-1 sr-1 nm-1 of a star or of the interstellar medium + star_radius : float, + star radius, in unit of solar radius + n_e : float, + electron density in cm-3 (n_e = n_h * 1.6e-4, n_h : hydrogene density) + parsec : float, + distance in parsec from the star + fc_pah : float, + fraction of cosmic carbon locked in PAHs (default: 0.1) + ISRF : bool, + Interstellar Radiation Field; if true, then the filename is a file for ISRF + if false, the radiation field of a star is studied (default: False) + ------- + + """ + + ''' input parameters ''' + self.filename = filename + self.t_gas = t_gas + self.n_e = n_e + self.n_c = n_c + self.parsec = parsec + self.star_radius = star_radius + self.fc_pah = fc_pah + self.ISRF = ISRF + + ''' parameters to be observed ''' + self.g_0 = None + self.distance = None + self.wavelength = None + self.wavelength_intensity = None + self.energy_intensity = None + self.energy = None + self.energy_range = None + self.energy_neutral = None + self.energy_charged = None + self.energy_double_charged = None + self.pah_cross_n = None + self.pah_cross_c = None + self.pah_cross_dc = None + self.ip_neutral = None + self.ip_charged = None + self.yield_of_first_photoionization = None + self.yield_of_second_photoionization = None + self.heating_efficiency = None + self.total_gas_heating = None + self.frac_neutral = None + self.frac_charged = None + self.frac_double_charged = None + + def parameters(self): + + ''' others parameters to be returned ''' + 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) + self.gamma = ( self.g_0 * np.sqrt(self.t_gas) ) / self.n_e #ionization parameter + + ''' Ionization Potential (IP) estimation ''' + a = (self.n_c/468)**(1/3) #molecule diameter, in nm + 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) + #IP to ionize the molecule : neutral to charged, in ev + 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) + #IP to ionize the charged molecule : charged 1 to charged 2, in ev + + '''==========================|building of the cross section|=======================''' + ''' derives a mean photoabsorption cross section of the molecule considered, in 3 size ranges''' + + ''' small size ''' + self.energy_neutral,crossn_1_case1 = np.loadtxt('./neutrals/ovalene_neutral.txt',unpack=True) #C32 + self.energy_charged,crossc_1_case1 = np.loadtxt('./cations/ovalene_cation.txt',unpack=True) #C32 + self.energy_double_charged,crossdc_1_case1 = np.loadtxt('./dications/ovalene_dication.txt',unpack=True) #C32 + + self.energy_neutral,crossn_2_case1 = np.loadtxt('./neutrals/tetrabenzocoronene_neutral.txt',unpack=True) #C36 + self.energy_charged,crossc_2_case1 = np.loadtxt('./cations/tetrabenzocoronene_cation.txt',unpack=True) #C36 + self.energy_double_charged,crossdc_2_case1 = np.loadtxt('./dications/tetrabenzocoronene_dication.txt',unpack=True) #C36 + + self.energy_neutral,crossn_3_case1 = np.loadtxt('./neutrals/circumbiphenyl_neutral.txt',unpack=True) #C38 + self.energy_charged,crossc_3_case1 = np.loadtxt('./cations/circumbiphenyl_cation.txt',unpack=True) #C38 + self.energy_double_charged,crossdc_3_case1 = np.loadtxt('./dications/circumbiphenyl_dication.txt',unpack=True) #C38 + + ''' medium size ''' + self.energy_neutral,crossn_1_case2=np.loadtxt('./neutrals/circumanthracene_neutral.txt',unpack=True) #C40 + self.energy_charged,crossc_1_case2=np.loadtxt('./cations/circumanthracene_cation.txt',unpack=True) #C40 + self.energy_double_charged,crossdc_1_case2=np.loadtxt('./dications/circumanthracene_dication.txt',unpack=True) #C40 + + self.energy_neutral,crossn_2_case2 = np.loadtxt('./neutrals/circumpyrene_neutral.txt',unpack=True) #C42 + self.energy_charged,crossc_2_case2 = np.loadtxt('./cations/circumpyrene_cation.txt',unpack=True) #C42 + self.energy_double_charged,crossdc_2_case2 = np.loadtxt('./dications/circumpyrene_dication.txt',unpack=True) #C42 + + self.energy_neutral,crossn_3_case2 = np.loadtxt('./neutrals/hexabenzocoronene_neutral.txt',unpack=True) #C42 + self.energy_charged,crossc_3_case2 = np.loadtxt('./cations/hexabenzocoronene_cation.txt',unpack=True) #C42 + self.energy_double_charged,crossdc_3_case2 = np.loadtxt('./dications/hexabenzocoronene_dication.txt',unpack=True) #C42 + + ''' large size ''' + self.energy_neutral,crossn_1_case3 = np.loadtxt('./neutrals/dicoronylene_neutral.txt',unpack=True) #C48 + self.energy_charged,crossc_1_case3 = np.loadtxt('./cations/dicoronylene_cation.txt',unpack=True) #C48 + self.energy_double_charged,crossdc_1_case3 = np.loadtxt('./dications/dicoronylene_dication.txt',unpack=True) #C48 + + self.energy_neutral,crossn_2_case3 = np.loadtxt('./neutrals/circumcoronene_neutral.txt',unpack=True) #C54 + self.energy_charged,crossc_2_case3 = np.loadtxt('./cations/circumcoronene_cation.txt',unpack=True) #C54 + self.energy_double_charged,crossdc_2_case3 = np.loadtxt('./dications/circumcoronene_dication.txt',unpack=True) #C54 + + self.energy_neutral,crossn_3_case3 = np.loadtxt('./neutrals/circumovalene_neutral.txt',unpack=True) #C66 + self.energy_charged,crossc_3_case3 = np.loadtxt('./cations/circumovalene_cation.txt',unpack=True) #C66 + self.energy_double_charged,crossdc_3_case3 = np.loadtxt('./dications/circumovalene_dication.txt',unpack=True) #C66 + #for each cross section for each state of the molecule, we have an associated energy + + self.energy_range = np.where(self.energy_neutral<13.6)[0] + #energy_neutral|charged|double_charged are the same + + self.pah_cross_n = ( ( (crossn_1_case1/32)+(crossn_2_case1/36)+(crossn_3_case1/38) +\ + (crossn_1_case2/40)+(crossn_2_case2/42)+(crossn_3_case2/42) +\ + (crossn_1_case3/48)+(crossn_2_case3/54)+(crossn_3_case3/66) )/9 ) * self.n_c + self.pah_cross_c = ( ( (crossc_1_case1/32)+(crossc_2_case1/36)+(crossc_3_case1/38) +\ + (crossc_1_case2/40)+(crossc_2_case2/42)+(crossc_3_case2/42) +\ + (crossc_1_case3/48)+(crossc_2_case3/54)+(crossc_3_case3/66) )/9 ) * self.n_c + self.pah_cross_dc = ( ((crossdc_1_case1/32)+(crossdc_2_case1/36)+(crossdc_3_case1/38) +\ + (crossdc_1_case2/40)+(crossdc_2_case2/42)+(crossdc_3_case2/42) +\ + (crossdc_1_case3/48)+(crossdc_2_case3/54)+(crossdc_3_case3/66) )/9 ) * self.n_c + #cross_n is the average cross section of a pah of all types of size + + ''' Ranges imposed ''' + self.energy_neutral = self.energy_neutral[self.energy_range] + self.energy_charged = self.energy_charged[self.energy_range] + self.energy_double_charged = self.energy_double_charged[self.energy_range] + self.pah_cross_n = self.pah_cross_n[self.energy_range] + self.pah_cross_c = self.pah_cross_c[self.energy_range] + self.pah_cross_dc = self.pah_cross_dc[self.energy_range] + + ''' yield from neutral to the first photoionization ''' + first_part = np.where(self.energy_neutral=self.ip_neutral) & (self.energy_neutral<(self.ip_neutral+9.2) ) )[0] + third_part = np.where((self.energy_neutral>self.ip_neutral+9.2))[0] + + y_1 = np.zeros([len(first_part)]) + y_2 = ( self.energy_neutral[second_part]-self.ip_neutral )/9.2 + y_3 = np.full(len(third_part), 1) + + self.yield_of_first_photoionization = np.concatenate( (y_1,y_2,y_3) ) + + ''' yield from the first photoionization to the second photoionization ''' + alpha = 0.3 #teepness coefficient, see Wenzel et al. 2020 + if (self.n_c >= 32) & (self.n_c < 50): + beta = 0.59 + 8.1e-3 * self.n_c + if self.n_c >= 50: + beta = 1 + + first_part = np.where( self.energy_charged=self.ip_charged ) & ( self.energy_charged<11.3 ) )[0] + third_part = np.where( ( self.energy_charged>=11.3 ) & ( self.energy_charged<12.9 ) )[0] + fourth_part = np.where( ( self.energy_charged>=12.9 ) & ( self.energy_charged<13.6 ) )[0] + + y_1 = np.zeros([len(first_part)]) + y_2 = ( alpha/(11.3-self.ip_charged) ) * (self.energy_charged[second_part]-self.ip_charged) + y_3 = np.full(len(third_part), alpha) + y_4 = ( (beta-alpha)/2.1 ) * (self.energy_charged[fourth_part]-12.9) + alpha + + self.yield_of_second_photoionization = np.concatenate( (y_1,y_2,y_3,y_4) ) + + ''' adaptating the cross section from 0 to 13.6eV''' + pah_crossn = np.interp(self.energy,self.energy_neutral,self.pah_cross_n)*mb #in cm2/Carbon (from Mb/C to cm2/C) + pah_crossc = np.interp(self.energy,self.energy_charged,self.pah_cross_c)*mb #in cm2/Carbon (from Mb/C to cm2/C) + pah_crossdc = np.interp(self.energy,self.energy_double_charged,self.pah_cross_dc)*mb #in cm2/Carbon (from Mb/C to cm2/C) + + yield_n = np.interp(self.energy,self.energy_neutral,self.yield_of_first_photoionization) + yield_c = np.interp(self.energy,self.energy_charged,self.yield_of_second_photoionization) + #interpolation + + '''===================== dust and gas heating calculation ===================''' + energy_range_power_absorbed = np.where(self.energy<=13.6)[0] + energy_range_neutral = np.where(np.logical_and(self.energy >= self.ip_neutral, self.energy <= 13.6))[0] + energy_range_charged = np.where(np.logical_and(self.energy >= self.ip_charged, self.energy <= 13.6))[0] + + ''' photoabsorption of the neutrals, cations and dications molecules ''' + 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 + photo_absorption_c = self.energy_intensity*pah_crossc*2*np.pi #erg s-1 eV-1 + photo_absorption_dc = self.energy_intensity*pah_crossdc*2*np.pi #erg s-1 eV-1 + + ''' power density absorbed for ionization ''' + ionization_absorption_n = yield_n*photo_absorption_n #erg s-1 eV-1 + ionization_absorption_c = yield_c*photo_absorption_c #erg s-1 eV-1 + + ''' number of ionizations ''' + number_ionization_absorption_n = ionization_absorption_n/(self.energy*ev_to_erg) + number_ionization_absorption_c = ionization_absorption_c/(self.energy*ev_to_erg) + #number of ionizations per s per eV for a charge state + + ''' photoemission rate ''' + kpe_neutral = np.trapz(number_ionization_absorption_n[energy_range_neutral], (self.energy-self.ip_neutral)[energy_range_neutral]) + kpe_charged = np.trapz(number_ionization_absorption_c[energy_range_charged], (self.energy-self.ip_charged)[energy_range_charged]) + #in s-1 + + ''' recombination rate ''' + phi = ( 1.85*1e5 )/( self.t_gas*np.sqrt(self.n_c) ) #dimensionless + krec_neutral = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi ) + krec_charged = self.n_e * 1.28e-10 * self.n_c * np.sqrt(self.t_gas) * ( 1 + phi * (1 + z_1) ) + #in s-1 + + ''' population fraction computation ''' + #cations + self.frac_charged = 1/( 1+( krec_neutral/kpe_neutral) + (kpe_charged/krec_charged) ) + #neutrals + self.frac_neutral = (1 - self.frac_charged * (kpe_charged/krec_charged) )/( 1 + (kpe_neutral/krec_neutral) ) + #dications + self.frac_double_charged = (1 - self.frac_neutral)/( 1 + (krec_charged/kpe_charged) ) + + ''' selection of the partition coefficient ''' + partition_coeff = 0.46 #PAH parameter, 0.46 + or - 0.06 + + ''' spectrum of the gas heating per charge state ''' + neutral_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_neutral) *\ + number_ionization_absorption_n * ev_to_erg #erg s-1 eV-1 + charged_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_charged) *\ + number_ionization_absorption_c * ev_to_erg #erg s-1 eV-1 + #partition_coeff * (E-IP) is the kinetic energy of the photoelectron following absorption of a UV photon of energy E + + ''' gas heating rate per charge state ''' + neutral_gas_heating_rate = np.trapz(neutral_heating_rate_spectrum[energy_range_neutral],\ + (self.energy-self.ip_neutral)[energy_range_neutral] ) #erg s-1 molecule-1 + charged_gas_heating_rate = np.trapz(charged_heating_rate_spectrum[energy_range_charged],\ + (self.energy-self.ip_charged)[energy_range_charged] ) #erg s-1 molecule-1 + #powers injected in the gas by photoelectrons ejected from + #neutral and cationic PAHs + + '''======================== heating efficiencies ========================''' + + ''' total power injected into the gas via the photoelectrons from pahs ''' + total_injected_power = self.frac_neutral * neutral_gas_heating_rate +\ + self.frac_charged * charged_gas_heating_rate + #erg s-1 /(molecule of size n_c) + + ''' heating rate of the molecule itself ''' + heating_pah_n = np.trapz(photo_absorption_n[energy_range_power_absorbed],\ + self.energy[energy_range_power_absorbed]) #erg s-1 + heating_pah_c = np.trapz(photo_absorption_c[energy_range_power_absorbed],\ + self.energy[energy_range_power_absorbed]) #erg s-1 + heating_pah_dc = np.trapz(photo_absorption_dc[energy_range_power_absorbed],\ + self.energy[energy_range_power_absorbed]) #erg s-1 + #power absorbed by each PAH charge state + + ''' total power of the radiation absorbed by PAHs ''' + total_absorbed_radiation_power = self.frac_neutral * heating_pah_n +\ + self.frac_charged * heating_pah_c +\ + self.frac_double_charged * heating_pah_dc + #erg s-1 /(molecule of size n_c) + + ''' heating efficiency ''' + self.heating_efficiency = total_injected_power/total_absorbed_radiation_power + + ''' gas heating ''' + self.total_gas_heating = total_injected_power * (self.fc_pah/self.n_c) * 2.7e-4 + #2.7e-4 : elemental abundance of C relative to H (Tielens 2021) + + return self.total_gas_heating, self.g_0 , self.gamma , self.t_gas, self.n_e, self.n_c \ No newline at end of file -- libgit2 0.21.2