# -*- coding: utf-8 -*-
"""
Created on Fri Jul 15 16:47:35 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_negative_charged = None
        self.energy_neutral = None
        self.energy_charged = None
        self.energy_double_charged = None
        self.pah_cross_a = 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.detachment_yield = None
        self.yield_of_first_photoionization = None
        self.yield_of_second_photoionization = None
        self.heating_efficiency = None
        self.total_gas_heating = None
        self.frac_anion = 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_negative_charged,cross_anion = np.loadtxt('./anions/coronene_anion.txt',unpack=True) #C24
        
        self.energy_negative_charged,crossa_1_case1 = np.loadtxt('./anions/ovalene_anion.txt',unpack=True) #C32
        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_negative_charged,crossa_2_case1 = np.loadtxt('./anions/tetrabenzocoronene_anion.txt',unpack=True) #C36
        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_negative_charged,crossa_3_case1 = np.loadtxt('./anions/circumbiphenyl_anion.txt',unpack=True) #C38
        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_negative_charged,crossa_1_case2 = np.loadtxt('./anions/circumanthracene_anion.txt',unpack=True) #C40
        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_negative_charged,crossa_2_case2 = np.loadtxt('./anions/circumpyrene_anion.txt',unpack=True) #C42
        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_negative_charged,crossa_3_case2 = np.loadtxt('./anions/hexabenzocoronene_anion.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_negative_charged,crossa_1_case3 = np.loadtxt('./anions/dicoronylene_anion.txt',unpack=True) #C48
        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_negative_charged,crossa_2_case3 = np.loadtxt('./anions/circumcoronene_anion.txt',unpack=True) #C54
        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_negative_charged,crossa_3_case3 = np.loadtxt('./anions/circumovalene_anion.txt',unpack=True) #C66
        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]
        
        self.pah_cross_a = ( ( (cross_anion/24)   +(crossa_1_case1/32)+(crossa_2_case1/36) +\
                               (crossa_3_case1/38)+(crossa_1_case2/40)+(crossa_2_case2/42) +\
                               (crossa_3_case2/42)+(crossa_1_case3/48)+(crossa_2_case3/54) +\
                               (crossa_3_case3/66)                                         )/10 ) * self.n_c
        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_negative_charged = self.energy_negative_charged[self.energy_range]
        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_a = self.pah_cross_a[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 anion to neutral '''
        part = np.where(self.energy_negative_charged >= 0)[0]

        self.detachment_yield = np.full(len(part), 1)
        #the anion being unstable, any energy is enough to detach the electron

        ''' yield from neutral to the first photoionization '''
        first_part = np.where(self.energy_neutral<self.ip_neutral)[0]

        second_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 )[0]
        second_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_crossa = np.interp(self.energy,self.energy_negative_charged,self.pah_cross_a)*mb #in cm2/Carbon (from Mb/C to cm2/C)
        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_a = np.interp(self.energy,self.energy_negative_charged,self.detachment_yield)
        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_anion = 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_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
        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_a = yield_a*photo_absorption_a #erg s-1 eV-1
        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_a = ionization_absorption_a/(self.energy*ev_to_erg)
        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
    
        ''' detachment rate '''
        kdet = np.trapz(number_ionization_absorption_a[energy_range_anion], self.energy[energy_range_anion])
        #in s-1
    
        ''' 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
        
        ''' attachment rate ''' #for coronene
        katt = self.n_e * 2.74e-9 * (self.t_gas/300)**(0.11) * np.exp(-(-1.12)/self.t_gas)
        #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 '''
        #anions
        self.frac_anion = 1/(1 + (kdet/katt) + (kdet*kpe_neutral)/(katt*krec_neutral) + (kdet*kpe_neutral*kpe_charged)/(katt*krec_neutral*krec_charged))
        #neutrals
        self.frac_neutral = 1/(1 + (katt/kdet) + (kpe_neutral/krec_neutral) + (kpe_neutral*kpe_charged)/(krec_neutral*krec_charged))
        #cations
        self.frac_charged = 1/(1 + (krec_neutral/kpe_neutral) + (krec_neutral*katt)/(kdet*kpe_neutral) + (kpe_charged/krec_charged))
        #dications
        self.frac_double_charged = 1/(1 + (krec_charged/kpe_charged) + (krec_neutral*krec_charged)/(kpe_neutral*kpe_charged) + (katt*krec_neutral*krec_charged)/(kdet*kpe_neutral*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 '''

        anion_heating_rate_spectrum = 1 * (self.energy-0) *\
                                      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
        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 '''
        anion_gas_heating_rate =   np.trapz(anion_heating_rate_spectrum[energy_range_anion],\
                                            self.energy[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],\
                                            (self.energy-self.ip_charged)[energy_range_charged] ) #erg s-1 molecule-1
        #powers injected in the gas by photoelectrons ejected from
        #anionic neutral and cationic PAHs
        
        '''======================== heating efficiencies ========================'''
        
        ''' total power injected into the gas via the photoelectrons from pahs '''
        total_injected_power = self.frac_anion   * anion_gas_heating_rate   +\
                               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_a = np.trapz(photo_absorption_a,\
                                        self.energy) #erg s-1
        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
               
        ''' total power of the radiation absorbed by pahs '''
        total_absorbed_radiation_power = self.frac_anion   * heating_pah_a +\
                                         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