# -*- 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()