plot_program.py 10 KB
# -*- 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()