addition of the ionization potential associated with the anion

Jalabert
1 parent 53761528
Showing 1 changed file with 13 additions and 11 deletions Show diff stats
four_levels_model.py
@@ -9,8 +9,7 @@ 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 programme computes the photoelectric heating rate and efficiency by PAHs
 This program works with the script "radiation_fields.py" 
 which calculates the scale factor of the radiation field G0
@@ -61,7 +60,7 @@ class HeatingGas:
     ----------
-    Returns total_gas_heating, heating efficiency, g_0 , gamma , t_gas, n_e , n_c
+    Returns total_gas_heating, heating_efficiency, g_0 , gamma , t_gas, n_e , n_c
     total_gas_heating : float,
         gas heating rate in erg s-1 H-1
@@ -131,6 +130,7 @@ class HeatingGas:
         self.pah_cross_n = None
         self.pah_cross_c = None
         self.pah_cross_dc = None
+        self.ip_negative_charged = None
         self.ip_neutral = None
         self.ip_charged = None
         self.detachment_yield = None
@@ -156,6 +156,8 @@ class HeatingGas:
         #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
+        self.ip_negative_charged = 6
+        #IP to detache the electron from a PAH anion, in ev
         '''==========================|building of the cross section|======================='''
         ''' derives a mean photoabsorption cross section of the molecule considered, in 3 size ranges'''
@@ -239,7 +241,7 @@ class HeatingGas:
         self.pah_cross_dc = self.pah_cross_dc[self.energy_range]
         ''' yield from anion to neutral '''
-        part = np.where(self.energy_negative_charged <=13.6)[0]
+        part = np.where(self.energy_negative_charged<=13.6)[0]
         self.detachment_yield = np.full(len(part), 1)
         #the anion being unstable, any energy is enough to detach the electron
@@ -288,10 +290,10 @@ class HeatingGas:
         '''===================== dust and gas heating calculation ==================='''
         energy_range_power_absorbed = np.where(self.energy<=13.6)[0]
-        energy_range_anion = np.where(self.energy<=13.6)[0]
+        energy_range_anion = np.where(np.logical_and(self.energy >= self.ip_negative_charged, 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_a = self.energy_intensity*pah_crossa*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_n = self.energy_intensity*pah_crossn*2*np.pi #erg s-1 eV-1
@@ -310,7 +312,7 @@ class HeatingGas:
         #number of ionizations per s per eV for a charge state
         ''' detachment rate '''
-        kdet = np.trapz(number_ionization_absorption_a[energy_range_anion], self.energy[energy_range_anion])
+        kdet = np.trapz(number_ionization_absorption_a[energy_range_anion], (self.energy-self.ip_negative_charged)[energy_range_anion])
         #in s-1
         ''' photoemission rate '''
@@ -323,7 +325,7 @@ class HeatingGas:
         #in s-1
         ''' recombination rate '''
-        phi = ( 1.85*1e5 )/( self.t_gas*np.sqrt(self.n_c) ) #dimensionless
+        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
@@ -343,7 +345,7 @@ class HeatingGas:
         ''' spectrum of the gas heating per charge state '''
-        anion_heating_rate_spectrum = 1 * (self.energy-0) *\
+        anion_heating_rate_spectrum = 1 * (self.energy-self.ip_negative_charged) *\
                                       number_ionization_absorption_a * ev_to_erg #erg s-1 eV-1
         neutral_heating_rate_spectrum = partition_coeff * (self.energy-self.ip_neutral) *\
                                         number_ionization_absorption_n * ev_to_erg #erg s-1 eV-1
@@ -353,7 +355,7 @@ class HeatingGas:
         ''' gas heating rate per charge state '''
         anion_gas_heating_rate =   np.trapz(anion_heating_rate_spectrum[energy_range_anion],\
-                                            self.energy[energy_range_anion] ) #erg s-1 molecule-1
+                                            (self.energy-self.ip_negative_charged)[energy_range_anion] ) #erg s-1 molecule-1
         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],\
@@ -394,5 +396,5 @@ class HeatingGas:
         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)  
         #total gas heating in erg s-1 H-1  
-        
+                
         return self.total_gas_heating, self.heating_efficiency, self.g_0 , self.gamma, self.t_gas, self.n_e, self.n_c
 \ No newline at end of file